Thelmholtz_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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
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_THELMHOLTZ_ELEMENTS_HEADER
32 #define OOMPH_THELMHOLTZ_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 "helmholtz_elements.h"
47 
48 namespace oomph
49 {
50 
51 /////////////////////////////////////////////////////////////////////////
52 /////////////////////////////////////////////////////////////////////////
53 // THelmholtzElement
54 ////////////////////////////////////////////////////////////////////////
55 ////////////////////////////////////////////////////////////////////////
56 
57 
58 
59 //======================================================================
60 /// THelmholtzElement<DIM,NNODE_1D> elements are isoparametric triangular
61 /// DIM-dimensional Helmholtz elements with NNODE_1D nodal points along each
62 /// element edge. Inherits from TElement and HelmholtzEquations
63 //======================================================================
64 template <unsigned DIM, unsigned NNODE_1D>
65 class THelmholtzElement : public TElement<DIM,NNODE_1D>,
66  public HelmholtzEquations<DIM>,
67  public virtual ElementWithZ2ErrorEstimator
68 {
69 
70  public:
71 
72  ///\short Constructor: Call constructors for TElement and
73  /// Helmholtz equations
74  THelmholtzElement() : TElement<DIM,NNODE_1D>(), HelmholtzEquations<DIM>()
75  { }
76 
77 
78  /// Broken copy constructor
80  {
81  BrokenCopy::broken_copy("THelmholtzElement");
82  }
83 
84  /// Broken assignment operator
85 //Commented out broken assignment operator because this can lead to a conflict warning
86 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
87 //realise that two separate implementations of the broken function are the same and so,
88 //quite rightly, it shouts.
89  /*void operator=(const THelmholtzElement<DIM,NNODE_1D>&)
90  {
91  BrokenCopy::broken_assign("THelmholtzElement");
92  }*/
93 
94  /// \short Access function for Nvalue: # of `values' (pinned or dofs)
95  /// at node n (always returns the same value at every node, 1)
96  inline unsigned required_nvalue(const unsigned &n) const
97  {return Initial_Nvalue;}
98 
99  /// \short Output function:
100  /// x,y,u or x,y,z,u
101  void output(std::ostream &outfile)
102  {
104  }
105 
106  /// \short Output function:
107  /// x,y,u or x,y,z,u at n_plot^DIM plot points
108  void output(std::ostream &outfile, const unsigned &n_plot)
109  {
110  HelmholtzEquations<DIM>::output(outfile,n_plot);
111  }
112 
113 
114  /// \short C-style output function:
115  /// x,y,u or x,y,z,u
116  void output(FILE* file_pt)
117  {
119  }
120 
121 
122  /// \short C-style output function:
123  /// x,y,u or x,y,z,u at n_plot^DIM plot points
124  void output(FILE* file_pt, const unsigned &n_plot)
125  {
126  HelmholtzEquations<DIM>::output(file_pt,n_plot);
127  }
128 
129 
130  /// \short Output function for an exact solution:
131  /// x,y,u_exact
132  void output_fct(std::ostream &outfile, const unsigned &n_plot,
134  {
135  HelmholtzEquations<DIM>::output_fct(outfile,n_plot,exact_soln_pt);
136  }
137 
138 
139  /// \short Output function for a time-dependent exact solution.
140  /// x,y,u_exact (calls the steady version)
141  void output_fct(std::ostream &outfile, const unsigned &n_plot,
142  const double& time,
144  {
145  HelmholtzEquations<DIM>::output_fct(outfile,n_plot,time,exact_soln_pt);
146  }
147 
148 protected:
149 
150 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
152  Shape &psi,
153  DShape &dpsidx,
154  Shape &test,
155  DShape &dtestdx) const;
156 
157 
158 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
159  inline double dshape_and_dtest_eulerian_at_knot_helmholtz(const unsigned &ipt,
160  Shape &psi,
161  DShape &dpsidx,
162  Shape &test,
163  DShape &dtestdx)
164  const;
165 
166 
167  /// \short Order of recovery shape functions for Z2 error estimation:
168  /// Same order as shape functions.
169  unsigned nrecovery_order() {return (NNODE_1D-1);}
170 
171  /// Number of 'flux' terms for Z2 error estimation
172  unsigned num_Z2_flux_terms() {return 2*DIM;}
173 
174  /// \short Get 'flux' for Z2 error recovery: Standard flux from
175  /// UnsteadyHeat equations
177  {
178  Vector<std::complex <double> > complex_flux(DIM);
179  this->get_flux(s,complex_flux);
180  unsigned count=0;
181  for (unsigned i=0;i<DIM;i++)
182  {
183  flux[count++]=complex_flux[i].real();
184  flux[count++]=complex_flux[i].imag();
185  }
186  }
187 
188  /// \short Number of vertex nodes in the element
189  unsigned nvertex_node() const
191 
192  /// \short Pointer to the j-th vertex node in the element
193  Node* vertex_node_pt(const unsigned& j) const
195 
196  private:
197 
198  /// Static unsigned that holds the (same) number of variables at every node
199  static const unsigned Initial_Nvalue;
200 
201 
202 };
203 
204 
205 
206 
207 //Inline functions:
208 
209 
210 //======================================================================
211 /// Define the shape functions and test functions and derivatives
212 /// w.r.t. global coordinates and return Jacobian of mapping.
213 ///
214 /// Galerkin: Test functions = shape functions
215 //======================================================================
216 template<unsigned DIM, unsigned NNODE_1D>
218  const Vector<double> &s,
219  Shape &psi,
220  DShape &dpsidx,
221  Shape &test,
222  DShape &dtestdx) const
223 {
224  unsigned n_node = this->nnode();
225 
226  //Call the geometrical shape functions and derivatives
227  double J = this->dshape_eulerian(s,psi,dpsidx);
228 
229  //Loop over the test functions and derivatives and set them equal to the
230  //shape functions
231  for(unsigned i=0;i<n_node;i++)
232  {
233  test[i] = psi[i];
234  dtestdx(i,0) = dpsidx(i,0);
235  dtestdx(i,1) = dpsidx(i,1);
236  }
237 
238  //Return the jacobian
239  return J;
240 }
241 
242 
243 
244 //======================================================================
245 /// Define the shape functions and test functions and derivatives
246 /// w.r.t. global coordinates and return Jacobian of mapping.
247 ///
248 /// Galerkin: Test functions = shape functions
249 //======================================================================
250 template<unsigned DIM, unsigned NNODE_1D>
253  const unsigned &ipt,
254  Shape &psi,
255  DShape &dpsidx,
256  Shape &test,
257  DShape &dtestdx) const
258 {
259 
260  //Call the geometrical shape functions and derivatives
261  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
262 
263  //Set the pointers of the test functions
264  test = psi;
265  dtestdx = dpsidx;
266 
267  //Return the jacobian
268  return J;
269 
270 }
271 
272 
273 
274 //=======================================================================
275 /// Face geometry for the THelmholtzElement elements: The spatial
276 /// dimension of the face elements is one lower than that of the
277 /// bulk element but they have the same number of points
278 /// along their 1D edges.
279 //=======================================================================
280 template<unsigned DIM, unsigned NNODE_1D>
281 class FaceGeometry<THelmholtzElement<DIM,NNODE_1D> >:
282  public virtual TElement<DIM-1,NNODE_1D>
283 {
284 
285  public:
286 
287  /// \short Constructor: Call the constructor for the
288  /// appropriate lower-dimensional TElement
289  FaceGeometry() : TElement<DIM-1,NNODE_1D>() {}
290 
291 };
292 
293 //=======================================================================
294 /// Face geometry for the 1D THelmholtzElement elements: Point elements
295 //=======================================================================
296 template<unsigned NNODE_1D>
297 class FaceGeometry<THelmholtzElement<1,NNODE_1D> >:
298  public virtual PointElement
299 {
300 
301  public:
302 
303  /// \short Constructor: Call the constructor for the
304  /// appropriate lower-dimensional TElement
306 
307 };
308 
309 
310 }
311 
312 #endif
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement. ...
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 output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,u_exact.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
unsigned nvertex_node() const
Number of vertex nodes in the element.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional TElement. ...
cstr elem_len * i
Definition: cfortran.h:607
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.
THelmholtzElement()
Constructor: Call constructors for TElement and Helmholtz equations.
static const unsigned Initial_Nvalue
Static unsigned that holds the (same) number of variables at every node.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
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 get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Standard flux from UnsteadyHeat equations.
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.
static char t char * s
Definition: cfortran.h:572
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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(* 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 (calls the steady version) ...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
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.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
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
THelmholtzElement(const THelmholtzElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
void output(std::ostream &outfile)
Output with default number of plot points.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.