pml_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 Helmholtz elements
31 #ifndef OOMPH_PML_HELMHOLTZ_ELEMENTS_HEADER
32 #define OOMPH_PML_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 //OOMPH-LIB headers
44 #include "../generic/nodes.h"
45 #include "../generic/Qelements.h"
46 #include "../generic/oomph_utilities.h"
47 #include "../generic/pml_meshes.h"
48 #include "../generic/projection.h"
49 #include "../generic/pml_mapping_functions.h"
50 
51 ////////////////////////////////////////////////////////////////////
52 ////////////////////////////////////////////////////////////////////
53 ////////////////////////////////////////////////////////////////////
54 
55 namespace oomph
56 {
57 
58 //=============================================================
59 /// A class for all isoparametric elements that solve the
60 /// Helmholtz equations with pml capabilities.
61 /// This contains the generic maths. Shape functions, geometric
62 /// mapping etc. must get implemented in derived class.
63 //=============================================================
64  template <unsigned DIM>
66  public virtual PMLElementBase<DIM>,
67  public virtual FiniteElement
68  {
69  public:
70 
71  /// \short Function pointer to source function fct(x,f(x)) --
72  /// x is a Vector!
73  typedef void (*PMLHelmholtzSourceFctPt)(const Vector<double>& x,
74  std::complex<double>& f);
75 
76  /// Constructor
78  K_squared_pt(0)
79  {
80  // Provide pointer to static method (Save memory)
83  }
84 
85 
86  /// Broken copy constructor
88  {
89  BrokenCopy::broken_copy("PMLHelmholtzEquations");
90  }
91 
92  /// Broken assignment operator
93 //Commented out broken assignment operator because this can lead to a conflict warning
94 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
95 //realise that two separate implementations of the broken function are the same and so,
96 //quite rightly, it shouts.
97  /*void operator=(const PMLHelmholtzEquations&)
98  {
99  BrokenCopy::broken_assign("PMLHelmholtzEquations");
100  }*/
101 
102  /// \short Return the index at which the unknown value
103  /// is stored.
104  virtual inline std::complex<unsigned> u_index_helmholtz() const
105  {return std::complex<unsigned>(0,1);}
106 
107  /// Get pointer to k_squared
108  double*& k_squared_pt(){ return K_squared_pt; }
109 
110 
111  /// Get the square of wavenumber
112  double k_squared()
113  {
114 #ifdef PARANOID
115  if (K_squared_pt==0)
116  {
117  throw OomphLibError(
118  "Please set pointer to k_squared using access fct to pointer!",
119  OOMPH_CURRENT_FUNCTION,
120  OOMPH_EXCEPTION_LOCATION);
121  }
122 #endif
123  return *K_squared_pt;
124  }
125 
126  /// Alpha, wavenumber complex shift
127  const double& alpha() const {return *Alpha_pt;}
128 
129 /// Pointer to Alpha, wavenumber complex shift
130  double* &alpha_pt() {return Alpha_pt;}
131 
132 
133  /// \short Number of scalars/fields output by this element. Reimplements
134  /// broken virtual function in base class.
135  unsigned nscalar_paraview() const
136  {
137  return 2;
138  }
139 
140  /// \short Write values of the i-th scalar field at the plot points. Needs
141  /// to be implemented for each new specific element type.
142  void scalar_value_paraview(std::ofstream& file_out,
143  const unsigned& i,
144  const unsigned& nplot) const
145  {
146 
147  //Vector of local coordinates
148  Vector<double> s(DIM);
149 
150  // Loop over plot points
151  unsigned num_plot_points=nplot_points_paraview(nplot);
152  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
153  {
154  // Get local coordinates of plot point
155  get_s_plot(iplot,nplot,s);
156  std::complex<double> u(interpolated_u_pml_helmholtz(s));
157 
158  // Paraview need to ouput the fields separately so it loops through all
159  // the elements twice
160  switch(i)
161  {
162  // Real part first
163  case 0:
164  file_out << u.real() << std::endl;
165  break;
166 
167  // Imaginary part second
168  case 1:
169  file_out << u.imag() << std::endl;
170  break;
171 
172  // Never get here
173  default:
174 #ifdef PARANOID
175  std::stringstream error_stream;
176  error_stream
177  << "PML Helmholtz elements only store 2 fields (real and imaginary) "
178  << "so i must be 0 or 1 rather "
179  << "than " << i << std::endl;
180  throw OomphLibError(
181  error_stream.str(),
182  OOMPH_CURRENT_FUNCTION,
183  OOMPH_EXCEPTION_LOCATION);
184 #endif
185  break;
186 
187  }
188  }
189  }
190 
191  /// \short Write values of the i-th scalar field at the plot points. Needs
192  /// to be implemented for each new specific element type.
193  void scalar_value_fct_paraview(std::ofstream& file_out,
194  const unsigned& i,
195  const unsigned& nplot,
197  exact_soln_pt) const
198  {
199  //Vector of local coordinates
200  Vector<double> s(DIM);
201 
202  // Vector for coordinates
203  Vector<double> x(DIM);
204 
205  // Exact solution Vector
206  Vector<double> exact_soln(2);
207 
208  // Loop over plot points
209  unsigned num_plot_points=nplot_points_paraview(nplot);
210  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
211  {
212  // Get local coordinates of plot point
213  get_s_plot(iplot,nplot,s);
214 
215  // Get x position as Vector
216  interpolated_x(s,x);
217 
218  // Get exact solution at this point
219  (*exact_soln_pt)(x,exact_soln);
220 
221  // Paraview need to output the fields separately so it loops through all
222  // the elements twice
223  switch(i)
224  {
225  // Real part first
226  case 0:
227  file_out << exact_soln[0] << std::endl;
228  break;
229 
230  // Imaginary part second
231  case 1:
232  file_out << exact_soln[1] << std::endl;
233  break;
234 
235  // Never get here
236  default:
237 #ifdef PARANOID
238  std::stringstream error_stream;
239  error_stream
240  << "PML Helmholtz elements only store 2 fields (real and imaginary) "
241  << "so i must be 0 or 1 rather "
242  << "than " << i << std::endl;
243  throw OomphLibError(
244  error_stream.str(),
245  OOMPH_CURRENT_FUNCTION,
246  OOMPH_EXCEPTION_LOCATION);
247 #endif
248  break;
249 
250  }
251  }
252  }
253 
254  /// \short Name of the i-th scalar field. Default implementation
255  /// returns V1 for the first one, V2 for the second etc. Can (should!)
256  /// be overloaded with more meaningful names in specific elements.
257  std::string scalar_name_paraview(const unsigned& i) const
258  {
259  switch(i)
260  {
261  case 0:
262  return "Real part";
263  break;
264 
265  case 1:
266  return "Imaginary part";
267  break;
268 
269  // Never get here
270  default:
271 #ifdef PARANOID
272  std::stringstream error_stream;
273  error_stream
274  << "PML Helmholtz elements only store 2 fields (real and imaginary) "
275  << "so i must be 0 or 1 rather "
276  << "than " << i << std::endl;
277  throw OomphLibError(
278  error_stream.str(),
279  OOMPH_CURRENT_FUNCTION,
280  OOMPH_EXCEPTION_LOCATION);
281 #endif
282 
283  // Dummy return for the default statement
284  return " ";
285  break;
286  }
287  }
288 
289  /// Output with default number of plot points
290  void output(std::ostream &outfile)
291  {
292  const unsigned n_plot=5;
293  output(outfile,n_plot);
294  }
295 
296  /// \short Output FE representation of soln: x,y,u_re,u_im or
297  /// x,y,z,u_re,u_im at n_plot^DIM plot points
298  void output(std::ostream &outfile, const unsigned &n_plot);
299 
300 
301 /// Output function for real part of full time-dependent solution
302 /// constructed by adding the scattered field
303 /// u = Re( (u_r +i u_i) exp(-i omega t)
304 /// at phase angle omega t = phi computed here, to the corresponding
305 /// incoming wave specified via the function pointer.
306  void output_total_real(
307  std::ostream &outfile,
308  FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt,
309  const double& phi,
310  const unsigned &nplot);
311 
312 
313  /// \short Output function for real part of full time-dependent solution
314  /// u = Re( (u_r +i u_i) exp(-i omega t))
315  /// at phase angle omega t = phi.
316  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
317  /// direction
318  void output_real(std::ostream &outfile, const double& phi,
319  const unsigned &n_plot);
320 
321  /// \short Output function for imaginary part of full time-dependent solution
322  /// u = Im( (u_r +i u_i) exp(-i omega t) )
323  /// at phase angle omega t = phi.
324  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
325  /// direction
326  void output_imag(std::ostream &outfile, const double& phi,
327  const unsigned &n_plot);
328 
329  /// C_style output with default number of plot points
330  void output(FILE* file_pt)
331  {
332  const unsigned n_plot=5;
333  output(file_pt,n_plot);
334  }
335 
336  /// \short C-style output FE representation of soln: x,y,u_re,u_im or
337  /// x,y,z,u_re,u_im at n_plot^DIM plot points
338  void output(FILE* file_pt, const unsigned &n_plot);
339 
340  /// Output exact soln: x,y,u_re_exact,u_im_exact
341  /// or x,y,z,u_re_exact,u_im_exact at n_plot^DIM plot points
342  void output_fct(std::ostream &outfile, const unsigned &n_plot,
344 
345  /// \short Output exact soln: (dummy time-dependent version to
346  /// keep intel compiler happy)
347  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
348  const double& time,
350  exact_soln_pt)
351  {
352  throw OomphLibError(
353  "There is no time-dependent output_fct() for PMLHelmholtz elements.",
354  OOMPH_CURRENT_FUNCTION,
355  OOMPH_EXCEPTION_LOCATION);
356  }
357 
358 
359 
360  /// \short Output function for real part of full time-dependent fct
361  /// u = Re( (u_r +i u_i) exp(-i omega t)
362  /// at phase angle omega t = phi.
363  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
364  /// direction
365  void output_real_fct(std::ostream &outfile,
366  const double& phi,
367  const unsigned &n_plot,
369 
370  /// \short Output function for imaginary part of full time-dependent fct
371  /// u = Im( (u_r +i u_i) exp(-i omega t))
372  /// at phase angle omega t = phi.
373  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
374  /// direction
375  void output_imag_fct(std::ostream &outfile,
376  const double& phi,
377  const unsigned &n_plot,
379 
380 
381  /// Get error against and norm of exact solution
382  void compute_error(std::ostream &outfile,
384  double& error, double& norm);
385 
386 
387  /// Dummy, time dependent error checker
388  void compute_error(std::ostream &outfile,
390  const double& time, double& error, double& norm)
391  {
392  throw OomphLibError(
393  "There is no time-dependent compute_error() for PMLHelmholtz elements.",
394  OOMPH_CURRENT_FUNCTION,
395  OOMPH_EXCEPTION_LOCATION);
396  }
397 
398 
399  /// Compute norm of fe solution
400  void compute_norm(double& norm);
401 
402  /// Access function: Pointer to source function
404 
405  /// Access function: Pointer to source function. Const version
407 
408 
409  /// Get source term at (Eulerian) position x. This function is
410  /// virtual to allow overloading in multi-physics problems where
411  /// the strength of the source function might be determined by
412  /// another system of equations.
413  inline virtual void get_source_helmholtz(const unsigned& ipt,
414  const Vector<double>& x,
415  std::complex<double>& source) const
416  {
417  //If no source function has been set, return zero
418  if(Source_fct_pt==0)
419  {
420  source = std::complex<double>(0.0,0.0);
421  }
422  else
423  {
424  // Get source strength
425  (*Source_fct_pt)(x,source);
426  }
427  }
428 
429  /// \short Pure virtual function in which we specify the
430  /// values to be pinned (and set to zero) on the outer edge of
431  /// the pml layer. All of them! Vector is resized internally.
433  values_to_pin)
434  {
435  values_to_pin.resize(2);
436 
437  values_to_pin[0]=0;
438  values_to_pin[1]=1;
439  }
440 
441 
442  /// Get flux: flux[i] = du/dx_i for real and imag part
443  void get_flux(const Vector<double>& s,
444  Vector<std::complex <double> >& flux) const
445  {
446  //Find out how many nodes there are in the element
447  const unsigned n_node = nnode();
448 
449 
450  //Set up memory for the shape and test functions
451  Shape psi(n_node);
452  DShape dpsidx(n_node,DIM);
453 
454  //Call the derivatives of the shape and test functions
455  dshape_eulerian(s,psi,dpsidx);
456 
457  //Initialise to zero
458  const std::complex<double> zero(0.0,0.0);
459  for(unsigned j=0;j<DIM;j++)
460  {
461  flux[j] = zero;
462  }
463 
464  // Loop over nodes
465  for(unsigned l=0;l<n_node;l++)
466  {
467  //Cache the complex value of the unknown
468  const std::complex<double> u_value(
469  this->nodal_value(l,u_index_helmholtz().real()),
470  this->nodal_value(l,u_index_helmholtz().imag()));
471  //Loop over derivative directions
472  for(unsigned j=0;j<DIM;j++)
473  {
474  flux[j] += u_value*dpsidx(l,j);
475  }
476  }
477  }
478 
479 
480  /// Add the element's contribution to its residual vector (wrapper)
482  {
483  //Call the generic residuals function with flag set to 0
484  //using a dummy matrix argument
487  }
488 
489 
490  /// \short Add the element's contribution to its residual vector and
491  /// element Jacobian matrix (wrapper)
493  DenseMatrix<double> &jacobian)
494  {
495  //Call the generic routine with the flag set to 1
497  }
498 
499 
500 
501  /// \short Return FE representation of function value u_helmholtz(s)
502  /// at local coordinate s
503  inline std::complex<double>
505  {
506  //Find number of nodes
507  const unsigned n_node = nnode();
508 
509  //Local shape function
510  Shape psi(n_node);
511 
512  //Find values of shape function
513  shape(s,psi);
514 
515  //Initialise value of u
516  std::complex<double> interpolated_u(0.0,0.0);
517 
518  //Get the index at which the helmholtz unknown is stored
519  const unsigned u_nodal_index_real = u_index_helmholtz().real();
520  const unsigned u_nodal_index_imag = u_index_helmholtz().imag();
521 
522  //Loop over the local nodes and sum
523  for(unsigned l=0;l<n_node;l++)
524  {
525  //Make a temporary complex number from the stored data
526  const std::complex<double> u_value(
527  this->nodal_value(l,u_nodal_index_real),
528  this->nodal_value(l,u_nodal_index_imag));
529  //Add to the interpolated value
530  interpolated_u += u_value*psi[l];
531  }
532  return interpolated_u;
533  }
534 
535 
536  /// \short Self-test: Return 0 for OK
537  unsigned self_test();
538 
539 
540 /// \short Compute pml coefficients at position x and integration point ipt.
541 /// pml_laplace_factor is used in the residual contribution from the laplace
542 /// operator, similarly pml_k_squared_factor is used in the contribution from
543 /// the k^2 of the Helmholtz operator.
545  const unsigned& ipt,
546  const Vector<double>& x,
547  Vector< std::complex<double> >& pml_laplace_factor,
548  std::complex<double>& pml_k_squared_factor)
549  {
550 
551  /// Vector which points from the inner boundary to x
552  Vector<double> nu(DIM);
553  for(unsigned k=0; k<DIM; k++)
554  {
555  nu[k] = x[k] - this->Pml_inner_boundary[k];
556  }
557 
558  /// Vector which points from the inner boundary to the edge of the boundary
559  Vector<double> pml_width(DIM);
560  for(unsigned k=0; k<DIM; k++)
561  {
562  pml_width[k] = this->Pml_outer_boundary[k] - this->Pml_inner_boundary[k];
563  }
564 
565  // Declare gamma_i vectors of complex numbers for PML weights
566  Vector<std::complex<double> > pml_gamma(DIM);
567 
568  if (this->Pml_is_enabled)
569  {
570  // Cache k_squared to pass into mapping function
571  double k_squared_local = k_squared();
572 
573  for(unsigned k=0; k<DIM; k++) {
574  // If PML is enabled in the respective direction
575  if (this->Pml_direction_active[k])
576  {
577  pml_gamma[k] = Pml_mapping_pt->gamma(nu[k],pml_width[k],
578  k_squared_local,alpha());
579  }
580  else
581  {
582  pml_gamma[k] = 1.0;
583  }
584  }
585 
586  /// \short for 2D, in order:
587  /// g_y/g_x, g_x/g_y for Laplace bit and g_x*g_y for Helmholtz bit
588  /// 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
589  /// and g_x*g_y*g_z for Helmholtz bit
590  if (DIM == 2)
591  {
592  pml_laplace_factor[0] = pml_gamma[1] / pml_gamma[0];
593  pml_laplace_factor[1] = pml_gamma[0] / pml_gamma[1];
594  pml_k_squared_factor = pml_gamma[0] * pml_gamma[1];
595  }
596  else // if (DIM == 3)
597  {
598  pml_laplace_factor[0] = pml_gamma[1] * pml_gamma[2] / pml_gamma[0];
599  pml_laplace_factor[1] = pml_gamma[0] * pml_gamma[2] / pml_gamma[1];
600  pml_laplace_factor[2] = pml_gamma[0] * pml_gamma[1] / pml_gamma[2];
601  pml_k_squared_factor = pml_gamma[0] * pml_gamma[1] * pml_gamma[2];
602  }
603 
604  }
605  else
606  {
607  /// \short The weights all default to 1.0 as if the propagation
608  /// medium is the physical domain
609  for(unsigned k=0; k<DIM; k++)
610  {
611  pml_laplace_factor[k] = std::complex<double> (1.0,0.0);
612  }
613 
614  pml_k_squared_factor = std::complex<double> (1.0,0.0);
615  }
616 
617  }
618 
619  /// Return a pointer to the PML Mapping object
621 
622  /// Return a pointer to the PML Mapping object (const version)
623  PMLMapping* const &pml_mapping_pt() const {return Pml_mapping_pt;}
624 
625  /// Static so that the class doesn't need to instantiate a new default
626  /// everytime it uses it
628 
629  /// \short The number of "DOF types" that degrees of freedom in this element
630  /// are sub-divided into: real and imaginary part
631  unsigned ndof_types() const
632  {
633  return 2;
634  }
635 
636  /// \short Create a list of pairs for all unknowns in this element,
637  /// so that the first entry in each pair contains the global equation
638  /// number of the unknown, while the second one contains the number
639  /// of the "DOF type" that this unknown is associated with.
640  /// (Function can obviously only be called if the equation numbering
641  /// scheme has been set up.) Real=0; Imag=1
643  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
644  {
645  // temporary pair (used to store dof lookup prior to being added to list)
646  std::pair<unsigned,unsigned> dof_lookup;
647 
648  // number of nodes
649  unsigned n_node = this->nnode();
650 
651  // loop over the nodes
652  for (unsigned n = 0; n < n_node; n++)
653  {
654  // determine local eqn number for real part
655  int local_eqn_number =
656  this->nodal_local_eqn(n, u_index_helmholtz().real());
657 
658  // ignore pinned values
659  if (local_eqn_number >= 0)
660  {
661  // store dof lookup in temporary pair: First entry in pair
662  // is global equation number; second entry is dof type
663  dof_lookup.first = this->eqn_number(local_eqn_number);
664  dof_lookup.second = 0;
665 
666  // add to list
667  dof_lookup_list.push_front(dof_lookup);
668  }
669 
670  // determine local eqn number for imag part
671  local_eqn_number =
672  this->nodal_local_eqn(n, u_index_helmholtz().imag());
673 
674  // ignore pinned values
675  if (local_eqn_number >= 0)
676  {
677  // store dof lookup in temporary pair: First entry in pair
678  // is global equation number; second entry is dof type
679  dof_lookup.first = this->eqn_number(local_eqn_number);
680  dof_lookup.second = 1;
681 
682  // add to list
683  dof_lookup_list.push_front(dof_lookup);
684  }
685  }
686  }
687 
688 
689  protected:
690 
691  /// \short Shape/test functions and derivs w.r.t. to global coords at
692  /// local coord. s; return Jacobian of mapping
694  Shape &psi,
695  DShape &dpsidx, Shape &test,
696  DShape &dtestdx) const=0;
697 
698 
699  /// \short Shape/test functions and derivs w.r.t. to global coords at
700  /// integration point ipt; return Jacobian of mapping
702  const unsigned &ipt,
703  Shape &psi,
704  DShape &dpsidx,
705  Shape &test,
706  DShape &dtestdx) const=0;
707 
708  /// \short Compute element residual Vector only (if flag=and/or element
709  /// Jacobian matrix
711  Vector<double> &residuals, DenseMatrix<double> &jacobian,
712  const unsigned& flag);
713 
714  /// Pointer to wavenumber complex shift
715  double *Alpha_pt;
716 
717  /// Pointer to source function:
719 
720  /// Pointer to wave number (must be set!)
721  double* K_squared_pt;
722 
723  /// Pointer to class which holds the pml mapping function, also known as gamma
725 
726  /// Static default value for the physical constants (initialised to zero)
728 
729  };
730 
731 
732 ///////////////////////////////////////////////////////////////////////////
733 ///////////////////////////////////////////////////////////////////////////
734 ///////////////////////////////////////////////////////////////////////////
735 
736 
737 
738 //======================================================================
739 /// \short QPMLHelmholtzElement elements are linear/quadrilateral/
740 /// brick-shaped PMLHelmholtz elements with isoparametric
741 /// interpolation for the function.
742 //======================================================================
743  template <unsigned DIM, unsigned NNODE_1D>
744  class QPMLHelmholtzElement : public virtual QElement<DIM,NNODE_1D>,
745  public virtual PMLHelmholtzEquations<DIM>
746  {
747 
748  private:
749 
750  /// \short Static int that holds the number of variables at
751  /// nodes: always the same
752  static const unsigned Initial_Nvalue;
753 
754  public:
755 
756 
757  ///\short Constructor: Call constructors for QElement and
758  /// PMLHelmholtz equations
760  QElement<DIM,NNODE_1D>(), PMLHelmholtzEquations<DIM>()
761  {}
762 
763  /// Broken copy constructor
766  {
767  BrokenCopy::broken_copy("QPMLHelmholtzElement");
768  }
769 
770  /// Broken assignment operator
771  /*void operator=(const QPMLHelmholtzElement<DIM,NNODE_1D>&)
772  {
773  BrokenCopy::broken_assign("QPMLHelmholtzElement");
774  }*/
775 
776 
777  /// \short Required # of `values' (pinned or dofs)
778  /// at node n
779  inline unsigned required_nvalue(const unsigned &n) const
780  {return Initial_Nvalue;}
781 
782  /// \short Output function:
783  /// x,y,u or x,y,z,u
784  void output(std::ostream &outfile)
786 
787 
788  /// \short Output function:
789  /// x,y,u or x,y,z,u at n_plot^DIM plot points
790  void output(std::ostream &outfile, const unsigned &n_plot)
791  {PMLHelmholtzEquations<DIM>::output(outfile,n_plot);}
792 
793  /// \short Output function for real part of full time-dependent solution
794  /// u = Re( (u_r +i u_i) exp(-i omega t)
795  /// at phase angle omega t = phi.
796  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
797  /// direction
798  void output_real(std::ostream &outfile, const double& phi,
799  const unsigned &n_plot)
800  {PMLHelmholtzEquations<DIM>::output_real(outfile,phi,n_plot);}
801 
802  /// \short Output function for imaginary part of full time-dependent solution
803  /// u = Im( (u_r +i u_i) exp(-i omega t))
804  /// at phase angle omega t = phi.
805  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
806  /// direction
807  void output_imag(std::ostream &outfile, const double& phi,
808  const unsigned &n_plot)
809  {PMLHelmholtzEquations<DIM>::output_imag(outfile,phi,n_plot);}
810 
811 
812  /// \short C-style output function:
813  /// x,y,u or x,y,z,u
814  void output(FILE* file_pt)
816 
817 
818  /// \short C-style output function:
819  /// x,y,u or x,y,z,u at n_plot^DIM plot points
820  void output(FILE* file_pt, const unsigned &n_plot)
821  {PMLHelmholtzEquations<DIM>::output(file_pt,n_plot);}
822 
823 
824  /// \short Output function for an exact solution:
825  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
826  void output_fct(std::ostream &outfile, const unsigned &n_plot,
829  n_plot,
830  exact_soln_pt);}
831 
832 
833  /// \short Output function for real part of full time-dependent fct
834  /// u = Re( (u_r +i u_i) exp(-i omega t)
835  /// at phase angle omega t = phi.
836  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
837  /// direction
838  void output_real_fct(std::ostream &outfile,
839  const double& phi,
840  const unsigned &n_plot,
842  {
844  phi,
845  n_plot,
846  exact_soln_pt);
847  }
848 
849  /// \short Output function for imaginary part of full time-dependent fct
850  /// u = Im( (u_r +i u_i) exp(-i omega t))
851  /// at phase angle omega t = phi.
852  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
853  /// direction
854  void output_imag_fct(std::ostream &outfile,
855  const double& phi,
856  const unsigned &n_plot,
858  {
860  phi,
861  n_plot,
862  exact_soln_pt);
863  }
864 
865 
866  /// \short Output function for a time-dependent exact solution.
867  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
868  /// (Calls the steady version)
869  void output_fct(std::ostream &outfile, const unsigned &n_plot,
870  const double& time,
873  n_plot,
874  time,
875  exact_soln_pt);}
876 
877  protected:
878 
879 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
881  const Vector<double> &s, Shape &psi, DShape &dpsidx,
882  Shape &test, DShape &dtestdx) const;
883 
884 
885  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
886  /// integration point ipt. Return Jacobian.
887  inline double dshape_and_dtest_eulerian_at_knot_helmholtz(const unsigned& ipt,
888  Shape &psi,
889  DShape &dpsidx,
890  Shape &test,
891  DShape &dtestdx)
892  const;
893 
894  };
895 
896 
897 
898 
899 //Inline functions:
900 
901 
902 //======================================================================
903 /// Define the shape functions and test functions and derivatives
904 /// w.r.t. global coordinates and return Jacobian of mapping.
905 ///
906 /// Galerkin: Test functions = shape functions
907 //======================================================================
908  template<unsigned DIM, unsigned NNODE_1D>
910  const Vector<double> &s,
911  Shape &psi,
912  DShape &dpsidx,
913  Shape &test,
914  DShape &dtestdx) const
915  {
916  //Call the geometrical shape functions and derivatives
917  const double J = this->dshape_eulerian(s,psi,dpsidx);
918 
919  //Set the test functions equal to the shape functions
920  test = psi;
921  dtestdx= dpsidx;
922 
923  //Return the jacobian
924  return J;
925  }
926 
927 
928 
929 
930 //======================================================================
931 /// Define the shape functions and test functions and derivatives
932 /// w.r.t. global coordinates and return Jacobian of mapping.
933 ///
934 /// Galerkin: Test functions = shape functions
935 //======================================================================
936  template<unsigned DIM, unsigned NNODE_1D>
939  const unsigned &ipt,
940  Shape &psi,
941  DShape &dpsidx,
942  Shape &test,
943  DShape &dtestdx) const
944  {
945  //Call the geometrical shape functions and derivatives
946  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
947 
948  //Set the pointers of the test functions
949  test = psi;
950  dtestdx = dpsidx;
951 
952  //Return the jacobian
953  return J;
954  }
955 
956 ////////////////////////////////////////////////////////////////////////
957 ////////////////////////////////////////////////////////////////////////
958 ////////////////////////////////////////////////////////////////////////
959 
960 
961 
962 //=======================================================================
963 /// Face geometry for the QPMLHelmholtzElement elements: The spatial
964 /// dimension of the face elements is one lower than that of the
965 /// bulk element but they have the same number of points
966 /// along their 1D edges.
967 //=======================================================================
968  template<unsigned DIM, unsigned NNODE_1D>
969  class FaceGeometry<QPMLHelmholtzElement<DIM,NNODE_1D> >:
970  public virtual QElement<DIM-1,NNODE_1D>
971  {
972 
973  public:
974 
975  /// \short Constructor: Call the constructor for the
976  /// appropriate lower-dimensional QElement
977  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
978 
979  };
980 
981 ////////////////////////////////////////////////////////////////////////
982 ////////////////////////////////////////////////////////////////////////
983 ////////////////////////////////////////////////////////////////////////
984 
985 
986 //=======================================================================
987 /// Face geometry for the 1D QPMLHelmholtzElement elements:
988 /// Point elements
989 //=======================================================================
990  template<unsigned NNODE_1D>
991  class FaceGeometry<QPMLHelmholtzElement<1,NNODE_1D> >:
992  public virtual PointElement
993  {
994 
995  public:
996 
997  /// \short Constructor: Call the constructor for the
998  /// appropriate lower-dimensional QElement
1000 
1001  };
1002 
1003 
1004 ////////////////////////////////////////////////////////////////////////
1005 ////////////////////////////////////////////////////////////////////////
1006 ////////////////////////////////////////////////////////////////////////
1007 
1008 
1009 //==========================================================
1010 /// PMLHelmholtz upgraded to become projectable
1011 //==========================================================
1012  template<class HELMHOLTZ_ELEMENT>
1014  public virtual ProjectableElement<HELMHOLTZ_ELEMENT>
1015  {
1016 
1017  public:
1018 
1019 
1020  /// \short Constructor [this was only required explicitly
1021  /// from gcc 4.5.2 onwards...]
1023 
1024  /// \short Specify the values associated with field fld.
1025  /// The information is returned in a vector of pairs which comprise
1026  /// the Data object and the value within it, that correspond to field fld.
1028  {
1029 
1030 #ifdef PARANOID
1031  if (fld>1)
1032  {
1033  std::stringstream error_stream;
1034  error_stream
1035  << "PMLHelmholtz elements only store 2 fields so fld = "
1036  << fld << " is illegal \n";
1037  throw OomphLibError(
1038  error_stream.str(),
1039  OOMPH_CURRENT_FUNCTION,
1040  OOMPH_EXCEPTION_LOCATION);
1041  }
1042 #endif
1043 
1044  // Create the vector
1045  unsigned nnod=this->nnode();
1046  Vector<std::pair<Data*,unsigned> > data_values(nnod);
1047 
1048  // Loop over all nodes
1049  for (unsigned j=0;j<nnod;j++)
1050  {
1051  // Add the data value associated field: The node itself
1052  data_values[j]=std::make_pair(this->node_pt(j),fld);
1053  }
1054 
1055  // Return the vector
1056  return data_values;
1057  }
1058 
1059  /// \short Number of fields to be projected: 2 (real and imag part)
1061  {
1062  return 2;
1063  }
1064 
1065  /// \short Number of history values to be stored for fld-th field.
1066  /// (Note: count includes current value!)
1067  unsigned nhistory_values_for_projection(const unsigned &fld)
1068  {
1069 #ifdef PARANOID
1070  if (fld>1)
1071  {
1072  std::stringstream error_stream;
1073  error_stream
1074  << "PMLHelmholtz elements only store two fields so fld = "
1075  << fld << " is illegal\n";
1076  throw OomphLibError(
1077  error_stream.str(),
1078  OOMPH_CURRENT_FUNCTION,
1079  OOMPH_EXCEPTION_LOCATION);
1080  }
1081 #endif
1082  return this->node_pt(0)->ntstorage();
1083  }
1084 
1085  ///\short Number of positional history values
1086  /// (Note: count includes current value!)
1088  {
1089  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
1090  }
1091 
1092  /// \short Return Jacobian of mapping and shape functions of field fld
1093  /// at local coordinate s
1094  double jacobian_and_shape_of_field(const unsigned &fld,
1095  const Vector<double> &s,
1096  Shape &psi)
1097  {
1098 #ifdef PARANOID
1099  if (fld>1)
1100  {
1101  std::stringstream error_stream;
1102  error_stream
1103  << "PMLHelmholtz elements only store two fields so fld = "
1104  << fld << " is illegal.\n";
1105  throw OomphLibError(
1106  error_stream.str(),
1107  OOMPH_CURRENT_FUNCTION,
1108  OOMPH_EXCEPTION_LOCATION);
1109  }
1110 #endif
1111  unsigned n_dim=this->dim();
1112  unsigned n_node=this->nnode();
1113  Shape test(n_node);
1114  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
1115  double J=this->dshape_and_dtest_eulerian_helmholtz(s,psi,dpsidx,
1116  test,dtestdx);
1117  return J;
1118  }
1119 
1120 
1121 
1122  /// \short Return interpolated field fld at local coordinate s, at time level
1123  /// t (t=0: present; t>0: history values)
1124  double get_field(const unsigned &t,
1125  const unsigned &fld,
1126  const Vector<double>& s)
1127  {
1128 #ifdef PARANOID
1129  if (fld>1)
1130  {
1131  std::stringstream error_stream;
1132  error_stream
1133  << "PMLHelmholtz elements only store two fields so fld = "
1134  << fld << " is illegal\n";
1135  throw OomphLibError(
1136  error_stream.str(),
1137  OOMPH_CURRENT_FUNCTION,
1138  OOMPH_EXCEPTION_LOCATION);
1139  }
1140 #endif
1141  //Find the index at which the variable is stored
1142  std::complex<unsigned> complex_u_nodal_index = this->u_index_helmholtz();
1143  unsigned u_nodal_index = 0;
1144  if (fld==0)
1145  {
1146  u_nodal_index = complex_u_nodal_index.real();
1147  }
1148  else
1149  {
1150  u_nodal_index = complex_u_nodal_index.imag();
1151  }
1152 
1153 
1154  //Local shape function
1155  unsigned n_node=this->nnode();
1156  Shape psi(n_node);
1157 
1158  //Find values of shape function
1159  this->shape(s,psi);
1160 
1161  //Initialise value of u
1162  double interpolated_u = 0.0;
1163 
1164  //Sum over the local nodes
1165  for(unsigned l=0;l<n_node;l++)
1166  {
1167  interpolated_u += this->nodal_value(t,l,u_nodal_index)*psi[l];
1168  }
1169  return interpolated_u;
1170  }
1171 
1172 
1173 
1174 
1175  ///Return number of values in field fld: One per node
1176  unsigned nvalue_of_field(const unsigned &fld)
1177  {
1178 #ifdef PARANOID
1179  if (fld>1)
1180  {
1181  std::stringstream error_stream;
1182  error_stream
1183  << "PMLHelmholtz elements only store two fields so fld = "
1184  << fld << " is illegal\n";
1185  throw OomphLibError(
1186  error_stream.str(),
1187  OOMPH_CURRENT_FUNCTION,
1188  OOMPH_EXCEPTION_LOCATION);
1189  }
1190 #endif
1191  return this->nnode();
1192  }
1193 
1194 
1195  ///Return local equation number of value j in field fld.
1196  int local_equation(const unsigned &fld,
1197  const unsigned &j)
1198  {
1199 #ifdef PARANOID
1200  if (fld>1)
1201  {
1202  std::stringstream error_stream;
1203  error_stream
1204  << "PMLHelmholtz elements only store two fields so fld = "
1205  << fld << " is illegal\n";
1206  throw OomphLibError(
1207  error_stream.str(),
1208  OOMPH_CURRENT_FUNCTION,
1209  OOMPH_EXCEPTION_LOCATION);
1210  }
1211 #endif
1212  std::complex<unsigned> complex_u_nodal_index = this->u_index_helmholtz();
1213  unsigned u_nodal_index = 0;
1214  if (fld==0)
1215  {
1216  u_nodal_index = complex_u_nodal_index.real();
1217  }
1218  else
1219  {
1220  u_nodal_index = complex_u_nodal_index.imag();
1221  }
1222  return this->nodal_local_eqn(j,u_nodal_index);
1223  }
1224 
1225  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
1226  /// n_plot^DIM plot points
1227  void output(std::ostream &outfile, const unsigned &nplot)
1228  {
1229  HELMHOLTZ_ELEMENT::output(outfile,nplot);
1230  }
1231 
1232 
1233  };
1234 
1235 
1236 //=======================================================================
1237 /// Face geometry for element is the same as that for the underlying
1238 /// wrapped element
1239 //=======================================================================
1240  template<class ELEMENT>
1242  : public virtual FaceGeometry<ELEMENT>
1243  {
1244  public:
1245  FaceGeometry() : FaceGeometry<ELEMENT>() {}
1246  };
1247 
1248 ////////////////////////////////////////////////////////////////////////
1249 ////////////////////////////////////////////////////////////////////////
1250 ////////////////////////////////////////////////////////////////////////
1251 
1252 //=======================================================================
1253 /// Policy class defining the elements to be used in the actual
1254 /// PML layers. Same!
1255 //=======================================================================
1256  template<unsigned NNODE_1D>
1258  QPMLHelmholtzElement<2,NNODE_1D> > :
1259  public virtual QPMLHelmholtzElement<2,NNODE_1D>
1260  {
1261 
1262  public:
1263 
1264  /// \short Constructor: Call the constructor for the
1265  /// appropriate QElement
1267  QPMLHelmholtzElement<2,NNODE_1D>()
1268  {}
1269 
1270  };
1271 
1272 } // End of Namespace oomph
1273 
1274 #endif
PMLHelmholtzSourceFctPt Source_fct_pt
Pointer to source function:
void output_total_real(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt, const double &phi, const unsigned &nplot)
virtual void get_source_helmholtz(const unsigned &ipt, const Vector< double > &x, std::complex< double > &source) const
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
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.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
const double & alpha() const
Alpha, wavenumber complex shift.
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
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.
PMLHelmholtzSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
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
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!) ...
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
PMLHelmholtz upgraded to become projectable.
cstr elem_len * i
Definition: cfortran.h:607
void output(std::ostream &outfile)
Output with default number of plot points.
PMLHelmholtzSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number...
Definition: elements.h:731
double dshape_and_dtest_eulerian_at_knot_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.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
ProjectablePMLHelmholtzElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
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
virtual double dshape_and_dtest_eulerian_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...
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
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
A general Finite Element class.
Definition: elements.h:1271
double * Alpha_pt
Pointer to wavenumber complex shift.
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
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(FILE *file_pt)
C_style output with default number of plot points.
PMLMapping * Pml_mapping_pt
Pointer to class which holds the pml mapping function, also known as gamma.
PMLHelmholtzEquations(const PMLHelmholtzEquations &dummy)
Broken copy constructor.
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
double k_squared()
Get the square of wavenumber.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
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...
static double Default_Physical_Constant_Value
Static default value for the physical constants (initialised to zero)
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
double dshape_and_dtest_eulerian_helmholtz(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
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) ...
void output_imag_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for imaginary part of full time-dependent fct u = Im( (u_r +i u_i) exp(-i omega t)) a...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
double *& k_squared_pt()
Get pointer to k_squared.
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.
Base class for elements with pml capabilities.
Definition: pml_meshes.h:65
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM 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 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...
std::complex< double > interpolated_u_pml_helmholtz(const Vector< double > &s) const
Return FE representation of function value u_helmholtz(s) at local coordinate s.
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...
static char t char * s
Definition: cfortran.h:572
PMLMapping *& pml_mapping_pt()
Return a pointer to the PML Mapping object.
bool Pml_is_enabled
Boolean indicating if element is used in pml mode.
Definition: pml_meshes.h:125
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
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, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot ...
std::vector< bool > Pml_direction_active
Coordinate direction along which pml boundary is constant; alternatively: coordinate direction in whi...
Definition: pml_meshes.h:130
unsigned nfields_for_projection()
Number of fields to be projected: 2 (real and imag part)
double * K_squared_pt
Pointer to wave number (must be set!)
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void compute_norm(double &norm)
Compute norm of fe solution.
void output(std::ostream &outfile)
Output with default number of plot points.
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
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 ...
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
static BermudezPMLMapping Default_pml_mapping
double *& alpha_pt()
Pointer to Alpha, wavenumber complex shift.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
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 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
virtual void fill_in_generic_residual_contribution_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points...
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
PMLMapping *const & pml_mapping_pt() const
Return a pointer to the PML Mapping object (const version)
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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...
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)
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
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
void output_imag_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for imaginary part of full time-dependent fct u = Im( (u_r +i u_i) exp(-i omega t)) a...
void scalar_value_fct_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
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...
QPMLHelmholtzElement()
Constructor: Call constructors for QElement and PMLHelmholtz equations.
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 output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
QPMLHelmholtzElement elements are linear/quadrilateral/ brick-shaped PMLHelmholtz elements with isopa...
void(* PMLHelmholtzSourceFctPt)(const Vector< double > &x, std::complex< double > &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
virtual double dshape_and_dtest_eulerian_at_knot_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 ...
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: real and imag...
PMLLayerElement()
Constructor: Call the constructor for the appropriate QElement.
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
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 ...
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
unsigned self_test()
Self-test: Return 0 for OK.
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...