refineable_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: 1282 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-01-16 08:27:53 +0000 (Mon, 16 Jan 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 refineable QPMLHelmholtzElement elements
31 
32 #ifndef OOMPH_REFINEABLE_PML_HELMHOLTZ_ELEMENTS_HEADER
33 #define OOMPH_REFINEABLE_PML_HELMHOLTZ_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 
41 //oomph-lib headers
42 #include "../generic/refineable_quad_element.h"
43 #include "../generic/refineable_brick_element.h"
44 #include "../generic/hp_refineable_elements.h"
45 #include "../generic/error_estimator.h"
46 #include "pml_helmholtz_elements.h"
47 
48 namespace oomph
49 {
50 
51 ///////////////////////////////////////////////////////////////////////////
52 ///////////////////////////////////////////////////////////////////////////
53 
54 
55 
56 //======================================================================
57 /// Refineable version of PMLHelmholtz equations
58 ///
59 ///
60 //======================================================================
61 template <unsigned DIM>
63  public virtual PMLHelmholtzEquations<DIM>,
64  public virtual RefineableElement,
65  public virtual ElementWithZ2ErrorEstimator
66  {
67  public:
68 
69  /// \short Constructor, simply call other constructors
72  { }
73 
74  /// Broken copy constructor
76  {
77  BrokenCopy::broken_copy("RefineablePMLHelmholtzEquations");
78  }
79 
80  /// Broken assignment operator
81 //Commented out broken assignment operator because this can lead to a conflict warning
82 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
83 //realise that two separate implementations of the broken function are the same and so,
84 //quite rightly, it shouts.
85  /*void operator=(const RefineablePMLHelmholtzEquations<DIM>&)
86  {
87  BrokenCopy::broken_assign("RefineablePMLHelmholtzEquations");
88  }*/
89 
90  /// Number of 'flux' terms for Z2 error estimation
91  unsigned num_Z2_flux_terms() {return 2*DIM;}
92 
93  /// \short Get 'flux' for Z2 error recovery: Complex flux from
94  /// PMLHelmholtz equations, strung together
96  {
97  Vector<std::complex <double> > complex_flux(DIM);
98  this->get_flux(s,complex_flux);
99  unsigned count=0;
100  for (unsigned i=0;i<DIM;i++)
101  {
102  flux[count]=complex_flux[i].real();
103  count++;
104  flux[count]=complex_flux[i].imag();
105  count++;
106  }
107  }
108 
109 
110 /// \short Get the function value u in Vector.
111 /// Note: Given the generality of the interface (this function
112 /// is usually called from black-box documentation or interpolation routines),
113 /// the values Vector sets its own size in here.
115  {
116  // Set size of Vector: u
117  values.resize(2);
118 
119  //Find number of nodes
120  unsigned n_node = nnode();
121 
122  //Local shape function
123  Shape psi(n_node);
124 
125  //Find values of shape function
126  shape(s,psi);
127 
128  //Initialise value of u
129  values[0] = 0.0;
130  values[1] = 0.0;
131 
132  //Find the index at which the pml_helmholtz unknown is stored
133  std::complex<unsigned> u_nodal_index = this->u_index_helmholtz();
134 
135  //Loop over the local nodes and sum up the values
136  for(unsigned l=0;l<n_node;l++)
137  {
138  values[0] += this->nodal_value(l,u_nodal_index.real())*psi[l];
139  values[1] += this->nodal_value(l,u_nodal_index.imag())*psi[l];
140  }
141  }
142 
143 
144  /// \short Get the function value u in Vector.
145  /// Note: Given the generality of the interface (this function
146  /// is usually called from black-box documentation or interpolation routines),
147  /// the values Vector sets its own size in here.
148  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
149  Vector<double>& values)
150  {
151  if (t!=0)
152  {
153  std::string error_message =
154  "Time-dependent version of get_interpolated_values() ";
155  error_message += "not implemented for this element \n";
156  throw
157  OomphLibError(error_message,
158  "RefineablePMLHelmholtzEquations::get_interpolated_values()",
159  OOMPH_EXCEPTION_LOCATION);
160  }
161  else
162  {
163  //Make sure that we call this particular object's steady
164  //get_interpolated_values (it could get overloaded lower down)
166  }
167  }
168 
169 
170  /// Further build: Copy source function pointer from father element
172  {
174  this->father_element_pt())->source_fct_pt();
175  }
176 
177 
178  private:
179 
180 
181 /// \short Add element's contribution to elemental residual vector and/or
182 /// Jacobian matrix
183 /// flag=1: compute both
184 /// flag=0: compute only residual vector
186  Vector<double> &residuals, DenseMatrix<double> &jacobian,
187  const unsigned& flag);
188 
189 };
190 
191 
192 //======================================================================
193 /// Refineable version of QPMLHelmholtzElement elements
194 ///
195 ///
196 //======================================================================
197 template <unsigned DIM, unsigned NNODE_1D>
199  public QPMLHelmholtzElement<DIM,NNODE_1D>,
200  public virtual RefineablePMLHelmholtzEquations<DIM>,
201  public virtual RefineableQElement<DIM>
202 {
203  public:
204 
205  /// \short Constructor, simply call the other constructors
209  RefineableQElement<DIM>(),
210  QPMLHelmholtzElement<DIM,NNODE_1D>()
211  {}
212 
213 
214  /// Broken copy constructor
216  dummy)
217  {
218  BrokenCopy::broken_copy("RefineableQuadPMLHelmholtzElement");
219  }
220 
221  /// Broken assignment operator
222  /*void operator=(const RefineableQPMLHelmholtzElement<DIM,NNODE_1D>&)
223  {
224  BrokenCopy::broken_assign("RefineableQuadPMLHelmholtzElement");
225  }*/
226 
227  /// Number of continuously interpolated values: 2
228  unsigned ncont_interpolated_values() const {return 2;}
229 
230  /// \short Number of vertex nodes in the element
231  unsigned nvertex_node() const
233 
234  /// \short Pointer to the j-th vertex node in the element
235  Node* vertex_node_pt(const unsigned& j) const
237 
238  /// Rebuild from sons: empty
239  void rebuild_from_sons(Mesh* &mesh_pt) {}
240 
241  /// \short Order of recovery shape functions for Z2 error estimation:
242  /// Same order as shape functions.
243  unsigned nrecovery_order() {return (NNODE_1D-1);}
244 
245  /// \short Perform additional hanging node procedures for variables
246  /// that are not interpolated by all nodes. Empty.
248 
249 };
250 
251 
252 
253 ////////////////////////////////////////////////////////////////////////
254 ////////////////////////////////////////////////////////////////////////
255 ////////////////////////////////////////////////////////////////////////
256 
257 
258 
259 //=======================================================================
260 /// Face geometry for the RefineableQuadPMLHelmholtzElement elements:
261 /// The spatial
262 /// dimension of the face elements is one lower than that of the
263 /// bulk element but they have the same number of points
264 /// along their 1D edges.
265 //=======================================================================
266 template<unsigned DIM, unsigned NNODE_1D>
268  public virtual QElement<DIM-1,NNODE_1D>
269 {
270 
271  public:
272 
273  /// \short Constructor: Call the constructor for the
274  /// appropriate lower-dimensional QElement
275  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
276 
277 };
278 
279 //=======================================================================
280 /// Policy class defining the elements to be used in the actual
281 /// PML layers. Same!
282 //=======================================================================
283 template<unsigned NNODE_1D>
285  public virtual RefineableQPMLHelmholtzElement<2,NNODE_1D>
286  {
287 
288  public:
289 
290  /// \short Constructor: Call the constructor for the
291  /// appropriate QElement
293  RefineableQPMLHelmholtzElement<2,NNODE_1D>()
294  {}
295 
296  };
297 
298 }
299 
300 #endif
301 
PMLHelmholtzSourceFctPt Source_fct_pt
Pointer to source function:
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
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.
void further_build()
Further build: Copy source function pointer from father element.
cstr elem_len * i
Definition: cfortran.h:607
PMLHelmholtzSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
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
char t
Definition: cfortran.h:572
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
RefineableQPMLHelmholtzElement()
Constructor, simply call the other constructors.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
RefineablePMLHelmholtzEquations(const RefineablePMLHelmholtzEquations< DIM > &dummy)
Broken copy constructor.
unsigned nvertex_node() const
Number of vertex nodes in the element.
static char t char * s
Definition: cfortran.h:572
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
RefineablePMLHelmholtzEquations()
Constructor, simply call other constructors.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Complex flux from PMLHelmholtz equations, strung together...
RefineableQPMLHelmholtzElement(const RefineableQPMLHelmholtzElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements...
Definition: elements.h:2371
unsigned num_Z2_flux_terms()
Broken assignment operator.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
PMLLayerElement()
Constructor: Call the constructor for the appropriate QElement.
unsigned ncont_interpolated_values() const
Broken assignment operator.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual unsigned nvertex_node() const
Definition: elements.h:2362
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
QPMLHelmholtzElement elements are linear/quadrilateral/ brick-shaped PMLHelmholtz elements with isopa...
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...
A general mesh class.
Definition: mesh.h:74
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
void fill_in_generic_residual_contribution_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...