pml_time_harmonic_linear_elasticity_elements.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for general linear elasticity elements
27
28// Include guards to prevent multiple inclusion of the header
29#ifndef OOMPH_PML_TIME_HARMONIC_LINEAR_ELASTICITY_ELEMENTS_HEADER
30#define OOMPH_PML_TIME_HARMONIC_LINEAR_ELASTICITY_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37
38#ifdef OOMPH_HAS_MPI
39#include "mpi.h"
40#endif
41
42#include <complex>
43#include "math.h"
44
45// OOMPH-LIB headers
46#include "../generic/Qelements.h"
47#include "../generic/mesh.h"
48#include "../generic/hermite_elements.h"
50#include "../generic/projection.h"
51#include "../generic/nodes.h"
52#include "../generic/oomph_utilities.h"
53#include "../generic/pml_meshes.h"
54#include "../generic/pml_mapping_functions.h"
55
56// The meshes (needed for building of pml meshes!)
57// Include template files to avoid unnecessary re-compilation
58// (*.template.h files get included indirectly).
59//#include "../meshes/triangle_mesh.template.cc"
60//#include "../meshes/rectangular_quadmesh.template.cc"
61
62// Why not just to include the *.h files, Just as all other files
63// #include "../meshes/triangle_mesh.template.h"
64// #include "../meshes/rectangular_quadmesh.template.h"
65
66namespace oomph
67{
68 //=======================================================================
69 /// A base class for elements that solve the equations of time-harmonic linear
70 /// elasticity in Cartesian coordinates.
71 /// Combines a few generic functions that are shared by
72 /// PMLTimeHarmonicLinearElasticityEquations
73 /// and PMLTimeHarmonicLinearElasticityEquationsWithPressure
74 /// (Note: The latter
75 /// don't exist yet but will be written as soon as somebody needs them...)
76 //=======================================================================
77 template<unsigned DIM>
79 : public virtual PMLElementBase<DIM>,
80 public virtual FiniteElement
81 {
82 public:
83 /// Constructor: Set null pointers for constitutive law and for
84 /// isotropic growth function. Set physical parameter values to
85 /// default values, and set body force to zero.
90 {
93 }
94
95 /// Return the index at which the i-th real or imag unknown
96 /// displacement component is stored. The default value is appropriate for
97 /// single-physics problems:
98 virtual inline std::complex<unsigned> u_index_time_harmonic_linear_elasticity(
99 const unsigned i) const
100 {
101 return std::complex<unsigned>(i, i + DIM);
102 }
103
104 /// Compute vector of FE interpolated displacement u at local coordinate s
106 const Vector<double>& s, Vector<std::complex<double>>& disp) const
107 {
108 // Find number of nodes
109 unsigned n_node = nnode();
110
111 // Local shape function
112 Shape psi(n_node);
113
114 // Find values of shape function
115 shape(s, psi);
116
117 for (unsigned i = 0; i < DIM; i++)
118 {
119 // Index at which the nodal value is stored
120 std::complex<unsigned> u_nodal_index =
122
123 // Initialise value of u
124 disp[i] = std::complex<double>(0.0, 0.0);
125
126 // Loop over the local nodes and sum
127 for (unsigned l = 0; l < n_node; l++)
128 {
129 const std::complex<double> u_value(
130 nodal_value(l, u_nodal_index.real()),
131 nodal_value(l, u_nodal_index.imag()));
132
133 disp[i] += u_value * psi[l];
134 }
135 }
136 }
137
138 /// Return FE interpolated displacement u[i] at local coordinate s
140 const Vector<double>& s, const unsigned& i) const
141 {
142 // Find number of nodes
143 unsigned n_node = nnode();
144
145 // Local shape function
146 Shape psi(n_node);
147
148 // Find values of shape function
149 shape(s, psi);
150
151 // Get nodal index at which i-th velocity is stored
152 std::complex<unsigned> u_nodal_index =
154
155 // Initialise value of u
156 std::complex<double> interpolated_u(0.0, 0.0);
157
158 // Loop over the local nodes and sum
159 for (unsigned l = 0; l < n_node; l++)
160 {
161 const std::complex<double> u_value(
162 nodal_value(l, u_nodal_index.real()),
163 nodal_value(l, u_nodal_index.imag()));
164
165 interpolated_u += u_value * psi[l];
166 }
167
168 return (interpolated_u);
169 }
170
171
172 /// Function pointer to function that specifies the body force
173 /// as a function of the Cartesian coordinates FCT(x,b) --
174 /// x and b are Vectors!
175 typedef void (*BodyForceFctPt)(const Vector<double>& x,
177
178 /// Return the pointer to the elasticity_tensor
180 {
182 }
183
184 /// Access function to the entries in the elasticity tensor
185 inline std::complex<double> E(const unsigned& i,
186 const unsigned& j,
187 const unsigned& k,
188 const unsigned& l) const
189 {
190 return (*Elasticity_tensor_pt)(i, j, k, l);
191 }
192
193 /// Access function to nu in the elasticity tensor
194 inline double nu() const
195 {
196 return Elasticity_tensor_pt->nu();
197 }
198
199 /// Access function for square of non-dim frequency
200 const double& omega_sq() const
201 {
202 return *Omega_sq_pt;
203 }
204
205 /// Access function for square of non-dim frequency
206 double*& omega_sq_pt()
207 {
208 return Omega_sq_pt;
209 }
210
211 /// Access function: Pointer to body force function
213 {
214 return Body_force_fct_pt;
215 }
216
217 /// Access function: Pointer to body force function (const version)
219 {
220 return Body_force_fct_pt;
221 }
222
223 /// Return the Cauchy stress tensor, as calculated
224 /// from the elasticity tensor at specified local coordinate
225 /// Virtual so separaete versions can (and must!) be provided
226 /// for displacement and pressure-displacement formulations.
227 virtual void get_stress(const Vector<double>& s,
228 DenseMatrix<std::complex<double>>& sigma) const = 0;
229
230 /// Return the strain tensor
231 void get_strain(const Vector<double>& s,
232 DenseMatrix<std::complex<double>>& strain) const;
233
234 /// Evaluate body force at Eulerian coordinate x
235 /// (returns zero vector if no body force function pointer has been set)
236 inline void body_force(const Vector<double>& x,
237 Vector<std::complex<double>>& b) const
238 {
239 // If no function has been set, return zero vector
240 if (Body_force_fct_pt == 0)
241 {
242 // Get spatial dimension of element
243 unsigned n = dim();
244 for (unsigned i = 0; i < n; i++)
245 {
246 b[i] = std::complex<double>(0.0, 0.0);
247 }
248 }
249 else
250 {
251 // Get body force
252 (*Body_force_fct_pt)(x, b);
253 }
254 }
255
256
257 /// Pure virtual function in which we specify the
258 /// values to be pinned (and set to zero) on the outer edge of
259 /// the pml layer. All of them! Vector is resized internally.
261 Vector<unsigned>& values_to_pin)
262 {
263 values_to_pin.resize(DIM * 2);
264 for (unsigned j = 0; j < DIM * 2; j++)
265 {
266 values_to_pin[j] = j;
267 }
268 }
269
270 /// The number of "DOF types" that degrees of freedom in this element
271 /// are sub-divided into: for now lump them all into one DOF type.
272 /// Can be adjusted later
273 unsigned ndof_types() const
274 {
275 return 1;
276 }
277
278 /// Create a list of pairs for all unknowns in this element,
279 /// so that the first entry in each pair contains the global equation
280 /// number of the unknown, while the second one contains the number
281 /// of the "DOF types" that this unknown is associated with.
282 /// (Function can obviously only be called if the equation numbering
283 /// scheme has been set up.)
285 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
286 {
287 // temporary pair (used to store dof lookup prior to being added
288 // to list)
289 std::pair<unsigned long, unsigned> dof_lookup;
290
291 // number of nodes
292 const unsigned n_node = this->nnode();
293
294 // Integer storage for local unknown
295 int local_unknown = 0;
296
297 // Loop over the nodes
298 for (unsigned n = 0; n < n_node; n++)
299 {
300 // Loop over dimension (real and imag)
301 for (unsigned i = 0; i < 2 * DIM; i++)
302 {
303 // If the variable is free
304 local_unknown = nodal_local_eqn(n, i);
305
306 // ignore pinned values
307 if (local_unknown >= 0)
308 {
309 // store dof lookup in temporary pair: First entry in pair
310 // is global equation number; second entry is dof type
311 dof_lookup.first = this->eqn_number(local_unknown);
312 dof_lookup.second = 0;
313
314 // add to list
315 dof_lookup_list.push_front(dof_lookup);
316 }
317 }
318 }
319 }
320
321 /// Compute pml coefficients at position x and integration point ipt.
322 /// pml_inverse_jacobian_diagonal contains the diagonal terms from the
323 /// inverse of the Jacobian of the PML transformation. These are used to
324 /// transform derivatives in real x to derivatives in transformed space
325 /// \f$\tilde x \f$. This can be interpreted as an anisotropic stiffness.
326 /// pml_jacobian_det is the determinant of the Jacobian of the PML
327 /// transformation, this allows us to transform volume integrals in
328 /// transformed space to real space.
329 /// This can be interpreted as a mass factor
330 /// If the PML is not enabled via enable_pml, both default to 1.0.
332 const unsigned& ipt,
333 const Vector<double>& x,
334 Vector<std::complex<double>>& pml_inverse_jacobian_diagonal,
335 std::complex<double>& pml_jacobian_det)
336 {
337 /// The factors all default to 1.0 if the propagation
338 /// medium is the physical domain (no PML transformation)
339 for (unsigned k = 0; k < DIM; k++)
340 {
341 pml_inverse_jacobian_diagonal[k] = std::complex<double>(1.0, 0.0);
342 }
343 pml_jacobian_det = std::complex<double>(1.0, 0.0);
344
345 // Only calculate PML factors if PML is enabled
346 if (this->Pml_is_enabled)
347 {
348 /// Vector which points from the inner boundary to x
349 Vector<double> nu(DIM);
350 for (unsigned k = 0; k < DIM; k++)
351 {
352 nu[k] = x[k] - this->Pml_inner_boundary[k];
353 }
354
355 /// Vector which points from the inner boundary to the edge of the
356 /// boundary
357 Vector<double> pml_width(DIM);
358 for (unsigned k = 0; k < DIM; k++)
359 {
360 pml_width[k] =
361 this->Pml_outer_boundary[k] - this->Pml_inner_boundary[k];
362 }
363
364#ifdef PARANOID
365 // Check if the Pml_mapping_pt is set
366 if (this->Pml_mapping_pt == 0)
367 {
368 std::ostringstream error_message_stream;
369 error_message_stream << "Pml_mapping_pt needs to be set "
370 << std::endl;
371
372 throw OomphLibError(error_message_stream.str(),
373 OOMPH_CURRENT_FUNCTION,
374 OOMPH_EXCEPTION_LOCATION);
375 }
376#endif
377 // Declare gamma_i vectors of complex numbers for PML weights
378 Vector<std::complex<double>> pml_gamma(DIM);
379
380 /// Calculate the square of the non-dimensional wavenumber
381 double wavenumber_squared = 2.0 * (1.0 + this->nu()) * this->omega_sq();
382
383 for (unsigned k = 0; k < DIM; k++)
384 {
385 // If PML is enabled in the respective direction
386 if (this->Pml_direction_active[k])
387 {
388 std::complex<double> pml_gamma =
389 Pml_mapping_pt->gamma(nu[k], pml_width[k], wavenumber_squared);
390
391 // The diagonals of the INVERSE of the PML transformation jacobian
392 // are 1/gamma
393 pml_inverse_jacobian_diagonal[k] = 1.0 / pml_gamma;
394 // To get the determinant, multiply all the diagonals together
395 pml_jacobian_det *= pml_gamma;
396 }
397 }
398 }
399 }
400
401 /// Return a pointer to the PML Mapping object
403 {
404 return Pml_mapping_pt;
405 }
406
407 /// Return a pointer to the PML Mapping object (const version)
409 {
410 return Pml_mapping_pt;
411 }
412
413 /// Static so that the class doesn't need to instantiate a new default
414 /// everytime it uses it
416
417 protected:
418 /// Pointer to the elasticity tensor
420
421 /// Square of nondim frequency
422 double* Omega_sq_pt;
423
424 /// Pointer to body force function
426
427 /// Static default value for square of frequency
429
430 /// Pointer to class which holds the pml mapping function, also known as
431 /// gamma
433 };
434
435
436 /// ////////////////////////////////////////////////////////////////////
437 /// ////////////////////////////////////////////////////////////////////
438 /// ////////////////////////////////////////////////////////////////////
439
440
441 //=======================================================================
442 /// A class for elements that solve the equations of linear elasticity
443 /// in cartesian coordinates.
444 //=======================================================================
445 template<unsigned DIM>
448 {
449 public:
450 /// Constructor
452
453 /// Number of values required at node n.
454 unsigned required_nvalue(const unsigned& n) const
455 {
456 return 2 * DIM;
457 }
458
459 /// Return the residuals for the solid equations (the discretised
460 /// principle of virtual displacements)
462 {
465 }
466
467 /// The jacobian is calculated by finite differences by default,
468 /// We need only to take finite differences w.r.t. positional variables
469 /// For this element
471 DenseMatrix<double>& jacobian)
472 {
473 // Add the contribution to the residuals
474 this
475 ->fill_in_generic_contribution_to_residuals_time_harmonic_linear_elasticity(
476 residuals, jacobian, 1);
477 }
478
479 /// Return the Cauchy stress tensor, as calculated
480 /// from the elasticity tensor at specified local coordinate
481 void get_stress(const Vector<double>& s,
482 DenseMatrix<std::complex<double>>& sigma) const;
483
484 /// Output exact solution x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i]
485 void output_fct(std::ostream& outfile,
486 const unsigned& nplot,
488
489 /// Output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i]
490 void output(std::ostream& outfile)
491 {
492 unsigned n_plot = 5;
493 output(outfile, n_plot);
494 }
495
496 /// Output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i]
497 void output(std::ostream& outfile, const unsigned& n_plot);
498
499
500 /// C-style output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i]
501 void output(FILE* file_pt)
502 {
503 unsigned n_plot = 5;
504 output(file_pt, n_plot);
505 }
506
507 /// Output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i]
508 void output(FILE* file_pt, const unsigned& n_plot);
509
510 /// Output function for real part of full time-dependent solution
511 /// constructed by adding the scattered field
512 /// u = Re( (u_r +i u_i) exp(-i omega t)
513 /// at phase angle omega t = phi computed here, to the corresponding
514 /// incoming wave specified via the function pointer.
516 std::ostream& outfile,
517 FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt,
518 const double& phi,
519 const unsigned& nplot);
520
521
522 /// Output function for real part of full time-dependent solution
523 /// u = Re( (u_r +i u_i) exp(-i omega t))
524 /// at phase angle omega t = phi.
525 /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
526 /// direction
527 void output_real(std::ostream& outfile,
528 const double& phi,
529 const unsigned& n_plot);
530
531 /// Output function for imaginary part of full time-dependent
532 /// solution u = Im( (u_r +i u_i) exp(-i omega t) ) at phase angle omega t =
533 /// phi. x,y,u or x,y,z,u at n_plot plot points in each coordinate
534 /// direction
535 void output_imag(std::ostream& outfile,
536 const double& phi,
537 const unsigned& n_plot);
538
539
540 /// Compute norm of solution: square of the L2 norm
541 void compute_norm(double& norm);
542
543 /// Get error against and norm of exact solution
544 void compute_error(std::ostream& outfile,
546 double& error,
547 double& norm);
548
549
550 /// Dummy, time dependent error checker
551 void compute_error(std::ostream& outfile,
553 const double& time,
554 double& error,
555 double& norm)
556 {
557 std::ostringstream error_stream;
558 error_stream << "There is no time-dependent compute_error() " << std::endl
559 << "for pml time harmonic linear elasticity elements"
560 << std::endl;
561 throw OomphLibError(
562 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
563 }
564
565 private:
566 /// Private helper function to compute residuals and (if requested
567 /// via flag) also the Jacobian matrix.
569 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
570 };
571
572 /// /////////////////////////////////////////////////////////////////////
573 /// /////////////////////////////////////////////////////////////////////
574 /// /////////////////////////////////////////////////////////////////////
575
576
577 //===========================================================================
578 /// An Element that solves the equations of linear elasticity
579 /// in Cartesian coordinates, using QElements for the geometry
580 //============================================================================
581 template<unsigned DIM, unsigned NNODE_1D>
583 : public virtual QElement<DIM, NNODE_1D>,
585 {
586 public:
587 /// Constructor
589 : QElement<DIM, NNODE_1D>(),
591 {
592 }
593
594 /// Output function
595 void output(std::ostream& outfile)
596 {
598 }
599
600 /// Output function
601 void output(std::ostream& outfile, const unsigned& n_plot)
602 {
604 }
605
606
607 /// C-style output function
608 void output(FILE* file_pt)
609 {
611 }
612
613 /// C-style output function
614 void output(FILE* file_pt, const unsigned& n_plot)
615 {
617 }
618 };
619
620
621 //============================================================================
622 /// FaceGeometry of a linear 2D
623 /// QPMLTimeHarmonicLinearElasticityElement element
624 //============================================================================
625 template<>
627 : public virtual QElement<1, 2>
628 {
629 public:
630 /// Constructor must call the constructor of the underlying solid element
631 FaceGeometry() : QElement<1, 2>() {}
632 };
633
634
635 //============================================================================
636 /// FaceGeometry of a quadratic 2D
637 /// QPMLTimeHarmonicLinearElasticityElement element
638 //============================================================================
639 template<>
641 : public virtual QElement<1, 3>
642 {
643 public:
644 /// Constructor must call the constructor of the underlying element
645 FaceGeometry() : QElement<1, 3>() {}
646 };
647
648
649 //============================================================================
650 /// FaceGeometry of a cubic 2D
651 /// QPMLTimeHarmonicLinearElasticityElement element
652 //============================================================================
653 template<>
655 : public virtual QElement<1, 4>
656 {
657 public:
658 /// Constructor must call the constructor of the underlying element
659 FaceGeometry() : QElement<1, 4>() {}
660 };
661
662
663 //============================================================================
664 /// FaceGeometry of a linear 3D
665 /// QPMLTimeHarmonicLinearElasticityElement element
666 //============================================================================
667 template<>
669 : public virtual QElement<2, 2>
670 {
671 public:
672 /// Constructor must call the constructor of the underlying element
673 FaceGeometry() : QElement<2, 2>() {}
674 };
675
676 //============================================================================
677 /// FaceGeometry of a quadratic 3D
678 /// QPMLTimeHarmonicLinearElasticityElement element
679 //============================================================================
680 template<>
682 : public virtual QElement<2, 3>
683 {
684 public:
685 /// Constructor must call the constructor of the underlying element
686 FaceGeometry() : QElement<2, 3>() {}
687 };
688
689
690 //============================================================================
691 /// FaceGeometry of a cubic 3D
692 /// QPMLTimeHarmonicLinearElasticityElement element
693 //============================================================================
694 template<>
696 : public virtual QElement<2, 4>
697 {
698 public:
699 /// Constructor must call the constructor of the underlying element
700 FaceGeometry() : QElement<2, 4>() {}
701 };
702
703 /// /////////////////////////////////////////////////////////////////
704 /// /////////////////////////////////////////////////////////////////
705 /// /////////////////////////////////////////////////////////////////
706
707
708 //==========================================================
709 /// Time-harmonic linear elasticity upgraded to become projectable
710 //==========================================================
711 template<class TIME_HARMONIC_LINEAR_ELAST_ELEMENT>
713 : public virtual ProjectableElement<TIME_HARMONIC_LINEAR_ELAST_ELEMENT>
714 {
715 public:
716 /// Constructor [this was only required explicitly
717 /// from gcc 4.5.2 onwards...]
719
720
721 /// Specify the values associated with field fld.
722 /// The information is returned in a vector of pairs which comprise
723 /// the Data object and the value within it, that correspond to field fld.
724 /// In the underlying time-harmonic linear elasticity elemements the
725 /// real and complex parts of the displacements are stored
726 /// at the nodal values
728 {
729 // Create the vector
731
732 // Loop over all nodes and extract the fld-th nodal value
733 unsigned nnod = this->nnode();
734 for (unsigned j = 0; j < nnod; j++)
735 {
736 // Add the data value associated with the velocity components
737 data_values.push_back(std::make_pair(this->node_pt(j), fld));
738 }
739
740 // Return the vector
741 return data_values;
742 }
743
744 /// Number of fields to be projected: 2*dim, corresponding to
745 /// real and imag parts of the displacement components
747 {
748 return 2 * this->dim();
749 }
750
751 /// Number of history values to be stored for fld-th field.
752 /// (includes present value!)
753 unsigned nhistory_values_for_projection(const unsigned& fld)
754 {
755#ifdef PARANOID
756 if (fld > 3)
757 {
758 std::stringstream error_stream;
759 error_stream << "Elements only store four fields so fld can't be"
760 << " " << fld << std::endl;
761 throw OomphLibError(
762 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
763 }
764#endif
765 return this->node_pt(0)->ntstorage();
766 }
767
768 /// Number of positional history values: Read out from
769 /// positional timestepper
770 /// (Note: count includes current value!)
772 {
773 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
774 }
775
776 /// Return Jacobian of mapping and shape functions of field fld
777 /// at local coordinate s
778 double jacobian_and_shape_of_field(const unsigned& fld,
779 const Vector<double>& s,
780 Shape& psi)
781 {
782 unsigned n_dim = this->dim();
783 unsigned n_node = this->nnode();
784 DShape dpsidx(n_node, n_dim);
785
786 // Call the derivatives of the shape functions and return
787 // the Jacobian
788 return this->dshape_eulerian(s, psi, dpsidx);
789 }
790
791
792 /// Return interpolated field fld at local coordinate s, at time
793 /// level t (t=0: present; t>0: history values)
794 double get_field(const unsigned& t,
795 const unsigned& fld,
796 const Vector<double>& s)
797 {
798 unsigned n_node = this->nnode();
799
800#ifdef PARANOID
801 unsigned n_dim = this->node_pt(0)->ndim();
802#endif
803
804 // Local shape function
805 Shape psi(n_node);
806
807 // Find values of shape function
808 this->shape(s, psi);
809
810 // Initialise value of u
811 double interpolated_u = 0.0;
812
813 // Sum over the local nodes at that time
814 for (unsigned l = 0; l < n_node; l++)
815 {
816#ifdef PARANOID
817 unsigned nvalue = this->node_pt(l)->nvalue();
818 if (nvalue != 2 * n_dim)
819 {
820 std::stringstream error_stream;
821 error_stream
822 << "Current implementation only works for non-resized nodes\n"
823 << "but nvalue= " << nvalue << "!= 2 dim = " << 2 * n_dim
824 << std::endl;
825 throw OomphLibError(error_stream.str(),
826 OOMPH_CURRENT_FUNCTION,
827 OOMPH_EXCEPTION_LOCATION);
828 }
829#endif
830 interpolated_u += this->nodal_value(t, l, fld) * psi[l];
831 }
832 return interpolated_u;
833 }
834
835
836 /// Return number of values in field fld
837 unsigned nvalue_of_field(const unsigned& fld)
838 {
839 return this->nnode();
840 }
841
842
843 /// Return local equation number of value j in field fld.
844 int local_equation(const unsigned& fld, const unsigned& j)
845 {
846#ifdef PARANOID
847 unsigned n_dim = this->node_pt(0)->ndim();
848 unsigned nvalue = this->node_pt(j)->nvalue();
849 if (nvalue != 2 * n_dim)
850 {
851 std::stringstream error_stream;
852 error_stream
853 << "Current implementation only works for non-resized nodes\n"
854 << "but nvalue= " << nvalue << "!= 2 dim = " << 2 * n_dim
855 << std::endl;
856 throw OomphLibError(
857 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
858 }
859#endif
860 return this->nodal_local_eqn(j, fld);
861 }
862 };
863
864
865 //=======================================================================
866 /// Face geometry for element is the same as that for the underlying
867 /// wrapped element
868 //=======================================================================
869 template<class ELEMENT>
871 : public virtual FaceGeometry<ELEMENT>
872 {
873 public:
874 FaceGeometry() : FaceGeometry<ELEMENT>() {}
875 };
876
877
878 //=======================================================================
879 /// Face geometry of the Face Geometry for element is the same as
880 /// that for the underlying wrapped element
881 //=======================================================================
882 template<class ELEMENT>
885 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
886 {
887 public:
889 };
890
891
892 /// ///////////////////////////////////////////////////////////////////
893 /// ///////////////////////////////////////////////////////////////////
894 /// ///////////////////////////////////////////////////////////////////
895
896
897 //=======================================================================
898 /// Policy class defining the elements to be used in the actual
899 /// PML layers. Same!
900 //=======================================================================
901 template<unsigned NNODE_1D>
903 : public virtual QPMLTimeHarmonicLinearElasticityElement<2, NNODE_1D>
904 {
905 public:
906 /// Constructor: Call the constructor for the
907 /// appropriate QElement
909 {
910 }
911 };
912
913
914} // namespace oomph
915
916
917#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A mapping function proposed by Bermudez et al, similar to the one above but is continuous across the ...
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:879
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
Definition: matrices.h:386
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
FaceGeometry()
Constructor must call the constructor of the underlying element.
FaceGeometry()
Constructor must call the constructor of the underlying element.
FaceGeometry()
Constructor must call the constructor of the underlying element.
FaceGeometry()
Constructor must call the constructor of the underlying element.
FaceGeometry()
Constructor must call the constructor of the underlying element.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
An OomphLibError object which should be thrown when an run-time error is encountered....
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: pml_meshes.h:60
bool Pml_is_enabled
Boolean indicating if element is used in pml mode.
Definition: pml_meshes.h:119
std::vector< bool > Pml_direction_active
Coordinate direction along which pml boundary is constant; alternatively: coordinate direction in whi...
Definition: pml_meshes.h:124
Vector< double > Pml_outer_boundary
Coordinate of outer pml boundary (Storage is provided for any coordinate direction; only the entries ...
Definition: pml_meshes.h:134
Vector< double > Pml_inner_boundary
Coordinate of inner pml boundary (Storage is provided for any coordinate direction; only the entries ...
Definition: pml_meshes.h:129
General definition of policy class defining the elements to be used in the actual PML layers....
Definition: pml_meshes.h:48
Class to hold the mapping function (gamma) for the Pml which defines how the coordinates are transfor...
virtual std::complex< double > gamma(const double &nu_i, const double &pml_width_i, const double &wavenumber_squared, const double &alpha_shift=0.0)=0
Pure virtual to return Pml mapping gamma, where gamma is the as function of where where h is the v...
An isotropic elasticity tensor defined in terms of Young's modulus and Poisson's ratio....
A base class for elements that solve the equations of time-harmonic linear elasticity in Cartesian co...
void interpolated_u_time_harmonic_linear_elasticity(const Vector< double > &s, Vector< std::complex< double > > &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
BodyForceFctPt body_force_fct_pt() const
Access function: Pointer to body force function (const version)
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: for now lump ...
PMLTimeHarmonicIsotropicElasticityTensor *& elasticity_tensor_pt()
Return the pointer to the elasticity_tensor.
void body_force(const Vector< double > &x, Vector< std::complex< double > > &b) const
Evaluate body force at Eulerian coordinate x (returns zero vector if no body force function pointer h...
PMLMapping *const & pml_mapping_pt() const
Return a pointer to the PML Mapping object (const version)
void(* BodyForceFctPt)(const Vector< double > &x, Vector< std::complex< double > > &b)
Function pointer to function that specifies the body force as a function of the Cartesian coordinates...
std::complex< double > E(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l) const
Access function to the entries in the elasticity tensor.
double nu() const
Access function to nu in the elasticity tensor.
static ContinuousBermudezPMLMapping Default_pml_mapping
Static so that the class doesn't need to instantiate a new default everytime it uses it.
const double & omega_sq() const
Access function for square of non-dim frequency.
void values_to_be_pinned_on_outer_pml_boundary(Vector< unsigned > &values_to_pin)
Pure virtual function in which we specify the values to be pinned (and set to zero) on the outer edge...
static double Default_omega_sq_value
Static default value for square of frequency.
PMLMapping * Pml_mapping_pt
Pointer to class which holds the pml mapping function, also known as gamma.
PMLTimeHarmonicIsotropicElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
virtual std::complex< unsigned > u_index_time_harmonic_linear_elasticity(const unsigned i) const
Return the index at which the i-th real or imag unknown displacement component is stored....
void compute_pml_coefficients(const unsigned &ipt, const Vector< double > &x, Vector< std::complex< double > > &pml_inverse_jacobian_diagonal, std::complex< double > &pml_jacobian_det)
Compute pml coefficients at position x and integration point ipt. pml_inverse_jacobian_diagonal conta...
PMLTimeHarmonicLinearElasticityEquationsBase()
Constructor: Set null pointers for constitutive law and for isotropic growth function....
std::complex< double > interpolated_u_time_harmonic_linear_elasticity(const Vector< double > &s, const unsigned &i) const
Return FE interpolated displacement u[i] at local coordinate s.
virtual void get_stress(const Vector< double > &s, DenseMatrix< std::complex< double > > &sigma) const =0
Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordina...
double *& omega_sq_pt()
Access function for square of non-dim frequency.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
PMLMapping *& pml_mapping_pt()
Return a pointer to the PML Mapping object.
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double > > &strain) const
Return the strain tensor.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
The jacobian is calculated by finite differences by default, We need only to take finite differences ...
void output_imag(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for imaginary part of full time-dependent solution u = Im( (u_r +i u_i) exp(-i omega ...
void get_stress(const Vector< double > &s, DenseMatrix< std::complex< double > > &sigma) const
Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordina...
void compute_norm(double &norm)
Compute norm of solution: square of the L2 norm.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals for the solid equations (the discretised principle of virtual displacements)
unsigned required_nvalue(const unsigned &n) const
Number of values required at node n.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
void output(std::ostream &outfile)
Output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i].
virtual void fill_in_generic_contribution_to_residuals_time_harmonic_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Private helper function to compute residuals and (if requested via flag) also the Jacobian matrix.
void output_total_real(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt, const double &phi, const unsigned &nplot)
Output function for real part of full time-dependent solution constructed by adding the scattered fie...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i].
void output_real(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for real part of full time-dependent solution u = Re( (u_r +i u_i) exp(-i omega t)) a...
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void output(FILE *file_pt)
C-style output: x,y,[z],u_r,v_r,[w_r],u_i,v_i,[w_i].
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (includes present value!)
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
unsigned nfields_for_projection()
Number of fields to be projected: 2*dim, corresponding to real and imag parts of the displacement com...
ProjectablePMLTimeHarmonicLinearElasticityElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values: Read out from positional timestepper (Note: count includes curre...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
//////////////////////////////////////////////////////////////////// ////////////////////////////////...