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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for elements that are used to apply prescribed flux
27 // boundary conditions to the Helmholtz equations
28 #ifndef OOMPH_HELMHOLTZ_FLUX_ELEMENTS_HEADER
29 #define OOMPH_HELMHOLTZ_FLUX_ELEMENTS_HEADER
30 
31 
32 // Config header generated by autoconfig
33 #ifdef HAVE_CONFIG_H
34 #include <oomph-lib-config.h>
35 #endif
36 
37 // oomph-lib ncludes
38 #include "../generic/Qelements.h"
39 #include "math.h"
40 #include <complex>
41 
42 namespace oomph
43 {
44  //======================================================================
45  /// A class for elements that allow the imposition of an
46  /// applied flux on the boundaries of Helmholtz elements.
47  /// The element geometry is obtained from the FaceGeometry<ELEMENT>
48  /// policy class.
49  //======================================================================
50  template<class ELEMENT>
51  class HelmholtzFluxElement : public virtual FaceGeometry<ELEMENT>,
52  public virtual FaceElement
53  {
54  public:
55  /// Function pointer to the prescribed-flux function fct(x,f(x)) --
56  /// x is a Vector and the flux is a complex
57 
58  typedef void (*HelmholtzPrescribedFluxFctPt)(const Vector<double>& x,
59  std::complex<double>& flux);
60 
61  /// Constructor, takes the pointer to the "bulk" element and the
62  /// index of the face to which the element is attached.
63  HelmholtzFluxElement(FiniteElement* const& bulk_el_pt,
64  const int& face_index);
65 
66  /// Broken empty constructor
68  {
69  throw OomphLibError(
70  "Don't call empty constructor for HelmholtzFluxElement",
71  OOMPH_CURRENT_FUNCTION,
72  OOMPH_EXCEPTION_LOCATION);
73  }
74 
75  /// Broken copy constructor
77 
78  /// Broken assignment operator
79  // Commented out broken assignment operator because this can lead to a
80  // conflict warning when used in the virtual inheritence hierarchy.
81  // Essentially the compiler doesn't realise that two separate
82  // implementations of the broken function are the same and so, quite
83  // rightly, it shouts.
84  /*void operator=(const HelmholtzFluxElement&) = delete;*/
85 
86 
87  /// Access function for the prescribed-flux function pointer
89  {
90  return Flux_fct_pt;
91  }
92 
93 
94  /// Add the element's contribution to its residual vector
96  {
97  // Call the generic residuals function with flag set to 0
98  // using a dummy matrix argument
100  residuals, GeneralisedElement::Dummy_matrix, 0);
101  }
102 
103 
104  /// Add the element's contribution to its residual vector and its
105  /// Jacobian matrix
107  DenseMatrix<double>& jacobian)
108  {
109  // Call the generic routine with the flag set to 1
111  residuals, jacobian, 1);
112  }
113 
114 
115  /// Specify the value of nodal zeta from the face geometry
116  /// The "global" intrinsic coordinate of the element when
117  /// viewed as part of a geometric object should be given by
118  /// the FaceElement representation, by default (needed to break
119  /// indeterminacy if bulk element is SolidElement)
120  double zeta_nodal(const unsigned& n,
121  const unsigned& k,
122  const unsigned& i) const
123  {
124  return FaceElement::zeta_nodal(n, k, i);
125  }
126 
127 
128  /// Output function -- forward to broken version in FiniteElement
129  /// until somebody decides what exactly they want to plot here...
130  void output(std::ostream& outfile)
131  {
132  FiniteElement::output(outfile);
133  }
134 
135  /// Output function -- forward to broken version in FiniteElement
136  /// until somebody decides what exactly they want to plot here...
137  void output(std::ostream& outfile, const unsigned& n_plot)
138  {
139  FiniteElement::output(outfile, n_plot);
140  }
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)
146  {
147  FiniteElement::output(file_pt);
148  }
149 
150  /// C-style output function -- forward to broken version in
151  /// FiniteElement until somebody decides what exactly they want to plot
152  /// here...
153  void output(FILE* file_pt, const unsigned& n_plot)
154  {
155  FiniteElement::output(file_pt, n_plot);
156  }
157 
158 
159  /// Return the index at which the unknown value
160  /// is stored.
161  virtual inline std::complex<unsigned> u_index_helmholtz() const
162  {
163  return std::complex<unsigned>(U_index_helmholtz.real(),
164  U_index_helmholtz.imag());
165  }
166 
167 
168  protected:
169  /// Function to compute the shape and test functions and to return
170  /// the Jacobian of mapping between local and global (Eulerian)
171  /// coordinates
172  inline double shape_and_test(const Vector<double>& s,
173  Shape& psi,
174  Shape& test) const
175  {
176  // Find number of nodes
177  unsigned n_node = nnode();
178 
179  // Get the shape functions
180  shape(s, psi);
181 
182  // Set the test functions to be the same as the shape functions
183  for (unsigned i = 0; i < n_node; i++)
184  {
185  test[i] = psi[i];
186  }
187 
188  // Return the value of the jacobian
189  return J_eulerian(s);
190  }
191 
192 
193  /// Function to compute the shape and test functions and to return
194  /// the Jacobian of mapping between local and global (Eulerian)
195  /// coordinates
196  inline double shape_and_test_at_knot(const unsigned& ipt,
197  Shape& psi,
198  Shape& test) const
199  {
200  // Find number of nodes
201  unsigned n_node = nnode();
202 
203  // Get the shape functions
204  shape_at_knot(ipt, psi);
205 
206  // Set the test functions to be the same as the shape functions
207  for (unsigned i = 0; i < n_node; i++)
208  {
209  test[i] = psi[i];
210  }
211 
212  // Return the value of the jacobian
213  return J_eulerian_at_knot(ipt);
214  }
215 
216 
217  /// Function to calculate the prescribed flux at a given spatial
218  /// position
219  void get_flux(const Vector<double>& x, std::complex<double>& flux)
220  {
221  // If the function pointer is zero return zero
222  if (Flux_fct_pt == 0)
223  {
224  flux = std::complex<double>(0.0, 0.0);
225  }
226  // Otherwise call the function
227  else
228  {
229  (*Flux_fct_pt)(x, flux);
230  }
231  }
232 
233 
234  /// The index at which the real and imag part of the unknown is
235  /// stored at the nodes
236  std::complex<unsigned> U_index_helmholtz;
237 
238 
239  /// Add the element's contribution to its residual vector.
240  /// flag=1(or 0): do (or don't) compute the contribution to the
241  /// Jacobian as well.
243  Vector<double>& residuals,
244  DenseMatrix<double>& jacobian,
245  const unsigned& flag);
246 
247 
248  /// Function pointer to the (global) prescribed-flux function
250 
251  /// The spatial dimension of the problem
252  unsigned Dim;
253  };
254 
255  /// ///////////////////////////////////////////////////////////////////
256  /// ///////////////////////////////////////////////////////////////////
257  /// ///////////////////////////////////////////////////////////////////
258 
259 
260  //===========================================================================
261  /// Constructor, takes the pointer to the "bulk" element, the
262  /// index of the fixed local coordinate and its value represented
263  /// by an integer (+/- 1), indicating that the face is located
264  /// at the max. or min. value of the "fixed" local coordinate
265  /// in the bulk element.
266  //===========================================================================
267  template<class ELEMENT>
269  FiniteElement* const& bulk_el_pt, const int& face_index)
270  : FaceGeometry<ELEMENT>(), FaceElement()
271  {
272  // Let the bulk element build the FaceElement, i.e. setup the pointers
273  // to its nodes (by referring to the appropriate nodes in the bulk
274  // element), etc.
275  bulk_el_pt->build_face_element(face_index, this);
276 
277 #ifdef PARANOID
278  {
279  // Check that the element is not a refineable 3d element
280  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
281  // If it's three-d
282  if (elem_pt->dim() == 3)
283  {
284  // Is it refineable
285  RefineableElement* ref_el_pt =
286  dynamic_cast<RefineableElement*>(elem_pt);
287  if (ref_el_pt != 0)
288  {
289  if (this->has_hanging_nodes())
290  {
291  throw OomphLibError("This flux element will not work correctly if "
292  "nodes are hanging\n",
293  OOMPH_CURRENT_FUNCTION,
294  OOMPH_EXCEPTION_LOCATION);
295  }
296  }
297  }
298  }
299 #endif
300 
301  // Initialise the prescribed-flux function pointer to zero
302  Flux_fct_pt = 0;
303 
304  // Extract the dimension of the problem from the dimension of
305  // the first node
306  Dim = this->node_pt(0)->ndim();
307 
308  // Set up U_index_helmholtz. Initialise to zero, which probably won't change
309  // in most cases, oh well, the price we pay for generality
310  U_index_helmholtz = std::complex<unsigned>(0, 1);
311 
312  // Cast to the appropriate HelmholtzEquation so that we can
313  // find the index at which the variable is stored
314  // We assume that the dimension of the full problem is the same
315  // as the dimension of the node, if this is not the case you will have
316  // to write custom elements, sorry
317  switch (Dim)
318  {
319  // One dimensional problem
320  case 1:
321  {
322  HelmholtzEquations<1>* eqn_pt =
323  dynamic_cast<HelmholtzEquations<1>*>(bulk_el_pt);
324  // If the cast has failed die
325  if (eqn_pt == 0)
326  {
327  std::string error_string =
328  "Bulk element must inherit from HelmholtzEquations.";
329  error_string +=
330  "Nodes are one dimensional, but cannot cast the bulk element to\n";
331  error_string += "HelmholtzEquations<1>\n.";
332  error_string += "If you desire this functionality, you must "
333  "implement it yourself\n";
334 
335  throw OomphLibError(
336  error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
337  }
338  // Otherwise read out the value
339  else
340  {
341  // Read the index from the (cast) bulk element
343  }
344  }
345  break;
346 
347  // Two dimensional problem
348  case 2:
349  {
350  HelmholtzEquations<2>* eqn_pt =
351  dynamic_cast<HelmholtzEquations<2>*>(bulk_el_pt);
352  // If the cast has failed die
353  if (eqn_pt == 0)
354  {
355  std::string error_string =
356  "Bulk element must inherit from HelmholtzEquations.";
357  error_string +=
358  "Nodes are two dimensional, but cannot cast the bulk element to\n";
359  error_string += "HelmholtzEquations<2>\n.";
360  error_string += "If you desire this functionality, you must "
361  "implement it yourself\n";
362 
363  throw OomphLibError(
364  error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
365  }
366  else
367  {
368  // Read the index from the (cast) bulk element
370  }
371  }
372 
373  break;
374 
375  // Three dimensional problem
376  case 3:
377  {
378  HelmholtzEquations<3>* eqn_pt =
379  dynamic_cast<HelmholtzEquations<3>*>(bulk_el_pt);
380  // If the cast has failed die
381  if (eqn_pt == 0)
382  {
383  std::string error_string =
384  "Bulk element must inherit from HelmholtzEquations.";
385  error_string += "Nodes are three dimensional, but cannot cast the "
386  "bulk element to\n";
387  error_string += "HelmholtzEquations<3>\n.";
388  error_string += "If you desire this functionality, you must "
389  "implement it yourself\n";
390 
391  throw OomphLibError(
392  error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
393  }
394  else
395  {
396  // Read the index from the (cast) bulk element
398  }
399  }
400  break;
401 
402  // Any other case is an error
403  default:
404  std::ostringstream error_stream;
405  error_stream << "Dimension of node is " << Dim
406  << ". It should be 1,2, or 3!" << std::endl;
407 
408  throw OomphLibError(
409  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
410  break;
411  }
412  }
413 
414 
415  //===========================================================================
416  /// Compute the element's residual vector and the (zero) Jacobian matrix.
417  //===========================================================================
418  template<class ELEMENT>
421  Vector<double>& residuals,
422  DenseMatrix<double>& jacobian,
423  const unsigned& flag)
424  {
425  // Find out how many nodes there are
426  const unsigned n_node = nnode();
427 
428  // Set up memory for the shape and test functions
429  Shape psif(n_node), testf(n_node);
430 
431  // Set the value of Nintpt
432  const unsigned n_intpt = integral_pt()->nweight();
433 
434  // Set the Vector to hold local coordinates
435  Vector<double> s(Dim - 1);
436 
437  // Integers to hold the local equation and unknown numbers
438  int local_eqn_real = 0, local_eqn_imag = 0;
439 
440  // Loop over the integration points
441  //--------------------------------
442  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
443  {
444  // Assign values of s
445  for (unsigned i = 0; i < (Dim - 1); i++)
446  {
447  s[i] = integral_pt()->knot(ipt, i);
448  }
449 
450  // Get the integral weight
451  double w = integral_pt()->weight(ipt);
452 
453  // Find the shape and test functions and return the Jacobian
454  // of the mapping
455  double J = shape_and_test(s, psif, testf);
456 
457  // Premultiply the weights and the Jacobian
458  double W = w * J;
459 
460  // Need to find position to feed into flux function, initialise to zero
461  Vector<double> interpolated_x(Dim, 0.0);
462 
463  // Calculate Eulerian position of integration point
464  for (unsigned l = 0; l < n_node; l++)
465  {
466  for (unsigned i = 0; i < Dim; i++)
467  {
468  interpolated_x[i] += nodal_position(l, i) * psif[l];
469  }
470  }
471 
472  // Get the imposed flux
473  std::complex<double> flux(0.0, 0.0);
474  get_flux(interpolated_x, flux);
475 
476  // Now add to the appropriate equations
477  // Loop over the test functions
478  for (unsigned l = 0; l < n_node; l++)
479  {
480  local_eqn_real = nodal_local_eqn(l, U_index_helmholtz.real());
481  /*IF it's not a boundary condition*/
482  if (local_eqn_real >= 0)
483  {
484  // Add the prescribed flux terms
485  residuals[local_eqn_real] -= flux.real() * testf[l] * W;
486 
487  // Imposed traction doesn't depend upon the solution,
488  // --> the Jacobian is always zero, so no Jacobian
489  // terms are required
490  }
491  local_eqn_imag = nodal_local_eqn(l, U_index_helmholtz.imag());
492  /*IF it's not a boundary condition*/
493  if (local_eqn_imag >= 0)
494  {
495  // Add the prescribed flux terms
496  residuals[local_eqn_imag] -= flux.imag() * testf[l] * W;
497 
498  // Imposed traction doesn't depend upon the solution,
499  // --> the Jacobian is always zero, so no Jacobian
500  // terms are required
501  }
502  }
503  }
504  }
505 
506 
507 } // namespace oomph
508 
509 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
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,...
Definition: elements.h:4497
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:5242
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:5328
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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:3050
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...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
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:5132
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:3220
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
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:227
A class for all isoparametric elements that solve the Helmholtz equations.
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
A class for elements that allow the imposition of an applied flux on the boundaries of Helmholtz elem...
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...
void output(std::ostream &outfile)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
void get_flux(const Vector< double > &x, std::complex< double > &flux)
Function to calculate the prescribed flux at a given spatial position.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
HelmholtzPrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
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 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...
unsigned Dim
The spatial dimension of the problem.
HelmholtzFluxElement()
Broken empty constructor.
HelmholtzFluxElement(const HelmholtzFluxElement &dummy)=delete
Broken copy constructor.
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 ...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void output(FILE *file_pt)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
HelmholtzPrescribedFluxFctPt & flux_fct_pt()
Broken assignment operator.
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.
std::complex< unsigned > U_index_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
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 std::complex< unsigned > u_index_helmholtz() const
Return the index at which the unknown value is stored.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
An OomphLibError object which should be thrown when an run-time error is encountered....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
//////////////////////////////////////////////////////////////////// ////////////////////////////////...