Tpml_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 THelmholtz elements
31 #ifndef OOMPH_TPML_HELMHOLTZ_ELEMENTS_HEADER
32 #define OOMPH_TPML_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 
41 //OOMPH-LIB headers
42 #include "../generic/nodes.h"
43 #include "../generic/oomph_utilities.h"
44 #include "../generic/Telements.h"
45 #include "../generic/error_estimator.h"
46 #include "pml_helmholtz_elements.h"
47 
48 namespace oomph
49 {
50 
51 /////////////////////////////////////////////////////////////////////////
52 /////////////////////////////////////////////////////////////////////////
53 // TPMLHelmholtzElement
54 ////////////////////////////////////////////////////////////////////////
55 ////////////////////////////////////////////////////////////////////////
56 
57 
58 
59 //======================================================================
60 /// TPMLHelmholtzElement<DIM,NNODE_1D> elements are
61 /// isoparametric triangular DIM-dimensional PMLHelmholtz elements
62 /// with NNODE_1D nodal points along each element edge. Inherits from
63 /// TElement and PMLHelmholtzEquations
64 //======================================================================
65 template <unsigned DIM, unsigned NNODE_1D>
66 class TPMLHelmholtzElement : public TElement<DIM,NNODE_1D>,
67  public PMLHelmholtzEquations<DIM>,
68  public virtual ElementWithZ2ErrorEstimator
69 {
70 
71  public:
72 
73  ///\short Constructor: Call constructors for TElement and
74  /// PMLHelmholtz equations
76  TElement<DIM,NNODE_1D>(), PMLHelmholtzEquations<DIM>()
77  { }
78 
79 
80  /// Broken copy constructor
83  {
84  BrokenCopy::broken_copy("TPMLHelmholtzElement");
85  }
86 
87  /// Broken assignment operator
88 //Commented out broken assignment operator because this can lead to a conflict warning
89 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
90 //realise that two separate implementations of the broken function are the same and so,
91 //quite rightly, it shouts.
92  /*void operator=(const TPMLHelmholtzElement<DIM,NNODE_1D>&)
93  {
94  BrokenCopy::broken_assign("TPMLHelmholtzElement");
95  }*/
96 
97  /// \short Access function for Nvalue: # of `values' (pinned or dofs)
98  /// at node n (always returns the same value at every node, 1)
99  inline unsigned required_nvalue(const unsigned &n) const
100  {return Initial_Nvalue;}
101 
102  /// \short Output function:
103  /// x,y,u or x,y,z,u
104  void output(std::ostream &outfile)
105  {
107  }
108 
109  /// \short Output function:
110  /// x,y,u or x,y,z,u at n_plot^DIM plot points
111  void output(std::ostream &outfile, const unsigned &n_plot)
112  {
113  PMLHelmholtzEquations<DIM>::output(outfile,n_plot);
114  }
115 
116 
117  /// \short C-style output function:
118  /// x,y,u or x,y,z,u
119  void output(FILE* file_pt)
120  {
122  }
123 
124 
125  /// \short C-style output function:
126  /// x,y,u or x,y,z,u at n_plot^DIM plot points
127  void output(FILE* file_pt, const unsigned &n_plot)
128  {
129  PMLHelmholtzEquations<DIM>::output(file_pt,n_plot);
130  }
131 
132 
133  /// \short Output function for an exact solution:
134  /// x,y,u_exact
135  void output_fct(std::ostream &outfile, const unsigned &n_plot,
137  {
138  PMLHelmholtzEquations<DIM>::output_fct(outfile,n_plot,exact_soln_pt);
139  }
140 
141 
142  /// \short Output function for a time-dependent exact solution.
143  /// x,y,u_exact (calls the steady version)
144  void output_fct(std::ostream &outfile, const unsigned &n_plot,
145  const double& time,
147  {
149  n_plot,
150  time,
151  exact_soln_pt);
152  }
153 
154 protected:
155 
156 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
158  Shape &psi,
159  DShape &dpsidx,
160  Shape &test,
161  DShape &dtestdx) const;
162 
163 
164 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
165  inline double dshape_and_dtest_eulerian_at_knot_helmholtz(const unsigned &ipt,
166  Shape &psi,
167  DShape &dpsidx,
168  Shape &test,
169  DShape &dtestdx)
170  const;
171 
172 
173  /// \short Order of recovery shape functions for Z2 error estimation:
174  /// Same order as shape functions.
175  unsigned nrecovery_order() {return (NNODE_1D-1);}
176 
177  /// Number of 'flux' terms for Z2 error estimation
178  unsigned num_Z2_flux_terms() {return 2*DIM;}
179 
180  /// \short Get 'flux' for Z2 error recovery: Standard flux from
181  /// UnsteadyHeat equations
183  {
184  Vector<std::complex <double> > complex_flux(DIM);
185  this->get_flux(s,complex_flux);
186  unsigned count=0;
187  for (unsigned i=0;i<DIM;i++)
188  {
189  flux[count++]=complex_flux[i].real();
190  flux[count++]=complex_flux[i].imag();
191  }
192  }
193 
194  /// \short Number of vertex nodes in the element
195  unsigned nvertex_node() const
197 
198  /// \short Pointer to the j-th vertex node in the element
199  Node* vertex_node_pt(const unsigned& j) const
201 
202  private:
203 
204  /// Static unsigned that holds the (same) number of variables at every node
205  static const unsigned Initial_Nvalue;
206 
207 
208 };
209 
210 
211 //!! Cleanup - this was not here before!
212 template<unsigned DIM, unsigned NNODE_1D>
214 
215 //Inline functions:
216 
217 
218 //======================================================================
219 /// Define the shape functions and test functions and derivatives
220 /// w.r.t. global coordinates and return Jacobian of mapping.
221 ///
222 /// Galerkin: Test functions = shape functions
223 //======================================================================
224 template<unsigned DIM, unsigned NNODE_1D>
225 double
227  (
228  const Vector<double> &s,
229  Shape &psi,
230  DShape &dpsidx,
231  Shape &test,
232  DShape &dtestdx) const
233 {
234  unsigned n_node = this->nnode();
235 
236  //Call the geometrical shape functions and derivatives
237  double J = this->dshape_eulerian(s,psi,dpsidx);
238 
239  //Loop over the test functions and derivatives and set them equal to the
240  //shape functions
241  for(unsigned i=0;i<n_node;i++)
242  {
243  test[i] = psi[i];
244  dtestdx(i,0) = dpsidx(i,0);
245  dtestdx(i,1) = dpsidx(i,1);
246  }
247 
248  //Return the jacobian
249  return J;
250 }
251 
252 
253 
254 //======================================================================
255 /// Define the shape functions and test functions and derivatives
256 /// w.r.t. global coordinates and return Jacobian of mapping.
257 ///
258 /// Galerkin: Test functions = shape functions
259 //======================================================================
260 template<unsigned DIM, unsigned NNODE_1D>
263  const unsigned &ipt,
264  Shape &psi,
265  DShape &dpsidx,
266  Shape &test,
267  DShape &dtestdx) const
268 {
269 
270  //Call the geometrical shape functions and derivatives
271  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
272 
273  //Set the pointers of the test functions
274  test = psi;
275  dtestdx = dpsidx;
276 
277  //Return the jacobian
278  return J;
279 
280 }
281 
282 
283 
284 //=======================================================================
285 /// Face geometry for the TPMLHelmholtzElement elements: The spatial
286 /// dimension of the face elements is one lower than that of the
287 /// bulk element but they have the same number of points
288 /// along their 1D edges.
289 //=======================================================================
290 template<unsigned DIM, unsigned NNODE_1D>
291 class FaceGeometry<TPMLHelmholtzElement<DIM,NNODE_1D> >:
292  public virtual TElement<DIM-1,NNODE_1D>
293 {
294 
295  public:
296 
297  /// \short Constructor: Call the constructor for the
298  /// appropriate lower-dimensional TElement
299  FaceGeometry() : TElement<DIM-1,NNODE_1D>() {}
300 
301 };
302 
303 //=======================================================================
304 /// Face geometry for the 1D TPMLHelmholtzElement elements:
305 /// Point elements
306 //=======================================================================
307 template<unsigned NNODE_1D>
308 class FaceGeometry<TPMLHelmholtzElement<1,NNODE_1D> >:
309  public virtual PointElement
310 {
311 
312  public:
313 
314  /// \short Constructor: Call the constructor for the
315  /// appropriate lower-dimensional TElement
317 
318 };
319 
320 
321 ////////////////////////////////////////////////////////////////////////
322 ////////////////////////////////////////////////////////////////////////
323 ////////////////////////////////////////////////////////////////////////
324 
325 //=======================================================================
326 /// Policy class defining the elements to be used in the actual
327 /// PML layers. It's the corresponding quads.
328 //=======================================================================
329  template<unsigned NNODE_1D>
330 class PMLLayerElement<TPMLHelmholtzElement<2,NNODE_1D> > :
331  public virtual QPMLHelmholtzElement<2,NNODE_1D>
332 {
333 
334  public:
335 
336  /// \short Constructor: Call the constructor for the
337  /// appropriate QElement
339  {}
340 
341 };
342 
343 //=======================================================================
344 /// Face geometry for the TPMLHelmholtzElement elements: The spatial
345 /// dimension of the face elements is one lower than that of the
346 /// bulk element but they have the same number of points
347 /// along their 1D edges.
348 //=======================================================================
349 template<unsigned DIM, unsigned NNODE_1D>
351  public virtual QElement<DIM-1,NNODE_1D>
352 {
353 
354  public:
355 
356  /// \short Constructor: Call the constructor for the
357  /// appropriate lower-dimensional TElement
358  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
359 
360 };
361 
362 ////////////////////////////////////////////////////////////////////////
363 ////////////////////////////////////////////////////////////////////////
364 ////////////////////////////////////////////////////////////////////////
365 
366 //=======================================================================
367 /// Policy class defining the elements to be used in the actual
368 /// PML layers. It's the corresponding quads.
369 //=======================================================================
370  template<unsigned NNODE_1D>
372  <TPMLHelmholtzElement<2,NNODE_1D> > > :
373  public virtual QPMLHelmholtzElement<2,NNODE_1D>
374 {
375 
376  public:
377 
378  /// \short Constructor: Call the constructor for the
379  /// appropriate QElement
381  {}
382 
383 };
384 
385 }
386 
387 
388 
389 #endif
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
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 output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
static const unsigned Initial_Nvalue
Static unsigned that holds the (same) number of variables at every node.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement. ...
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement. ...
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.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement. ...
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 (calls the steady version) ...
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
static char t char * s
Definition: cfortran.h:572
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.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
TPMLHelmholtzElement()
Constructor: Call constructors for TElement and PMLHelmholtz equations.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
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.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
PMLLayerElement()
Constructor: Call the constructor for the appropriate QElement.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
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 get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Standard flux from UnsteadyHeat equations.
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.
QPMLHelmholtzElement elements are linear/quadrilateral/ brick-shaped PMLHelmholtz elements with isopa...
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. Return Jacobian.