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