pml_fourier_decomposed_helmholtz_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: 1299 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-10-06 07:40:18 +0100 (Fri, 06 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 Fourier-decomposed Helmholtz elements
31 #ifndef OOMPH_PML_FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENTS_HEADER
32 #define OOMPH_PML_FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 #include "math.h"
41 #include <complex>
42 
43 
44 //OOMPH-LIB headers
45 #include "../generic/projection.h"
46 #include "../generic/nodes.h"
47 #include "../generic/Qelements.h"
48 #include "../generic/oomph_utilities.h"
49 #include "../generic/pml_meshes.h"
50 #include "../generic/projection.h"
51 #include "../generic/oomph_definitions.h"
52 
53 
54 namespace oomph
55 {
56 
57 
58 //========================================================================
59 /// Helper namespace for functions required for Helmholtz computations
60 //========================================================================
61  namespace Legendre_functions_helper
62  {
63 
64  /// Factorial
65  extern double factorial(const unsigned& l);
66 
67  /// Legendre polynomials depending on one parameter
68  extern double plgndr1(const unsigned& n, const double& x);
69 
70  /// Legendre polynomials depending on two parameters
71  extern double plgndr2(const unsigned& l, const unsigned& m, const double& x);
72 
73  } // end namespace
74 
75 
76 ///////////////////////////////////////////////////////////////////////
77 ///////////////////////////////////////////////////////////////////////
78 ///////////////////////////////////////////////////////////////////////
79 
80 //=======================================================================
81 /// Class to hold the mapping function for the PML
82 ///
83 //=======================================================================
85 {
86 
87 public:
88 
89  /// Default constructor (empty)
91 
92  /// \short Pure virtual to return PML mapping gamma, where gamma is the
93  /// \f$d\tilde x / d x\f$ as function of \f$\nu\f$ where \f$\nu = x - h\f$ where h is
94  /// the vector from the origin to the start of the PML
95  virtual std::complex<double> gamma(const double& nu_i,
96  const double& pml_width_i,
97  const double& k_squared) = 0;
98 
99  /// \short Pure virtual to return PML transformed coordinate, also known as
100  /// \f$d\tilde x \f$ as function of \f$\nu\f$ where \f$\nu = x - h\f$ where h
101  /// is the vector from the origin to the start of the PML
102  virtual std::complex<double> transformed_coordinate(const double& nu_i,
103  const double& pml_width_i,
104  const double& pml_inner_boundary,
105  const double& k_squared) = 0;
106 
107 };
108 
109 
110 //=======================================================================
111 /// The mapping function propsed by Bermudez et al, appears to be the best
112 /// and so this will be the default mapping (see definition of
113 /// PMLHelmholtzEquations)
114 //=======================================================================
117 {
118 
119  public:
120 
121  /// Default constructor (empty)
123 
124  /// \short Overwrite the pure PML mapping coefficient function to return the
125  /// mapping function proposed by Bermudez et al
126  std::complex<double> gamma(const double& nu_i,
127  const double& pml_width_i,
128  const double& k_squared)
129  {
130  /// \short return \f$\gamma=1 + (1/k)(i/|outer_boundary - x|)\f$ or more
131  /// abstractly \f$\gamma = 1 + \frac i {k\delta_{pml}}(1/|1-\bar\nu|)\f$
132  return 1.0 + std::complex<double> (1.0 / sqrt(k_squared), 0)
133  * std::complex<double> (0.0, 1.0/(std::fabs(pml_width_i - nu_i)));
134  }
135 
136  /// \short Overwrite the pure PML mapping coefficient function to return the
137  /// transformed coordinate proposed by Bermudez et al
138  std::complex<double> transformed_coordinate(const double& nu_i,
139  const double& pml_width_i,
140  const double& pml_inner_boundary,
141  const double& k_squared)
142  {
143  /// \short return \f$\tilde x = h + \nu + \log(1-|\nu / \delta|)\f$
144  double log_arg = 1.0 - std::fabs(nu_i / pml_width_i);
145  return std::complex<double> (pml_inner_boundary + nu_i,
146  -log(log_arg)/sqrt(k_squared) );
147  }
148 
149 };
150 
151 
152 ////////////////////////////////////////////////////////////////////////
153 ////////////////////////////////////////////////////////////////////////
154 ////////////////////////////////////////////////////////////////////////
155 
156 
157 
158 //=============================================================
159 /// A class for all isoparametric elements that solve the
160 /// Helmholtz equations with pml capabilities.
161 /// in Fourier decomposed form (cylindrical polars):
162 /// \f[
163 /// U(r,\varphi,z) = \Re( u^{(n)}(r,z) \exp(-i n \varphi))
164 /// \f]
165 /// We are solving for \f$ u^{(n)}(r,z)\f$ for given parameters
166 /// \f$ k^2 \f$ and \f$ n \f$ .
167 /// This contains the generic maths. Shape functions, geometric
168 /// mapping etc. must get implemented in derived class.
169 //=============================================================
171  public virtual PMLElementBase<2>,
172  public virtual FiniteElement
173  {
174 
175  public:
176 
177  /// \short Function pointer to source function fct(x,f(x)) --
178  /// x is a Vector!
180  const Vector<double>& x,
181  std::complex<double>& f);
182 
183  /// Constructor
186  {
187  // Provide pointer to static method (Save memory)
191 
193  }
194 
195 
196  /// Broken copy constructor
199  {
200  BrokenCopy::broken_copy("PMLFourierDecomposedHelmholtzEquations");
201  }
202 
203  /// Broken assignment operator
204 //Commented out broken assignment operator because this can lead to a conflict warning
205 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
206 //realise that two separate implementations of the broken function are the same and so,
207 //quite rightly, it shouts.
208  /*void operator=(const PMLFourierDecomposedHelmholtzEquations&)
209  {
210  BrokenCopy::broken_assign
211  ("PMLFourierDecomposedHelmholtzEquations");
212  }*/
213 
214  /// \short Return the index at which the unknown value
215  /// is stored: Real/imag part of index contains (real) index of
216  /// real/imag part.
217  virtual inline std::complex<unsigned>
219  const
220  {return std::complex<unsigned>(0,1);}
221 
222 
223  /// Get pointer to frequency
224  double*& k_squared_pt()
225  {
226  return K_squared_pt;
227  }
228 
229 
230  /// Get k squared
231  double k_squared()
232  {
233  #ifdef PARANOID
234  if (K_squared_pt==0)
235  {
236  throw OomphLibError(
237  "Please set pointer to k_squared using access fct to pointer!",
238  OOMPH_CURRENT_FUNCTION,
239  OOMPH_EXCEPTION_LOCATION);
240  }
241  #endif
242  return *K_squared_pt;
243  }
244 
245  /// Get pointer to complex shift
246  double*& alpha_pt()
247  {
248  return Alpha_pt;
249  }
250 
251 
252  /// Get complex shift
253  double alpha()
254  {
255  return *Alpha_pt;
256  }
257 
258  /// Get pointer to Fourier wavenumber
260  {
261  return N_pml_fourier_pt;
262  }
263 
264  /// Get the Fourier wavenumber
266  {
267  if (N_pml_fourier_pt==0)
268  {
269  return 0;
270  }
271  else
272  {
273  return *N_pml_fourier_pt;
274  }
275  }
276 
277 
278  /// Output with default number of plot points
279  void output(std::ostream &outfile)
280  {
281  const unsigned n_plot=5;
282  output(outfile,n_plot);
283  }
284 
285  /// \short Output FE representation of soln: x,y,u_re,u_im or
286  /// x,y,z,u_re,u_im at n_plot^2 plot points
287  void output(std::ostream &outfile, const unsigned &n_plot);
288 
289  /// \short Output function for real part of full time-dependent solution
290  /// u = Re( (u_r +i u_i) exp(-i omega t)
291  /// at phase angle omega t = phi.
292  /// r,z,u at n_plot plot points in each coordinate
293  /// direction
294  void output_real(std::ostream &outfile, const double& phi,
295  const unsigned &n_plot);
296 
297  /// C_style output with default number of plot points
298  void output(FILE* file_pt)
299  {
300  const unsigned n_plot=5;
301  output(file_pt,n_plot);
302  }
303 
304  /// \short C-style output FE representation of soln: r,z,u_re,u_im or
305  /// at n_plot^2 plot points
306  void output(FILE* file_pt, const unsigned &n_plot);
307 
308  /// Output exact soln: r,z,u_re_exact,u_im_exact
309  /// at n_plot^2 plot points
310  void output_fct(std::ostream &outfile, const unsigned &n_plot,
312 
313  /// \short Output exact soln: (dummy time-dependent version to
314  /// keep intel compiler happy)
315  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
316  const double& time,
318  exact_soln_pt)
319  {
320  throw OomphLibError(
321  "There is no time-dependent output_fct() for PMLFourierDecomposedHelmholtz elements ",
322  OOMPH_CURRENT_FUNCTION,
323  OOMPH_EXCEPTION_LOCATION);
324  }
325 
326 
327  /// \short Output function for real part of full time-dependent fct
328  /// u = Re( (u_r +i u_i) exp(-i omega t)
329  /// at phase angle omega t = phi.
330  /// r,z,u at n_plot plot points in each coordinate
331  /// direction
332  void output_real_fct(std::ostream &outfile,
333  const double& phi,
334  const unsigned &n_plot,
336 
337 
338  /// Get error against and norm of exact solution
339  void compute_error(std::ostream &outfile,
341  double& error, double& norm);
342 
343 
344  /// Dummy, time dependent error checker
345  void compute_error(std::ostream &outfile,
347  const double& time, double& error, double& norm)
348  {
349  throw OomphLibError(
350  "There is no time-dependent compute_error() for PMLFourierDecomposedHelmholtz elements",
351  OOMPH_CURRENT_FUNCTION,
352  OOMPH_EXCEPTION_LOCATION);
353  }
354 
355  /// Compute norm of fe solution
356  void compute_norm(double& norm);
357 
358 
359  /// Access function: Pointer to source function
361  {return Source_fct_pt;}
362 
363 
364  /// Access function: Pointer to source function. Const version
366  {return Source_fct_pt;}
367 
368 
369  /// Get source term at (Eulerian) position x. This function is
370  /// virtual to allow overloading in multi-physics problems where
371  /// the strength of the source function might be determined by
372  /// another system of equations.
374  const unsigned& ipt,
375  const Vector<double>& x,
376  std::complex<double>& source) const
377  {
378  //If no source function has been set, return zero
379  if(Source_fct_pt==0)
380  {
381  source = std::complex<double>(0.0,0.0);
382  }
383  else
384  {
385  // Get source strength
386  (*Source_fct_pt)(x,source);
387  }
388  }
389 
390  /// \short Pure virtual function in which we specify the
391  /// values to be pinned (and set to zero) on the outer edge of
392  /// the pml layer. All of them! Vector is resized internally.
394  values_to_pin)
395  {
396  values_to_pin.resize(2);
397  for (unsigned j=0;j<2;j++)
398  {
399  values_to_pin[j]=j;
400  }
401  }
402 
403 
404  /// Get flux: flux[i] = du/dx_i for real and imag part
405  void get_flux(const Vector<double>& s,
406  Vector<std::complex <double> >& flux) const
407  {
408  //Find out how many nodes there are in the element
409  const unsigned n_node = nnode();
410 
411  //Set up memory for the shape and test functions
412  Shape psi(n_node);
413  DShape dpsidx(n_node,2);
414 
415  //Call the derivatives of the shape and test functions
416  dshape_eulerian(s,psi,dpsidx);
417 
418  //Initialise to zero
419  const std::complex<double> zero(0.0,0.0);
420  for(unsigned j=0;j<2;j++)
421  {
422  flux[j] = zero;
423  }
424 
425  // Loop over nodes
426  for(unsigned l=0;l<n_node;l++)
427  {
428  //Cache the complex value of the unknown
429  const std::complex<double> u_value(
430  this->nodal_value
432  this->nodal_value
434 
435  //Loop over derivative directions
436  for(unsigned j=0;j<2;j++)
437  {
438  flux[j] += u_value*dpsidx(l,j);
439  }
440  }
441  }
442 
443 
444  /// Add the element's contribution to its residual vector (wrapper)
446  {
447  //Call the generic residuals function with flag set to 0
448  //using a dummy matrix argument
451  }
452 
453 
454 
455  /// \short Add the element's contribution to its residual vector and
456  /// element Jacobian matrix (wrapper)
458  DenseMatrix<double> &jacobian)
459  {
460  //Call the generic routine with the flag set to 1
462  residuals,jacobian,1);
463  }
464 
465 
466  /// \short Return FE representation of function value u(s)
467  /// at local coordinate s
468  inline std::complex<double>
470  const Vector<double> &s) const
471  {
472  //Find number of nodes
473  const unsigned n_node = nnode();
474 
475  //Local shape function
476  Shape psi(n_node);
477 
478  //Find values of shape function
479  shape(s,psi);
480 
481  //Initialise value of u
482  std::complex<double> interpolated_u(0.0,0.0);
483 
484  //Get the index at which the helmholtz unknown is stored
485  const unsigned u_nodal_index_real =
487  const unsigned u_nodal_index_imag =
489 
490  //Loop over the local nodes and sum
491  for(unsigned l=0;l<n_node;l++)
492  {
493  //Make a temporary complex number from the stored data
494  const std::complex<double> u_value(
495  this->nodal_value(l,u_nodal_index_real),
496  this->nodal_value(l,u_nodal_index_imag));
497  //Add to the interpolated value
498  interpolated_u += u_value*psi[l];
499  }
500  return interpolated_u;
501  }
502 
503 
504  /// \short Self-test: Return 0 for OK
505  unsigned self_test();
506 
507 
508 protected:
509 
510 
511 
512  /// \short Compute pml coefficients at position x and integration point ipt.
513  /// pml_laplace_factor is used in the residual contribution from the laplace
514  /// operator, similarly pml_k_squared_factor is used in the contribution from
515  /// the k^2 of the Helmholtz operator.
517  const unsigned& ipt,
518  const Vector<double>& x,
519  Vector< std::complex<double> >& pml_laplace_factor,
520  std::complex<double>& pml_k_squared_factor)
521  {
522 
523  /// Vector which points from the inner boundary to x
524  Vector<double> nu(2);
525  for(unsigned k=0; k<2; k++)
526  {
527  nu[k] = x[k] - this->Pml_inner_boundary[k];
528  }
529 
530  /// Vector which points from the inner boundary to the edge of the boundary
531  Vector<double> pml_width(2);
532  for(unsigned k=0; k<2; k++)
533  {
534  pml_width[k] = this->Pml_outer_boundary[k] - this->Pml_inner_boundary[k];
535  }
536 
537  // Declare gamma_i vectors of complex numbers for PML weights
538  Vector<std::complex<double> > pml_gamma(2);
539 
540  if (this->Pml_is_enabled)
541  {
542  // Cache k_squared to pass into mapping function
543  double k_squared_local = k_squared();
544 
545  for(unsigned k=0; k<2; k++)
546  {
547  // If PML is enabled in the respective direction
548  if (this->Pml_direction_active[k])
549  {
551  gamma(nu[k], pml_width[k], k_squared_local);
552  }
553  else
554  {
555  pml_gamma[k] = 1.0;
556  }
557  }
558 
559  /// \short for 2D, in order:
560  /// g_y/g_x, g_x/g_y for Laplace bit and g_x*g_y for Helmholtz bit
561  /// for 3D, in order: g_y*g_x/g_x, g*x*g_z/g_y, g_x*g_y/g_z for Laplace bit
562  /// and g_x*g_y*g_z for Helmholtz factor
563  pml_laplace_factor[0] = pml_gamma[1] / pml_gamma[0];
564  pml_laplace_factor[1] = pml_gamma[0] / pml_gamma[1];
565  pml_k_squared_factor = pml_gamma[0] * pml_gamma[1];
566 
567 
568  }
569  else
570  {
571  /// \short The weights all default to 1.0 as if the propagation
572  /// medium is the physical domain
573  for(unsigned k=0; k<2; k++)
574  {
575  pml_laplace_factor[k] = std::complex<double> (1.0,0.0);
576  }
577 
578  pml_k_squared_factor = std::complex<double> (1.0,0.0);
579  }
580 
581  }
582 
583 
584 /// \short Compute complex variable r at position x[0] and
585 /// integration point ipt
586  void compute_complex_r(const unsigned& ipt,
587  const Vector<double>& x,
588  std::complex<double>& complex_r)
589  {
590  // Cache current position r
591  double r = x[0];
592 
593  /// The complex r variable is only imaginary on two
594  /// conditions: First, the decaying nature of the
595  /// pml layers is active. Secondly, the
596  /// integration point is contained in the right pml
597  /// layer or the two corner pml layers.
598 
599  // If the complex r variable is imaginary
600  if (this->Pml_is_enabled && (this->Pml_direction_active[0]))
601  {
602  double nu = x[0] - Pml_inner_boundary[0];
603  double pml_width = Pml_outer_boundary[0] - Pml_inner_boundary[0];
604  double k_squared_local = k_squared();
605 
606  // Determine the complex r variable
608  transformed_coordinate(nu,pml_width,Pml_inner_boundary[0],k_squared_local);
609  }
610  else
611  {
612  // The complex r variable is infact purely real, and
613  // is equal to x[0]
614  complex_r = std::complex<double> (r,0.0);
615  }
616 
617  } // end of compute_complex_r
618 
619  /// Return a pointer to the PML Mapping object
621  {
623  }
624 
625  /// Return a pointer to the PML Mapping object (const version)
628  {
630  }
631 
632  /// Static so that the class doesn't need to instantiate a new default
633  /// everytime it uses it
636 
637 
638  /// \short Shape/test functions and derivs w.r.t. to global coords at
639  /// local coord. s; return Jacobian of mapping
640  virtual double
642  const Vector<double> &s,
643  Shape &psi,
644  DShape &dpsidx, Shape &test,
645  DShape &dtestdx) const=0;
646 
647 
648  /// \short Shape/test functions and derivs w.r.t. to global coords at
649  /// integration point ipt; return Jacobian of mapping
650  virtual double
652  const unsigned &ipt,
653  Shape &psi,
654  DShape &dpsidx,
655  Shape &test,
656  DShape &dtestdx) const=0;
657 
658  /// \short Compute element residual Vector only (if flag=and/or element
659  /// Jacobian matrix
660  virtual void
662  Vector<double> &residuals, DenseMatrix<double> &jacobian,
663  const unsigned& flag);
664 
665  /// Pointer to source function:
667 
668  /// Pointer to k^2 (wavenumber squared)
669  double* K_squared_pt;
670 
671  /// \short Pointer to class which holds the pml mapping function (also known
672  /// as gamma) and the associated transformed coordinate
674 
675  /// Pointer to wavenumber complex shift
676  double *Alpha_pt;
677 
678  /// Static default value for the physical constants (initialised to zero)
680 
681  /// Pointer to Fourier wave number
683 
684  };
685 
686 
687 
688 
689 ///////////////////////////////////////////////////////////////////////////
690 ///////////////////////////////////////////////////////////////////////////
691 ///////////////////////////////////////////////////////////////////////////
692 
693 
694 
695 //======================================================================
696 /// QPMLFourierDecomposedHelmholtzElement elements are
697 /// linear/quadrilateral/brick-shaped PMLFourierDecomposedHelmholtz
698 /// elements with isoparametric interpolation for the function.
699 //======================================================================
700  template <unsigned NNODE_1D>
702  public virtual QElement<2,NNODE_1D>,
704  {
705 
706  private:
707 
708  /// \short Static int that holds the number of variables at
709  /// nodes: always the same
710  static const unsigned Initial_Nvalue;
711 
712  public:
713 
714 
715  ///\short Constructor: Call constructors for QElement and
716  /// PMLFourierDecomposedHelmholtz equations
719  {}
720 
721  /// Broken copy constructor
724  {
725  BrokenCopy::broken_copy("QPMLFourierDecomposedHelmholtzElement");
726  }
727 
728  /// Broken assignment operator
729  /*void operator=(const
730  QPMLFourierDecomposedHelmholtzElement<NNODE_1D>&)
731  {
732  BrokenCopy::broken_assign
733  ("QPMLFourierDecomposedHelmholtzElement");
734  }*/
735 
736 
737  /// \short Required # of `values' (pinned or dofs)
738  /// at node n
739  inline unsigned required_nvalue(const unsigned &n) const
740  {return Initial_Nvalue;}
741 
742  /// \short Output function: r,z,u
743  void output(std::ostream &outfile)
745 
746  /// \short Output function:
747  /// r,z,u at n_plot^2 plot points
748  void output(std::ostream &outfile, const unsigned &n_plot)
750 
751  /// \short Output function for real part of full time-dependent solution
752  /// u = Re( (u_r +i u_i) exp(-i omega t)
753  /// at phase angle omega t = phi.
754  /// r,z,u at n_plot plot points in each coordinate
755  /// direction
756  void output_real(std::ostream &outfile, const double& phi,
757  const unsigned &n_plot)
759  (outfile,phi,n_plot);}
760 
761  /// \short C-style output function: r,z,u
762  void output(FILE* file_pt)
764 
765  /// \short C-style output function:
766  /// r,z,u at n_plot^2 plot points
767  void output(FILE* file_pt, const unsigned &n_plot)
769 
770  /// \short Output function for an exact solution:
771  /// r,z,u_exact at n_plot^2 plot points
772  void output_fct(std::ostream &outfile, const unsigned &n_plot,
774  {
776  exact_soln_pt);
777  }
778 
779  /// \short Output function for real part of full time-dependent fct
780  /// u = Re( (u_r +i u_i) exp(-i omega t)
781  /// at phase angle omega t = phi.
782  /// r,z,u at n_plot plot points in each coordinate
783  /// direction
784  void output_real_fct(std::ostream &outfile,
785  const double& phi,
786  const unsigned &n_plot,
788  {
790  (outfile,phi,
791  n_plot,exact_soln_pt);
792  }
793 
794 
795  /// \short Output function for a time-dependent exact solution.
796  /// r,z,u_exact at n_plot^2 plot points
797  /// (Calls the steady version)
798  void output_fct(std::ostream &outfile, const unsigned &n_plot,
799  const double& time,
801  {
803  (outfile,n_plot,time,exact_soln_pt);
804  }
805 
806  protected:
807 
808  /// Shape, test functions & derivs. w.r.t. to global coords.
809  /// Return Jacobian.
810  inline double
812  const Vector<double> &s, Shape &psi, DShape &dpsidx,
813  Shape &test, DShape &dtestdx) const;
814 
815 
816  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
817  /// integration point ipt. Return Jacobian.
818  inline double
820  const unsigned& ipt,
821  Shape &psi,
822  DShape &dpsidx,
823  Shape &test,
824  DShape &dtestdx) const;
825 
826  };
827 
828 
829 
830 
831 //Inline functions:
832 
833 
834 //======================================================================
835 /// Define the shape functions and test functions and derivatives
836 /// w.r.t. global coordinates and return Jacobian of mapping.
837 ///
838 /// Galerkin: Test functions = shape functions
839 //======================================================================
840  template<unsigned NNODE_1D>
843  const Vector<double> &s,
844  Shape &psi,
845  DShape &dpsidx,
846  Shape &test,
847  DShape &dtestdx) const
848  {
849  //Call the geometrical shape functions and derivatives
850  const double J = this->dshape_eulerian(s,psi,dpsidx);
851 
852  //Set the test functions equal to the shape functions
853  test = psi;
854  dtestdx= dpsidx;
855 
856  //Return the jacobian
857  return J;
858  }
859 
860 
861 
862 
863 //======================================================================
864 /// Define the shape functions and test functions and derivatives
865 /// w.r.t. global coordinates and return Jacobian of mapping.
866 ///
867 /// Galerkin: Test functions = shape functions
868 //======================================================================
869  template<unsigned NNODE_1D>
872  const unsigned &ipt,
873  Shape &psi,
874  DShape &dpsidx,
875  Shape &test,
876  DShape &dtestdx) const
877  {
878  //Call the geometrical shape functions and derivatives
879  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
880 
881  //Set the pointers of the test functions
882  test = psi;
883  dtestdx = dpsidx;
884 
885  //Return the jacobian
886  return J;
887  }
888 
889 ////////////////////////////////////////////////////////////////////////
890 ////////////////////////////////////////////////////////////////////////
891 ////////////////////////////////////////////////////////////////////////
892 
893 
894 
895 //=======================================================================
896 /// Face geometry for the QPMLFourierDecomposedHelmholtzElement
897 /// elements:
898 /// The spatial dimension of the face elements is one lower than that of the
899 /// bulk element but they have the same number of points
900 /// along their 1D edges.
901 //=======================================================================
902  template<unsigned NNODE_1D>
904  public virtual QElement<1,NNODE_1D>
905  {
906 
907  public:
908 
909  /// \short Constructor: Call the constructor for the
910  /// appropriate lower-dimensional QElement
911  FaceGeometry() : QElement<1,NNODE_1D>() {}
912  };
913 
914 
915 ////////////////////////////////////////////////////////////////////////
916 ////////////////////////////////////////////////////////////////////////
917 ////////////////////////////////////////////////////////////////////////
918 
919 
920 
921 //==========================================================
922 /// Fourier decomposed Helmholtz upgraded to become projectable
923 //==========================================================
924  template<class FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENT>
926  public virtual ProjectableElement<FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENT>
927  {
928 
929  public:
930 
931 
932  /// \short Constructor [this was only required explicitly
933  /// from gcc 4.5.2 onwards...]
935 
936  /// \short Specify the values associated with field fld.
937  /// The information is returned in a vector of pairs which comprise
938  /// the Data object and the value within it, that correspond to field fld.
940  {
941 
942 #ifdef PARANOID
943  if (fld>1)
944  {
945  std::stringstream error_stream;
946  error_stream
947  << "Fourier decomposed Helmholtz elements only store 2 fields so fld = "
948  << fld << " is illegal \n";
949  throw OomphLibError(
950  error_stream.str(),
951  OOMPH_CURRENT_FUNCTION,
952  OOMPH_EXCEPTION_LOCATION);
953  }
954 #endif
955 
956  // Create the vector
957  unsigned nnod=this->nnode();
958  Vector<std::pair<Data*,unsigned> > data_values(nnod);
959 
960  // Loop over all nodes
961  for (unsigned j=0;j<nnod;j++)
962  {
963  // Add the data value associated field: The node itself
964  data_values[j]=std::make_pair(this->node_pt(j),fld);
965  }
966 
967  // Return the vector
968  return data_values;
969  }
970 
971  /// \short Number of fields to be projected: 2 (real and imag part)
973  {
974  return 2;
975  }
976 
977  /// \short Number of history values to be stored for fld-th field.
978  /// (Note: count includes current value!)
979  unsigned nhistory_values_for_projection(const unsigned &fld)
980  {
981 #ifdef PARANOID
982  if (fld>1)
983  {
984  std::stringstream error_stream;
985  error_stream
986  << "Helmholtz elements only store two fields so fld = "
987  << fld << " is illegal\n";
988  throw OomphLibError(
989  error_stream.str(),
990  OOMPH_CURRENT_FUNCTION,
991  OOMPH_EXCEPTION_LOCATION);
992  }
993 #endif
994  return this->node_pt(0)->ntstorage();
995  }
996 
997  ///\short Number of positional history values
998  /// (Note: count includes current value!)
1000  {
1001  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
1002  }
1003 
1004  /// \short Return Jacobian of mapping and shape functions of field fld
1005  /// at local coordinate s
1006  double jacobian_and_shape_of_field(const unsigned &fld,
1007  const Vector<double> &s,
1008  Shape &psi)
1009  {
1010 #ifdef PARANOID
1011  if (fld>1)
1012  {
1013  std::stringstream error_stream;
1014  error_stream
1015  << "Helmholtz elements only store two fields so fld = "
1016  << fld << " is illegal.\n";
1017  throw OomphLibError(
1018  error_stream.str(),
1019  OOMPH_CURRENT_FUNCTION,
1020  OOMPH_EXCEPTION_LOCATION);
1021  }
1022 #endif
1023  unsigned n_dim=this->dim();
1024  unsigned n_node=this->nnode();
1025  Shape test(n_node);
1026  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
1027  double J=this->
1028  dshape_and_dtest_eulerian_pml_fourier_decomposed_helmholtz(
1029  s,psi,dpsidx,
1030  test,dtestdx);
1031  return J;
1032  }
1033 
1034 
1035 
1036  /// \short Return interpolated field fld at local coordinate s, at time level
1037  /// t (t=0: present; t>0: history values)
1038  double get_field(const unsigned &t,
1039  const unsigned &fld,
1040  const Vector<double>& s)
1041  {
1042 #ifdef PARANOID
1043  if (fld>1)
1044  {
1045  std::stringstream error_stream;
1046  error_stream
1047  << "Helmholtz elements only store two fields so fld = "
1048  << fld << " is illegal\n";
1049  throw OomphLibError(
1050  error_stream.str(),
1051  OOMPH_CURRENT_FUNCTION,
1052  OOMPH_EXCEPTION_LOCATION);
1053  }
1054 #endif
1055  //Find the index at which the variable is stored
1056  std::complex<unsigned> complex_u_nodal_index =
1057  this->u_index_pml_fourier_decomposed_helmholtz();
1058  unsigned u_nodal_index = 0;
1059  if (fld==0)
1060  {
1061  u_nodal_index = complex_u_nodal_index.real();
1062  }
1063  else
1064  {
1065  u_nodal_index = complex_u_nodal_index.imag();
1066  }
1067 
1068 
1069  //Local shape function
1070  unsigned n_node=this->nnode();
1071  Shape psi(n_node);
1072 
1073  //Find values of shape function
1074  this->shape(s,psi);
1075 
1076  //Initialise value of u
1077  double interpolated_u = 0.0;
1078 
1079  //Sum over the local nodes
1080  for(unsigned l=0;l<n_node;l++)
1081  {
1082  interpolated_u += this->nodal_value(t,l,u_nodal_index)*psi[l];
1083  }
1084  return interpolated_u;
1085  }
1086 
1087 
1088 
1089 
1090  ///Return number of values in field fld: One per node
1091  unsigned nvalue_of_field(const unsigned &fld)
1092  {
1093 #ifdef PARANOID
1094  if (fld>1)
1095  {
1096  std::stringstream error_stream;
1097  error_stream
1098  << "Helmholtz elements only store two fields so fld = "
1099  << fld << " is illegal\n";
1100  throw OomphLibError(
1101  error_stream.str(),
1102  OOMPH_CURRENT_FUNCTION,
1103  OOMPH_EXCEPTION_LOCATION);
1104  }
1105 #endif
1106  return this->nnode();
1107  }
1108 
1109 
1110  ///Return local equation number of value j in field fld.
1111  int local_equation(const unsigned &fld,
1112  const unsigned &j)
1113  {
1114 #ifdef PARANOID
1115  if (fld>1)
1116  {
1117  std::stringstream error_stream;
1118  error_stream
1119  << "Helmholtz elements only store two fields so fld = "
1120  << fld << " is illegal\n";
1121  throw OomphLibError(
1122  error_stream.str(),
1123  OOMPH_CURRENT_FUNCTION,
1124  OOMPH_EXCEPTION_LOCATION);
1125  }
1126 #endif
1127  std::complex<unsigned> complex_u_nodal_index =
1128  this->u_index_pml_fourier_decomposed_helmholtz();
1129  unsigned u_nodal_index = 0;
1130  if (fld==0)
1131  {
1132  u_nodal_index = complex_u_nodal_index.real();
1133  }
1134  else
1135  {
1136  u_nodal_index = complex_u_nodal_index.imag();
1137  }
1138  return this->nodal_local_eqn(j,u_nodal_index);
1139  }
1140 
1141 
1142 
1143 
1144  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
1145  /// n_plot^DIM plot points
1146  void output(std::ostream &outfile, const unsigned &nplot)
1147  {
1149  }
1150 
1151 
1152  };
1153 
1154 
1155 //=======================================================================
1156 /// Face geometry for element is the same as that for the underlying
1157 /// wrapped element
1158 //=======================================================================
1159  template<class ELEMENT>
1162  : public virtual FaceGeometry<ELEMENT>
1163  {
1164  public:
1165  FaceGeometry() : FaceGeometry<ELEMENT>() {}
1166  };
1167 
1168 
1169 //=======================================================================
1170 /// Face geometry of the Face Geometry for element is the same as
1171 /// that for the underlying wrapped element
1172 //=======================================================================
1173  template<class ELEMENT>
1176  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
1177  {
1178  public:
1180  };
1181 
1182 
1183 
1184 ////////////////////////////////////////////////////////////////////////
1185 ////////////////////////////////////////////////////////////////////////
1186 ////////////////////////////////////////////////////////////////////////
1187 
1188 //=======================================================================
1189 /// Policy class defining the elements to be used in the actual
1190 /// PML layers. Same!
1191 //=======================================================================
1192  template<unsigned NNODE_1D>
1195  public virtual QPMLFourierDecomposedHelmholtzElement<NNODE_1D>
1196 {
1197 
1198  public:
1199 
1200  /// \short Constructor: Call the constructor for the
1201  /// appropriate QElement
1204  {}
1205 
1206 };
1207 
1208 }
1209 
1210 #endif
void compute_norm(double &norm)
Compute norm of fe solution.
PMLFourierDecomposedHelmholtzSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
QPMLFourierDecomposedHelmholtzElement(const QPMLFourierDecomposedHelmholtzElement< NNODE_1D > &dummy)
Broken copy constructor.
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!) ...
std::complex< double > gamma(const double &nu_i, const double &pml_width_i, const double &k_squared)
Overwrite the pure PML mapping coefficient function to return the mapping function proposed by Bermud...
PMLMappingAndTransformedCoordinate * Pml_mapping_and_transformed_coordinate_pt
Pointer to class which holds the pml mapping function (also known as gamma) and the associated transf...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
double factorial(const unsigned &l)
Factorial.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
Fourier decomposed Helmholtz upgraded to become projectable.
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
void output_real_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for real part of full time-dependent fct u = Re( (u_r +i u_i) exp(-i omega t) at phas...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
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
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) at...
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
A general Finite Element class.
Definition: elements.h:1271
void output_real_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for real part of full time-dependent fct u = Re( (u_r +i u_i) exp(-i omega t) at phas...
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and element Jacobian matrix (wrapper) ...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
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.
PMLFourierDecomposedHelmholtzSourceFctPt Source_fct_pt
Pointer to source function:
double plgndr1(const unsigned &n, const double &x)
Legendre polynomials depending on one parameter.
std::complex< double > interpolated_u_pml_fourier_decomposed_helmholtz(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
void(* PMLFourierDecomposedHelmholtzSourceFctPt)(const Vector< double > &x, std::complex< double > &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
static BermudezPMLMappingAndTransformedCoordinate Default_pml_mapping_and_transformed_coordinate
PMLFourierDecomposedHelmholtzSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
virtual double dshape_and_dtest_eulerian_at_knot_pml_fourier_decomposed_helmholtz(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
QPMLFourierDecomposedHelmholtzElement()
Constructor: Call constructors for QElement and PMLFourierDecomposedHelmholtz equations.
Base class for elements with pml capabilities.
Definition: pml_meshes.h:65
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...
PMLMappingAndTransformedCoordinate *& pml_mapping_and_transformed_coordinate_pt()
Return a pointer to the PML Mapping object.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: r,z,u_exact at n_plot^2 plot points.
static char t char * s
Definition: cfortran.h:572
PMLMappingAndTransformedCoordinate *const & pml_mapping_and_transformed_coordinate_pt() const
Return a pointer to the PML Mapping object (const version)
void output(FILE *file_pt)
C-style output function: r,z,u.
void output(FILE *file_pt)
C_style output with default number of plot points.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
bool Pml_is_enabled
Boolean indicating if element is used in pml mode.
Definition: pml_meshes.h:125
std::complex< double > transformed_coordinate(const double &nu_i, const double &pml_width_i, const double &pml_inner_boundary, const double &k_squared)
Overwrite the pure PML mapping coefficient function to return the transformed coordinate proposed by ...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
virtual void get_source_pml_fourier_decomposed_helmholtz(const unsigned &ipt, const Vector< double > &x, std::complex< double > &source) const
std::vector< bool > Pml_direction_active
Coordinate direction along which pml boundary is constant; alternatively: coordinate direction in whi...
Definition: pml_meshes.h:130
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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 with default number of plot points.
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 output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
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
double dshape_and_dtest_eulerian_at_knot_pml_fourier_decomposed_helmholtz(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
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
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) at...
virtual void fill_in_generic_residual_contribution_pml_fourier_decomposed_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
virtual std::complex< unsigned > u_index_pml_fourier_decomposed_helmholtz() const
Broken assignment operator.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
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 ...
double dshape_and_dtest_eulerian_pml_fourier_decomposed_helmholtz(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
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
unsigned nfields_for_projection()
Number of fields to be projected: 2 (real and imag part)
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
virtual std::complex< double > gamma(const double &nu_i, const double &pml_width_i, const double &k_squared)=0
Pure virtual to return PML mapping gamma, where gamma is the as function of where where h is the v...
virtual std::complex< double > transformed_coordinate(const double &nu_i, const double &pml_width_i, const double &pml_inner_boundary, const double &k_squared)=0
Pure virtual to return PML transformed coordinate, also known as as function of where where h is t...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
void output(std::ostream &outfile)
Output with default number of plot points.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
void output(std::ostream &outfile)
Output function: r,z,u.
virtual double dshape_and_dtest_eulerian_pml_fourier_decomposed_helmholtz(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
int *& pml_fourier_wavenumber_pt()
Get pointer to Fourier wavenumber.
void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. r,z,u_exact at n_plot^2 plot points (Calls the s...
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: (dummy time-dependent version to keep intel compiler happy)
PMLFourierDecomposedHelmholtzEquations(const PMLFourierDecomposedHelmholtzEquations &dummy)
Broken copy constructor.
void compute_pml_coefficients(const unsigned &ipt, const Vector< double > &x, Vector< std::complex< double > > &pml_laplace_factor, std::complex< double > &pml_k_squared_factor)
Compute pml coefficients at position x and integration point ipt. pml_laplace_factor is used in the r...
static double Default_Physical_Constant_Value
Static default value for the physical constants (initialised to zero)
ProjectablePMLFourierDecomposedHelmholtzElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
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...
void get_flux(const Vector< double > &s, Vector< std::complex< double > > &flux) const
Get flux: flux[i] = du/dx_i for real and imag part.
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
double plgndr2(const unsigned &l, const unsigned &m, const double &x)
Legendre polynomials depending on two parameters.
void output(std::ostream &outfile, const unsigned &nplot)
Output FE representation of soln: x,y,u or x,y,z,u at n_plot^DIM plot points.
void compute_complex_r(const unsigned &ipt, const Vector< double > &x, std::complex< double > &complex_r)
Compute complex variable r at position x[0] and integration point ipt.