poisson_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 Poisson equations
28#ifndef OOMPH_POISSON_FLUX_ELEMENTS_HEADER
29#define OOMPH_POISSON_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 includes
38#include "../generic/Qelements.h"
39
40namespace oomph
41{
42 //======================================================================
43 /// A class for elements that allow the imposition of an
44 /// applied flux on the boundaries of Poisson elements.
45 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
46 /// policy class.
47 //======================================================================
48 template<class ELEMENT>
49 class PoissonFluxElement : public virtual FaceGeometry<ELEMENT>,
50 public virtual FaceElement
51 {
52 public:
53 /// Function pointer to the prescribed-flux function fct(x,f(x)) --
54 /// x is a Vector!
56 double& flux);
57
58 /// Constructor, takes the pointer to the "bulk" element and the
59 /// index of the face to which the element is attached.
60 PoissonFluxElement(FiniteElement* const& bulk_el_pt, const int& face_index);
61
62 /// Broken empty constructor
64 {
65 throw OomphLibError("Don't call empty constructor for PoissonFluxElement",
66 OOMPH_CURRENT_FUNCTION,
67 OOMPH_EXCEPTION_LOCATION);
68 }
69
70 /// Broken copy constructor
71 PoissonFluxElement(const PoissonFluxElement& dummy) = delete;
72
73 /// Broken assignment operator
74 void operator=(const PoissonFluxElement&) = delete;
75
76 /// Specify the value of nodal zeta from the face geometry
77 /// The "global" intrinsic coordinate of the element when
78 /// viewed as part of a geometric object should be given by
79 /// the FaceElement representation, by default (needed to break
80 /// indeterminacy if bulk element is SolidElement)
81 double zeta_nodal(const unsigned& n,
82 const unsigned& k,
83 const unsigned& i) const
84 {
85 return FaceElement::zeta_nodal(n, k, i);
86 }
87
88 /// Access function for the prescribed-flux function pointer
90 {
91 return Flux_fct_pt;
92 }
93
94
95 /// Add the element's contribution to its residual vector
97 {
98 // Call the generic residuals function with flag set to 0
99 // using a dummy matrix argument
102 }
103
104
105 /// Add the element's contribution to its residual vector and its
106 /// Jacobian matrix
108 DenseMatrix<double>& jacobian)
109 {
110 // Call the generic routine with the flag set to 1
112 residuals, jacobian, 1);
113 }
114
115 /// Output function
116 void output(std::ostream& outfile)
117 {
118 const unsigned n_plot = 5;
119 output(outfile, n_plot);
120 }
121
122 /// Output function
123 void output(std::ostream& outfile, const unsigned& nplot)
124 {
125 // Dimension of element
126 unsigned el_dim = dim();
127
128 // Vector of local coordinates
129 Vector<double> s(el_dim);
130
131 // Tecplot header info
132 outfile << tecplot_zone_string(nplot);
133
134 // Loop over plot points
135 unsigned num_plot_points = nplot_points(nplot);
136 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
137 {
138 // Get local coordinates of plot point
139 get_s_plot(iplot, nplot, s);
140
141 Vector<double> x(el_dim + 1);
142 for (unsigned i = 0; i < el_dim + 1; i++)
143 {
144 x[i] = interpolated_x(s, i);
145 outfile << x[i] << " ";
146 }
147 double flux = 0.0;
148 get_flux(x, flux);
149 outfile << flux << std::endl;
150 }
151
152 // Write tecplot footer (e.g. FE connectivity lists)
153 write_tecplot_zone_footer(outfile, nplot);
154 }
155
156
157 /// C-style output function -- forward to broken version in FiniteElement
158 /// until somebody decides what exactly they want to plot here...
159 void output(FILE* file_pt)
160 {
161 FiniteElement::output(file_pt);
162 }
163
164 /// C-style output function -- forward to broken version in
165 /// FiniteElement until somebody decides what exactly they want to plot
166 /// here...
167 void output(FILE* file_pt, const unsigned& n_plot)
168 {
169 FiniteElement::output(file_pt, n_plot);
170 }
171
172
173 protected:
174 /// Function to compute the shape and test functions and to return
175 /// the Jacobian of mapping between local and global (Eulerian)
176 /// coordinates
177 inline double shape_and_test(const Vector<double>& s,
178 Shape& psi,
179 Shape& test) const
180 {
181 // Find number of nodes
182 unsigned n_node = nnode();
183
184 // Get the shape functions
185 shape(s, psi);
186
187 // Set the test functions to be the same as the shape functions
188 for (unsigned i = 0; i < n_node; i++)
189 {
190 test[i] = psi[i];
191 }
192
193 // Return the value of the jacobian
194 return J_eulerian(s);
195 }
196
197
198 /// Function to compute the shape and test functions and to return
199 /// the Jacobian of mapping between local and global (Eulerian)
200 /// coordinates
201 inline double shape_and_test_at_knot(const unsigned& ipt,
202 Shape& psi,
203 Shape& test) const
204 {
205 // Find number of nodes
206 unsigned n_node = nnode();
207
208 // Get the shape functions
209 shape_at_knot(ipt, psi);
210
211 // Set the test functions to be the same as the shape functions
212 for (unsigned i = 0; i < n_node; i++)
213 {
214 test[i] = psi[i];
215 }
216
217 // Return the value of the jacobian
218 return J_eulerian_at_knot(ipt);
219 }
220
221
222 /// Function to calculate the prescribed flux at a given spatial
223 /// position
224 void get_flux(const Vector<double>& x, double& flux)
225 {
226 // If the function pointer is zero return zero
227 if (Flux_fct_pt == 0)
228 {
229 flux = 0.0;
230 }
231 // Otherwise call the function
232 else
233 {
234 (*Flux_fct_pt)(x, flux);
235 }
236 }
237
238 private:
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 /// The index at which the unknown is stored at the nodes
256 };
257
258 /// ///////////////////////////////////////////////////////////////////
259 /// ///////////////////////////////////////////////////////////////////
260 /// ///////////////////////////////////////////////////////////////////
261
262
263 //===========================================================================
264 /// Constructor, takes the pointer to the "bulk" element, the
265 /// index of the fixed local coordinate and its value represented
266 /// by an integer (+/- 1), indicating that the face is located
267 /// at the max. or min. value of the "fixed" local coordinate
268 /// in the bulk element.
269 //===========================================================================
270 template<class ELEMENT>
272 FiniteElement* const& bulk_el_pt, const int& face_index)
273 : FaceGeometry<ELEMENT>(), FaceElement()
274 {
275 // Let the bulk element build the FaceElement, i.e. setup the pointers
276 // to its nodes (by referring to the appropriate nodes in the bulk
277 // element), etc.
278 bulk_el_pt->build_face_element(face_index, this);
279
280#ifdef PARANOID
281 {
282 // Check that the element is not a refineable 3d element
283 ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
284 // If it's three-d
285 if (elem_pt->dim() == 3)
286 {
287 // Is it refineable
288 RefineableElement* ref_el_pt =
289 dynamic_cast<RefineableElement*>(elem_pt);
290 if (ref_el_pt != 0)
291 {
292 if (this->has_hanging_nodes())
293 {
294 throw OomphLibError("This flux element will not work correctly if "
295 "nodes are hanging\n",
296 OOMPH_CURRENT_FUNCTION,
297 OOMPH_EXCEPTION_LOCATION);
298 }
299 }
300 }
301 }
302#endif
303
304 // Initialise the prescribed-flux function pointer to zero
305 Flux_fct_pt = 0;
306
307 // Extract the dimension of the problem from the dimension of
308 // the first node
309 Dim = this->node_pt(0)->ndim();
310
311 // Set up U_index_poisson. Initialise to zero, which probably won't change
312 // in most cases, oh well, the price we pay for generality
313 U_index_poisson = 0;
314
315 // Cast to the appropriate PoissonEquation so that we can
316 // find the index at which the variable is stored
317 // We assume that the dimension of the full problem is the same
318 // as the dimension of the node, if this is not the case you will have
319 // to write custom elements, sorry
320 switch (Dim)
321 {
322 // One dimensional problem
323 case 1:
324 {
325 PoissonEquations<1>* eqn_pt =
326 dynamic_cast<PoissonEquations<1>*>(bulk_el_pt);
327 // If the cast has failed die
328 if (eqn_pt == 0)
329 {
330 std::string error_string =
331 "Bulk element must inherit from PoissonEquations.";
332 error_string +=
333 "Nodes are one dimensional, but cannot cast the bulk element to\n";
334 error_string += "PoissonEquations<1>\n.";
335 error_string += "If you desire this functionality, you must "
336 "implement it yourself\n";
337
338 throw OomphLibError(
339 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
340 }
341 // Otherwise read out the value
342 else
343 {
344 // Read the index from the (cast) bulk element
346 }
347 }
348 break;
349
350 // Two dimensional problem
351 case 2:
352 {
353 PoissonEquations<2>* eqn_pt =
354 dynamic_cast<PoissonEquations<2>*>(bulk_el_pt);
355 // If the cast has failed die
356 if (eqn_pt == 0)
357 {
358 std::string error_string =
359 "Bulk element must inherit from PoissonEquations.";
360 error_string +=
361 "Nodes are two dimensional, but cannot cast the bulk element to\n";
362 error_string += "PoissonEquations<2>\n.";
363 error_string += "If you desire this functionality, you must "
364 "implement it yourself\n";
365
366 throw OomphLibError(
367 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
368 }
369 else
370 {
371 // Read the index from the (cast) bulk element.
373 }
374 }
375 break;
376
377 // Three dimensional problem
378 case 3:
379 {
380 PoissonEquations<3>* eqn_pt =
381 dynamic_cast<PoissonEquations<3>*>(bulk_el_pt);
382 // If the cast has failed die
383 if (eqn_pt == 0)
384 {
385 std::string error_string =
386 "Bulk element must inherit from PoissonEquations.";
387 error_string += "Nodes are three dimensional, but cannot cast the "
388 "bulk element to\n";
389 error_string += "PoissonEquations<3>\n.";
390 error_string += "If you desire this functionality, you must "
391 "implement it yourself\n";
392
393 throw OomphLibError(
394 error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
395 }
396 else
397 {
398 // Read the index from the (cast) bulk element.
400 }
401 }
402 break;
403
404 // Any other case is an error
405 default:
406 std::ostringstream error_stream;
407 error_stream << "Dimension of node is " << Dim
408 << ". It should be 1,2, or 3!" << std::endl;
409
410 throw OomphLibError(
411 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
412 break;
413 }
414 }
415
416
417 //===========================================================================
418 /// Compute the element's residual vector and the (zero) Jacobian matrix.
419 //===========================================================================
420 template<class ELEMENT>
423 Vector<double>& residuals,
424 DenseMatrix<double>& jacobian,
425 const unsigned& flag)
426 {
427 // Find out how many nodes there are
428 const unsigned n_node = nnode();
429
430 // Set up memory for the shape and test functions
431 Shape psif(n_node), testf(n_node);
432
433 // Set the value of Nintpt
434 const unsigned n_intpt = integral_pt()->nweight();
435
436 // Set the Vector to hold local coordinates
437 Vector<double> s(Dim - 1);
438
439 // Integers to hold the local equation and unknown numbers
440 int local_eqn = 0;
441
442 // Locally cache the index at which the variable is stored
443 const unsigned u_index_poisson = U_index_poisson;
444
445 // Loop over the integration points
446 //--------------------------------
447 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
448 {
449 // Assign values of s
450 for (unsigned i = 0; i < (Dim - 1); i++)
451 {
452 s[i] = integral_pt()->knot(ipt, i);
453 }
454
455 // Get the integral weight
456 double w = integral_pt()->weight(ipt);
457
458 // Find the shape and test functions and return the Jacobian
459 // of the mapping
460 double J = shape_and_test(s, psif, testf);
461
462 // Premultiply the weights and the Jacobian
463 double W = w * J;
464
465 // Need to find position to feed into flux function, initialise to zero
466 Vector<double> interpolated_x(Dim, 0.0);
467
468 // Calculate velocities and derivatives
469 for (unsigned l = 0; l < n_node; l++)
470 {
471 // Loop over velocity components
472 for (unsigned i = 0; i < Dim; i++)
473 {
474 interpolated_x[i] += nodal_position(l, i) * psif[l];
475 }
476 }
477
478 // Get the imposed flux
479 double flux;
480 get_flux(interpolated_x, flux);
481
482 // Now add to the appropriate equations
483
484 // Loop over the test functions
485 for (unsigned l = 0; l < n_node; l++)
486 {
487 local_eqn = nodal_local_eqn(l, u_index_poisson);
488 /*IF it's not a boundary condition*/
489 if (local_eqn >= 0)
490 {
491 // Add the prescribed flux terms
492 residuals[local_eqn] -= flux * testf[l] * W;
493
494 // Imposed traction doesn't depend upon the solution,
495 // --> the Jacobian is always zero, so no Jacobian
496 // terms are required
497 }
498 }
499 }
500 }
501
502
503} // namespace oomph
504
505#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 interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
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 std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
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 dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
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 write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
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....
A class for all isoparametric elements that solve the Poisson equations.
virtual unsigned u_index_poisson() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
A class for elements that allow the imposition of an applied flux on the boundaries of Poisson elemen...
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 output(FILE *file_pt)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
PoissonPrescribedFluxFctPt & flux_fct_pt()
Access function for the prescribed-flux function pointer.
PoissonFluxElement()
Broken empty constructor.
PoissonFluxElement(const PoissonFluxElement &dummy)=delete
Broken copy constructor.
unsigned Dim
The spatial dimension of the problem.
void operator=(const PoissonFluxElement &)=delete
Broken assignment operator.
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 ...
void get_flux(const Vector< double > &x, double &flux)
Function to calculate the prescribed flux at a given spatial position.
void(* PoissonPrescribedFluxFctPt)(const Vector< double > &x, double &flux)
Function pointer to the prescribed-flux function fct(x,f(x)) – x is a Vector!
void fill_in_generic_residual_contribution_poisson_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 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.
unsigned U_index_poisson
The index at which the unknown is stored at the nodes.
void output(std::ostream &outfile)
Output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
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...
PoissonPrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...