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