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-2022 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
42namespace 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
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
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
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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...
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...
HelmholtzPrescribedFluxFctPt & flux_fct_pt()
Broken assignment operator.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...