helmholtz_flux_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 elements that are used to apply prescribed flux
31 // boundary conditions to the Helmholtz equations
32 #ifndef OOMPH_HELMHOLTZ_FLUX_ELEMENTS_HEADER
33 #define OOMPH_HELMHOLTZ_FLUX_ELEMENTS_HEADER
34 
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 // oomph-lib ncludes
42 #include "../generic/Qelements.h"
43 #include "math.h"
44 #include <complex>
45 
46 namespace oomph
47 {
48 
49 //======================================================================
50 /// \short A class for elements that allow the imposition of an
51 /// applied flux on the boundaries of Helmholtz elements.
52 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
53 /// policy class.
54 //======================================================================
55 template <class ELEMENT>
56 class HelmholtzFluxElement : public virtual FaceGeometry<ELEMENT>,
57 public virtual FaceElement
58 {
59 
60 public:
61 
62  /// \short Function pointer to the prescribed-flux function fct(x,f(x)) --
63  /// x is a Vector and the flux is a complex
64 
65  typedef void (*HelmholtzPrescribedFluxFctPt)
66  (const Vector<double>& x, std::complex<double>& flux);
67 
68  /// \short Constructor, takes the pointer to the "bulk" element and the
69  /// index of the face to which the element is attached.
70  HelmholtzFluxElement(FiniteElement* const &bulk_el_pt,
71  const int& face_index);
72 
73  ///\short Broken empty constructor
75  {
76  throw OomphLibError(
77  "Don't call empty constructor for HelmholtzFluxElement",
78  OOMPH_CURRENT_FUNCTION,
79  OOMPH_EXCEPTION_LOCATION);
80  }
81 
82  /// Broken copy constructor
84  {
85  BrokenCopy::broken_copy("HelmholtzFluxElement");
86  }
87 
88  /// Broken assignment operator
89 //Commented out broken assignment operator because this can lead to a conflict warning
90 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
91 //realise that two separate implementations of the broken function are the same and so,
92 //quite rightly, it shouts.
93  /*void operator=(const HelmholtzFluxElement&)
94  {
95  BrokenCopy::broken_assign("HelmholtzFluxElement");
96  }*/
97 
98 
99  /// Access function for the prescribed-flux function pointer
101 
102 
103  /// Add the element's contribution to its residual vector
105  {
106  //Call the generic residuals function with flag set to 0
107  //using a dummy matrix argument
110  }
111 
112 
113  /// \short Add the element's contribution to its residual vector and its
114  /// Jacobian matrix
116  DenseMatrix<double> &jacobian)
117  {
118  //Call the generic routine with the flag set to 1
120  }
121 
122 
123  /// \short Specify the value of nodal zeta from the face geometry
124  /// The "global" intrinsic coordinate of the element when
125  /// viewed as part of a geometric object should be given by
126  /// the FaceElement representation, by default (needed to break
127  /// indeterminacy if bulk element is SolidElement)
128  double zeta_nodal(const unsigned &n, const unsigned &k,
129  const unsigned &i) const
130  {return FaceElement::zeta_nodal(n,k,i);}
131 
132 
133  /// Output function -- forward to broken version in FiniteElement
134  /// until somebody decides what exactly they want to plot here...
135  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
136 
137  /// \short Output function -- forward to broken version in FiniteElement
138  /// until somebody decides what exactly they want to plot here...
139  void output(std::ostream &outfile, const unsigned &n_plot)
140  {FiniteElement::output(outfile,n_plot);}
141 
142 
143  /// C-style output function -- forward to broken version in FiniteElement
144  /// until somebody decides what exactly they want to plot here...
145  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
146 
147  /// \short C-style output function -- forward to broken version in
148  /// FiniteElement until somebody decides what exactly they want to plot
149  /// here...
150  void output(FILE* file_pt, const unsigned &n_plot)
151  {FiniteElement::output(file_pt,n_plot);}
152 
153 
154  /// \short Return the index at which the unknown value
155  /// is stored.
156  virtual inline std::complex<unsigned> u_index_helmholtz() const
157  {return std::complex<unsigned>(U_index_helmholtz.real(),U_index_helmholtz.imag());}
158 
159 
160 protected:
161 
162  /// \short Function to compute the shape and test functions and to return
163  /// the Jacobian of mapping between local and global (Eulerian)
164  /// coordinates
165  inline double shape_and_test(const Vector<double> &s, Shape &psi, Shape &test)
166  const
167  {
168  //Find number of nodes
169  unsigned n_node = nnode();
170 
171  //Get the shape functions
172  shape(s,psi);
173 
174  //Set the test functions to be the same as the shape functions
175  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
176 
177  //Return the value of the jacobian
178  return J_eulerian(s);
179  }
180 
181 
182  /// \short Function to compute the shape and test functions and to return
183  /// the Jacobian of mapping between local and global (Eulerian)
184  /// coordinates
185  inline double shape_and_test_at_knot(const unsigned &ipt,
186  Shape &psi, Shape &test)
187  const
188  {
189  //Find number of nodes
190  unsigned n_node = nnode();
191 
192  //Get the shape functions
193  shape_at_knot(ipt,psi);
194 
195  //Set the test functions to be the same as the shape functions
196  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
197 
198  //Return the value of the jacobian
199  return J_eulerian_at_knot(ipt);
200  }
201 
202 
203  /// Function to calculate the prescribed flux at a given spatial
204  /// position
205  void get_flux(const Vector<double>& x, std::complex<double>& flux)
206  {
207  //If the function pointer is zero return zero
208  if(Flux_fct_pt == 0)
209  {
210  flux = std::complex<double>(0.0,0.0);
211  }
212  //Otherwise call the function
213  else
214  {
215  (*Flux_fct_pt)(x,flux);
216  }
217  }
218 
219 
220  /// \short The index at which the real and imag part of the unknown is stored
221  /// at the nodes
222  std::complex<unsigned> U_index_helmholtz;
223 
224 
225  /// \short Add the element's contribution to its residual vector.
226  /// flag=1(or 0): do (or don't) compute the contribution to the
227  /// Jacobian as well.
229  Vector<double> &residuals, DenseMatrix<double> &jacobian,
230  const unsigned& flag);
231 
232 
233  /// Function pointer to the (global) prescribed-flux function
235 
236  ///The spatial dimension of the problem
237  unsigned Dim;
238 
239 
240 };
241 
242 //////////////////////////////////////////////////////////////////////
243 //////////////////////////////////////////////////////////////////////
244 //////////////////////////////////////////////////////////////////////
245 
246 
247 
248 //===========================================================================
249 /// Constructor, takes the pointer to the "bulk" element, the
250 /// index of the fixed local coordinate and its value represented
251 /// by an integer (+/- 1), indicating that the face is located
252 /// at the max. or min. value of the "fixed" local coordinate
253 /// in the bulk element.
254 //===========================================================================
255 template<class ELEMENT>
258  const int &face_index) :
259  FaceGeometry<ELEMENT>(), FaceElement()
260  {
261 
262  // Let the bulk element build the FaceElement, i.e. setup the pointers
263  // to its nodes (by referring to the appropriate nodes in the bulk
264  // element), etc.
265  bulk_el_pt->build_face_element(face_index,this);
266 
267 #ifdef PARANOID
268  {
269  //Check that the element is not a refineable 3d element
270  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
271  //If it's three-d
272  if(elem_pt->dim()==3)
273  {
274  //Is it refineable
275  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
276  if(ref_el_pt!=0)
277  {
278  if (this->has_hanging_nodes())
279  {
280  throw OomphLibError(
281  "This flux element will not work correctly if nodes are hanging\n",
282  OOMPH_CURRENT_FUNCTION,
283  OOMPH_EXCEPTION_LOCATION);
284  }
285  }
286  }
287  }
288 #endif
289 
290  // Initialise the prescribed-flux function pointer to zero
291  Flux_fct_pt = 0;
292 
293  // Extract the dimension of the problem from the dimension of
294  // the first node
295  Dim = this->node_pt(0)->ndim();
296 
297  //Set up U_index_helmholtz. Initialise to zero, which probably won't change
298  //in most cases, oh well, the price we pay for generality
299  U_index_helmholtz = std::complex<unsigned>(0,1);
300 
301  //Cast to the appropriate HelmholtzEquation so that we can
302  //find the index at which the variable is stored
303  //We assume that the dimension of the full problem is the same
304  //as the dimension of the node, if this is not the case you will have
305  //to write custom elements, sorry
306  switch(Dim)
307  {
308  //One dimensional problem
309  case 1:
310  {
311  HelmholtzEquations<1>* eqn_pt =
312  dynamic_cast<HelmholtzEquations<1>*>(bulk_el_pt);
313  //If the cast has failed die
314  if(eqn_pt==0)
315  {
316  std::string error_string =
317  "Bulk element must inherit from HelmholtzEquations.";
318  error_string +=
319  "Nodes are one dimensional, but cannot cast the bulk element to\n";
320  error_string += "HelmholtzEquations<1>\n.";
321  error_string +=
322  "If you desire this functionality, you must implement it yourself\n";
323 
324  throw OomphLibError(error_string,
325  OOMPH_CURRENT_FUNCTION,
326  OOMPH_EXCEPTION_LOCATION);
327  }
328  //Otherwise read out the value
329  else
330  {
331  //Read the index from the (cast) bulk element
333  }
334  }
335  break;
336 
337  //Two dimensional problem
338  case 2:
339  {
340  HelmholtzEquations<2>* eqn_pt =
341  dynamic_cast<HelmholtzEquations<2>*>(bulk_el_pt);
342  //If the cast has failed die
343  if(eqn_pt==0)
344  {
345  std::string error_string =
346  "Bulk element must inherit from HelmholtzEquations.";
347  error_string +=
348  "Nodes are two dimensional, but cannot cast the bulk element to\n";
349  error_string += "HelmholtzEquations<2>\n.";
350  error_string +=
351  "If you desire this functionality, you must implement it yourself\n";
352 
353  throw OomphLibError(error_string,
354  OOMPH_CURRENT_FUNCTION,
355  OOMPH_EXCEPTION_LOCATION);
356  }
357  else
358  {
359  //Read the index from the (cast) bulk element
361  }
362  }
363 
364  break;
365 
366  //Three dimensional problem
367  case 3:
368  {
369  HelmholtzEquations<3>* eqn_pt =
370  dynamic_cast<HelmholtzEquations<3>*>(bulk_el_pt);
371  //If the cast has failed die
372  if(eqn_pt==0)
373  {
374  std::string error_string =
375  "Bulk element must inherit from HelmholtzEquations.";
376  error_string +=
377  "Nodes are three dimensional, but cannot cast the bulk element to\n";
378  error_string += "HelmholtzEquations<3>\n.";
379  error_string +=
380  "If you desire this functionality, you must implement it yourself\n";
381 
382  throw OomphLibError(error_string,
383  OOMPH_CURRENT_FUNCTION,
384  OOMPH_EXCEPTION_LOCATION);
385  }
386  else
387  {
388  //Read the index from the (cast) bulk element
390  }
391  }
392  break;
393 
394  //Any other case is an error
395  default:
396  std::ostringstream error_stream;
397  error_stream << "Dimension of node is " << Dim
398  << ". It should be 1,2, or 3!" << std::endl;
399 
400  throw OomphLibError(error_stream.str(),
401  OOMPH_CURRENT_FUNCTION,
402  OOMPH_EXCEPTION_LOCATION);
403  break;
404  }
405  }
406 
407 
408 //===========================================================================
409 /// Compute the element's residual vector and the (zero) Jacobian matrix.
410 //===========================================================================
411 template<class ELEMENT>
414  Vector<double> &residuals, DenseMatrix<double> &jacobian,
415  const unsigned& flag)
416 {
417  //Find out how many nodes there are
418  const unsigned n_node = nnode();
419 
420  //Set up memory for the shape and test functions
421  Shape psif(n_node), testf(n_node);
422 
423  //Set the value of Nintpt
424  const unsigned n_intpt = integral_pt()->nweight();
425 
426  //Set the Vector to hold local coordinates
427  Vector<double> s(Dim-1);
428 
429  //Integers to hold the local equation and unknown numbers
430  int local_eqn_real=0 ,local_eqn_imag=0;
431 
432  //Loop over the integration points
433  //--------------------------------
434  for(unsigned ipt=0;ipt<n_intpt;ipt++)
435  {
436 
437  //Assign values of s
438  for(unsigned i=0;i<(Dim-1);i++) {s[i] = integral_pt()->knot(ipt,i);}
439 
440  //Get the integral weight
441  double w = integral_pt()->weight(ipt);
442 
443  //Find the shape and test functions and return the Jacobian
444  //of the mapping
445  double J = shape_and_test(s,psif,testf);
446 
447  //Premultiply the weights and the Jacobian
448  double W = w*J;
449 
450  //Need to find position to feed into flux function, initialise to zero
451  Vector<double> interpolated_x(Dim,0.0);
452 
453  //Calculate Eulerian position of integration point
454  for(unsigned l=0;l<n_node;l++)
455  {
456  for(unsigned i=0;i<Dim;i++)
457  {
458  interpolated_x[i] += nodal_position(l,i)*psif[l];
459  }
460  }
461 
462  //Get the imposed flux
463  std::complex<double> flux(0.0,0.0);
464  get_flux(interpolated_x,flux);
465 
466  //Now add to the appropriate equations
467  //Loop over the test functions
468  for(unsigned l=0;l<n_node;l++)
469  {
470  local_eqn_real = nodal_local_eqn(l,U_index_helmholtz.real());
471  /*IF it's not a boundary condition*/
472  if(local_eqn_real >= 0)
473  {
474  //Add the prescribed flux terms
475  residuals[local_eqn_real] -= flux.real()*testf[l]*W;
476 
477  // Imposed traction doesn't depend upon the solution,
478  // --> the Jacobian is always zero, so no Jacobian
479  // terms are required
480  }
481  local_eqn_imag = nodal_local_eqn(l,U_index_helmholtz.imag());
482  /*IF it's not a boundary condition*/
483  if(local_eqn_imag >= 0)
484  {
485  //Add the prescribed flux terms
486  residuals[local_eqn_imag] -= flux.imag()*testf[l]*W;
487 
488  // Imposed traction doesn't depend upon the solution,
489  // --> the Jacobian is always zero, so no Jacobian
490  // terms are required
491  }
492  }
493  }
494 }
495 
496 
497 }
498 
499 #endif
500 
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
HelmholtzPrescribedFluxFctPt & flux_fct_pt()
Broken assignment operator.
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5121
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
HelmholtzFluxElement(const HelmholtzFluxElement &dummy)
Broken copy constructor.
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
std::complex< unsigned > U_index_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
A general Finite Element class.
Definition: elements.h:1271
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
HelmholtzFluxElement()
Broken empty constructor.
virtual std::complex< unsigned > u_index_helmholtz() const
Return the index at which the unknown value is stored.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
void get_flux(const Vector< double > &x, std::complex< double > &flux)
void(* HelmholtzPrescribedFluxFctPt)(const Vector< double > &x, std::complex< double > &flux)
Function pointer to the prescribed-flux function fct(x,f(x)) – x is a Vector and the flux is a comple...
virtual void fill_in_generic_residual_contribution_helmholtz_flux(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the contri...
void output(std::ostream &outfile)
static char t char * s
Definition: cfortran.h:572
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3158
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:5033
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
HelmholtzPrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned Dim
The spatial dimension of the problem.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
A class for elements that allow the imposition of an applied flux on the boundaries of Helmholtz elem...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
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