advection_diffusion_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 Advection Diffusion equations
28#ifndef OOMPH_ADV_DIFF_FLUX_ELEMENTS_HEADER
29#define OOMPH_ADV_DIFF_FLUX_ELEMENTS_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36// oomph-lib ncludes
37#include "../generic/Qelements.h"
38
39namespace oomph
40{
41 /// /////////////////////////////////////////////////////////////////////
42 /// /////////////////////////////////////////////////////////////////////
43 /// /////////////////////////////////////////////////////////////////////
44
45
46 //======================================================================
47 /// A class for elements that allow the imposition of an
48 /// applied flux on the boundaries of Advection Diffusion elements.
49 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
50 /// policy class.
51 //======================================================================
52 template<class ELEMENT>
53 class AdvectionDiffusionFluxElement : public virtual FaceGeometry<ELEMENT>,
54 public virtual FaceElement
55 {
56 public:
57 /// Function pointer to the prescribed-flux function fct(x,f(x)) --
58 /// x is a Vector!
60 const Vector<double>& x, double& flux);
61
62
63 /// Constructor, takes the pointer to the "bulk" element
64 /// and the index of the face to be created
66 const int& face_index);
67
68
69 /// Broken empty constructor
71 {
72 throw OomphLibError(
73 "Don't call empty constructor for AdvectionDiffusionFluxElement",
74 OOMPH_CURRENT_FUNCTION,
75 OOMPH_EXCEPTION_LOCATION);
76 }
77
78 /// Broken copy constructor
80 delete;
81
82 /// Broken assignment operator
84
85 /// Access function for the prescribed-flux function pointer
87 {
88 return Flux_fct_pt;
89 }
90
91
92 /// Add the element's contribution to its residual vector
94 {
95 // Call the generic residuals function with flag set to 0
96 // using a dummy matrix
99 }
100
101
102 /// Add the element's contribution to its residual vector and
103 /// its Jacobian matrix
105 DenseMatrix<double>& jacobian)
106 {
107 // Call the generic routine with the flag set to 1
109 residuals, jacobian, 1);
110 }
111
112 /// Specify the value of nodal zeta from the face geometry
113 /// The "global" intrinsic coordinate of the element when
114 /// viewed as part of a geometric object should be given by
115 /// the FaceElement representation, by default (needed to break
116 /// indeterminacy if bulk element is SolidElement)
117 double zeta_nodal(const unsigned& n,
118 const unsigned& k,
119 const unsigned& i) const
120 {
121 return FaceElement::zeta_nodal(n, k, i);
122 }
123
124 /// Output function -- forward to broken version in FiniteElement
125 /// until somebody decides what exactly they want to plot here...
126 void output(std::ostream& outfile)
127 {
128 FiniteElement::output(outfile);
129 }
130
131 /// Output function -- forward to broken version in FiniteElement
132 /// until somebody decides what exactly they want to plot here...
133 void output(std::ostream& outfile, const unsigned& nplot)
134 {
135 FiniteElement::output(outfile, nplot);
136 }
137
138
139 protected:
140 /// Function to compute the shape and test functions and to return
141 /// the Jacobian of mapping between local and global (Eulerian)
142 /// coordinates
143 inline double shape_and_test(const Vector<double>& s,
144 Shape& psi,
145 Shape& test) const
146 {
147 // Find number of nodes
148 unsigned n_node = nnode();
149
150 // Get the shape functions
151 shape(s, psi);
152
153 // Set the test functions to be the same as the shape functions
154 for (unsigned i = 0; i < n_node; i++)
155 {
156 test[i] = psi[i];
157 }
158
159 // Return the value of the jacobian
160 return J_eulerian(s);
161 }
162
163
164 /// Function to compute the shape and test functions and to return
165 /// the Jacobian of mapping between local and global (Eulerian)
166 /// coordinates
167 inline double shape_and_test_at_knot(const unsigned& ipt,
168 Shape& psi,
169 Shape& test) const
170 {
171 // Find number of nodes
172 unsigned n_node = nnode();
173
174 // Get the shape functions
175 shape_at_knot(ipt, psi);
176
177 // Set the test functions to be the same as the shape functions
178 for (unsigned i = 0; i < n_node; i++)
179 {
180 test[i] = psi[i];
181 }
182
183 // Return the value of the jacobian
184 return J_eulerian_at_knot(ipt);
185 }
186
187
188 /// Function to calculate the prescribed flux at a given spatial
189 /// position
190 void get_flux(const Vector<double>& x, double& flux)
191 {
192 // If the function pointer is zero return zero
193 if (Flux_fct_pt == 0)
194 {
195 flux = 0.0;
196 }
197 // Otherwise call the function
198 else
199 {
200 (*Flux_fct_pt)(x, flux);
201 }
202 }
203
204 private:
205 /// Add the element's contribution to its residual vector.
206 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
208 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
209
210
211 /// Function pointer to the (global) prescribed-flux function
213
214 /// The spatial dimension of the problem
215 unsigned Dim;
216
217 /// The index at which the unknown is stored at the nodes
219 };
220
221
222 /// ////////////////////////////////////////////////////////////////////
223 /// ////////////////////////////////////////////////////////////////////
224 /// ////////////////////////////////////////////////////////////////////
225
226
227 //===========================================================================
228 /// Constructor, takes the pointer to the "bulk" element and the index
229 /// of the face to be created
230 //===========================================================================
231 template<class ELEMENT>
233 FiniteElement* const& bulk_el_pt, const int& face_index)
234 : FaceGeometry<ELEMENT>(), FaceElement()
235 {
236 // Let the bulk element build the FaceElement, i.e. setup the pointers
237 // to its nodes (by referring to the appropriate nodes in the bulk
238 // element), etc.
239 bulk_el_pt->build_face_element(face_index, this);
240
241#ifdef PARANOID
242 {
243 // Check that the element is not a refineable 3d element
244 ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
245 // If it's three-d
246 if (elem_pt->dim() == 3)
247 {
248 // Is it refineable
249 RefineableElement* ref_el_pt =
250 dynamic_cast<RefineableElement*>(elem_pt);
251 if (ref_el_pt != 0)
252 {
253 if (this->has_hanging_nodes())
254 {
255 throw OomphLibError("This flux element will not work correctly if "
256 "nodes are hanging\n",
257 OOMPH_CURRENT_FUNCTION,
258 OOMPH_EXCEPTION_LOCATION);
259 }
260 }
261 }
262 }
263#endif
264
265 // Initialise the prescribed-flux function pointer to zero
266 Flux_fct_pt = 0;
267
268 // Extract the dimension of the problem from the dimension of
269 // the first node
270 Dim = this->node_pt(0)->ndim();
271
272
273 // Set up U_index_adv_diff. Initialise to zero, which probably won't change
274 // in most cases, oh well, the price we pay for generality
276
277 // Cast to the appropriate AdvectionDiffusionEquation so that we can
278 // find the index at which the variable is stored
279 // We assume that the dimension of the full problem is the same
280 // as the dimension of the node, if this is not the case you will have
281 // to write custom elements, sorry
282 switch (Dim)
283 {
284 // One dimensional problem
285 case 1:
286 {
288 dynamic_cast<AdvectionDiffusionEquations<1>*>(bulk_el_pt);
289 // If the cast has failed die
290 if (eqn_pt == 0)
291 {
292 std::string error_string =
293 "Bulk element must inherit from AdvectionDiffusionEquations.";
294 error_string +=
295 "Nodes are one dimensional, but cannot cast the bulk element to\n";
296 error_string += "AdvectionDiffusionEquations<1>\n.";
297 error_string += "If you desire this functionality, you must "
298 "implement it yourself\n";
299
300 throw OomphLibError(
301 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
302 }
303 // Otherwise read out the value
304 else
305 {
306 // Read the index from the (cast) bulk element
308 }
309 }
310 break;
311
312 // Two dimensional problem
313 case 2:
314 {
316 dynamic_cast<AdvectionDiffusionEquations<2>*>(bulk_el_pt);
317 // If the cast has failed die
318 if (eqn_pt == 0)
319 {
320 std::string error_string =
321 "Bulk element must inherit from AdvectionDiffusionEquations.";
322 error_string +=
323 "Nodes are two dimensional, but cannot cast the bulk element to\n";
324 error_string += "AdvectionDiffusionEquations<2>\n.";
325 error_string += "If you desire this functionality, you must "
326 "implement it yourself\n";
327
328 throw OomphLibError(
329 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
330 }
331 else
332 {
333 // Read the index from the (cast) bulk element.
335 }
336 }
337 break;
338
339 // Three dimensional problem
340 case 3:
341 {
343 dynamic_cast<AdvectionDiffusionEquations<3>*>(bulk_el_pt);
344 // If the cast has failed die
345 if (eqn_pt == 0)
346 {
347 std::string error_string =
348 "Bulk element must inherit from AdvectionDiffusionEquations.";
349 error_string += "Nodes are three dimensional, but cannot cast the "
350 "bulk element to\n";
351 error_string += "AdvectionDiffusionEquations<3>\n.";
352 error_string += "If you desire this functionality, you must "
353 "implement it yourself\n";
354
355 throw OomphLibError(
356 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
357 }
358 else
359 {
360 // Read the index from the (cast) bulk element.
362 }
363 }
364 break;
365
366 // Any other case is an error
367 default:
368 std::ostringstream error_stream;
369 error_stream << "Dimension of node is " << Dim
370 << ". It should be 1,2, or 3!" << std::endl;
371
372 throw OomphLibError(
373 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
374 break;
375 }
376 }
377
378
379 //===========================================================================
380 /// Compute the element's residual vector and the (zero) Jacobian matrix.
381 //===========================================================================
382 template<class ELEMENT>
385 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
386 {
387 // Find out how many nodes there are
388 const unsigned n_node = nnode();
389
390 // Set up memory for the shape and test functions
391 Shape psif(n_node), testf(n_node);
392
393 // Set the value of n_intpt
394 const unsigned n_intpt = integral_pt()->nweight();
395
396 // Set the Vector to hold local coordinates
397 Vector<double> s(Dim - 1);
398
399 // Integers used to store the local equation number and local unknown
400 // indices for the residuals and jacobians
401 int local_eqn = 0;
402
403 // Locally cache the index at which the variable is stored
404 const unsigned u_index_adv_diff = U_index_adv_diff;
405
406
407 // Loop over the integration points
408 //--------------------------------
409 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
410 {
411 // Assign values of s
412 for (unsigned i = 0; i < (Dim - 1); i++)
413 {
414 s[i] = integral_pt()->knot(ipt, i);
415 }
416
417 // Get the integral weight
418 double w = integral_pt()->weight(ipt);
419
420 // Find the shape and test functions and return the Jacobian
421 // of the mapping
422 double J = shape_and_test(s, psif, testf);
423
424 // Premultiply the weights and the Jacobian
425 double W = w * J;
426
427 // Need to find position to feed into flux function
428 Vector<double> interpolated_x(Dim, 0.0);
429
430 // Calculate position
431 for (unsigned l = 0; l < n_node; l++)
432 {
433 // Loop over coordinate directions
434 for (unsigned i = 0; i < Dim; i++)
435 {
436 interpolated_x[i] += nodal_position(l, i) * psif[l];
437 }
438 }
439
440 // Get the imposed flux
441 double flux;
442 get_flux(interpolated_x, flux);
443
444 // Now add to the appropriate equations
445
446 // Loop over the test functions
447 for (unsigned l = 0; l < n_node; l++)
448 {
449 // Set the local equation number
450 local_eqn = nodal_local_eqn(l, u_index_adv_diff);
451 /*IF it's not a boundary condition*/
452 if (local_eqn >= 0)
453 {
454 // Add the prescribed flux terms
455 residuals[local_eqn] += flux * testf[l] * W;
456
457 // Imposed traction doesn't depend upon the solution,
458 // --> the Jacobian is always zero, so no Jacobian
459 // terms are required
460 }
461 }
462 }
463 }
464
465
466} // namespace oomph
467
468#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class for all elements that solve the Advection Diffusion equations using isoparametric elements.
virtual unsigned u_index_adv_diff() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void output(std::ostream &outfile)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
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 ...
AdvectionDiffusionPrescribedFluxFctPt & flux_fct_pt()
Access function for the prescribed-flux function pointer.
AdvectionDiffusionFluxElement(const AdvectionDiffusionFluxElement &dummy)=delete
Broken copy constructor.
void operator=(const AdvectionDiffusionFluxElement &)=delete
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.
void output(std::ostream &outfile, const unsigned &nplot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
unsigned U_index_adv_diff
The index at which the unknown is stored at the nodes.
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.
AdvectionDiffusionPrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
unsigned Dim
The spatial dimension of the problem.
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 get_flux(const Vector< double > &x, double &flux)
Function to calculate the prescribed flux at a given spatial position.
void fill_in_generic_residual_contribution_adv_diff_flux(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the Jacobi...
void(* AdvectionDiffusionPrescribedFluxFctPt)(const Vector< double > &x, double &flux)
Function pointer to the prescribed-flux function fct(x,f(x)) – x is a Vector!
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
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...