refineable_advection_diffusion_reaction_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 solve the advection diffusion equation
27// and that can be refined.
28
29#ifndef OOMPH_REFINEABLE_ADVECTION_DIFFUSION_REACTION_ELEMENTS_HEADER
30#define OOMPH_REFINEABLE_ADVECTION_DIFFUSION_REACTION_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37// oomph-lib headers
38#include "../generic/refineable_quad_element.h"
39#include "../generic/refineable_brick_element.h"
40#include "../generic/error_estimator.h"
42
43namespace oomph
44{
45 //======================================================================
46 /// A version of the Advection Diffusion Reaction equations that can be
47 /// used with non-uniform mesh refinement. In essence, the class overloads
48 /// the fill_in_generic_residual_contribution_adv_diff_react()
49 /// function so that contributions
50 /// from hanging nodes (or alternatively in-compatible function values)
51 /// are taken into account.
52 //======================================================================
53 template<unsigned NREAGENT, unsigned DIM>
55 : public virtual AdvectionDiffusionReactionEquations<NREAGENT, DIM>,
56 public virtual RefineableElement,
57 public virtual ElementWithZ2ErrorEstimator
58 {
59 public:
60 /// Empty Constructor
62 : AdvectionDiffusionReactionEquations<NREAGENT, DIM>(),
65 {
66 }
67
68
69 /// Broken copy constructor
72 dummy) = delete;
73
74 /// Broken assignment operator
77 delete;
78
79 /// Number of 'flux' terms for Z2 error estimation
81 {
82 return NREAGENT * DIM;
83 }
84
85 /// Get 'flux' for Z2 error recovery:
86 /// Standard flux.from AdvectionDiffusionReaction equations
88 {
89 this->get_flux(s, flux);
90 }
91
92
93 /// Get the function values c in Vector.
94 /// Note: Given the generality of the interface (this function
95 /// is usually called from black-box documentation or interpolation
96 /// routines), the values Vector sets its own size in here.
98 Vector<double>& values)
99 {
100 // Set size of Vector: c
101 values.resize(NREAGENT);
102
103 // Find number of nodes
104 unsigned n_node = nnode();
105
106 // Local shape function
107 Shape psi(n_node);
108
109 // Find values of shape function
110 shape(s, psi);
111
112 // Loop over the unknowns
113 for (unsigned r = 0; r < NREAGENT; r++)
114 {
115 unsigned c_nodal_index = this->c_index_adv_diff_react(r);
116
117 // Initialise value of c
118 values[r] = 0.0;
119
120 // Loop over the local nodes and sum
121 for (unsigned l = 0; l < n_node; l++)
122 {
123 values[r] += this->nodal_value(l, c_nodal_index) * psi[l];
124 }
125 }
126 }
127
128 /// Get the function values c in Vector.
129 /// Note: Given the generality of the interface (this function
130 /// is usually called from black-box documentation or interpolation
131 /// routines), the values Vector sets its own size in here.
132 void get_interpolated_values(const unsigned& t,
133 const Vector<double>& s,
134 Vector<double>& values)
135 {
136 // Set size of Vector:
137 values.resize(NREAGENT);
138
139 // Find out how many nodes there are
140 const unsigned n_node = nnode();
141
142 // Shape functions
143 Shape psi(n_node);
144 shape(s, psi);
145
146 // Loop over the reagents
147 for (unsigned r = 0; r < NREAGENT; r++)
148 {
149 // Find the nodal index at which the unknown is stored
150 unsigned c_nodal_index = this->c_index_adv_diff_react(r);
151
152 // Initialise
153 values[r] = 0.0;
154
155 // Calculate value
156 for (unsigned l = 0; l < n_node; l++)
157 {
158 values[r] += this->nodal_value(t, l, c_nodal_index) * psi[l];
159 }
160 }
161 }
162
163
164 /// Further build: Copy all pointers from the father
165 /// element
167 {
169 cast_father_element_pt = dynamic_cast<
171 this->father_element_pt());
172
173 // Set the values of the pointers from the father
174 this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
175 this->Wind_fct_pt = cast_father_element_pt->wind_fct_pt();
176 this->Reaction_fct_pt = cast_father_element_pt->reaction_fct_pt();
178 cast_father_element_pt->reaction_deriv_fct_pt();
179
180 this->Diff_pt = cast_father_element_pt->diff_pt();
181 this->Tau_pt = cast_father_element_pt->tau_pt();
182
183 // Set the value of the ALE flag
184 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
185 }
186
187 protected:
188 /// Add the element's contribution to the elemental residual vector
189 /// and/or Jacobian matrix
190 /// flag=1: compute both
191 /// flag=0: compute only residual vector
193 Vector<double>& residuals,
194 DenseMatrix<double>& jacobian,
195 DenseMatrix<double>& mass_matrix,
196 unsigned flag);
197 };
198
199
200 //======================================================================
201 /// Refineable version of QAdvectionDiffusionReactionElement.
202 /// Inherit from the standard QAdvectionDiffusionReactionElement and the
203 /// appropriate refineable geometric element and the refineable equations.
204 //======================================================================
205 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
207 : public QAdvectionDiffusionReactionElement<NREAGENT, DIM, NNODE_1D>,
208 public virtual RefineableAdvectionDiffusionReactionEquations<NREAGENT,
209 DIM>,
210 public virtual RefineableQElement<DIM>
211 {
212 public:
213 /// Empty Constructor:
217 RefineableQElement<DIM>(),
218 QAdvectionDiffusionReactionElement<NREAGENT, DIM, NNODE_1D>()
219 {
220 }
221
222
223 /// Broken copy constructor
226 DIM,
227 NNODE_1D>& dummy) =
228 delete;
229
230 /// Broken assignment operator
233 DIM,
234 NNODE_1D>&) = delete;
235
236 /// Number of continuously interpolated values: NREAGENT
238 {
239 return NREAGENT;
240 }
241
242 /// Number of vertex nodes in the element
243 unsigned nvertex_node() const
244 {
247 }
248
249 /// Pointer to the j-th vertex node in the element
250 Node* vertex_node_pt(const unsigned& j) const
251 {
254 }
255
256 /// Rebuild from sons: empty
257 void rebuild_from_sons(Mesh*& mesh_pt) {}
258
259 /// Order of recovery shape functions for Z2 error estimation:
260 /// Same order as shape functions.
262 {
263 return (NNODE_1D - 1);
264 }
265
266 /// Perform additional hanging node procedures for variables
267 /// that are not interpolated by all nodes. Empty.
269 };
270
271 /// /////////////////////////////////////////////////////////////////////
272 /// /////////////////////////////////////////////////////////////////////
273 /// /////////////////////////////////////////////////////////////////////
274
275
276 //=======================================================================
277 /// Face geometry for the RefineableQuadAdvectionDiffusionReactionElement
278 /// elements: The spatial dimension of the face elements is one lower than
279 /// that of the bulk element but they have the same number of points along
280 /// their 1D edges.
281 //=======================================================================
282 template<unsigned NREAGENT, unsigned DIM, unsigned NNODE_1D>
284 RefineableQAdvectionDiffusionReactionElement<NREAGENT, DIM, NNODE_1D>>
285 : public virtual QElement<DIM - 1, NNODE_1D>
286 {
287 public:
288 /// Constructor: Call the constructor for the
289 /// appropriate lower-dimensional QElement
290 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
291 };
292
293} // namespace oomph
294
295#endif
static char t char * s
Definition: cfortran.h:568
char t
Definition: cfortran.h:568
A class for all elements that solve the Advection Diffusion Reaction equations using isoparametric el...
Vector< double > *& diff_pt()
Pointer to vector of diffusion coefficients.
AdvectionDiffusionReactionSourceFctPt Source_fct_pt
Pointer to source function:
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
AdvectionDiffusionReactionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
AdvectionDiffusionReactionReactionDerivFctPt Reaction_deriv_fct_pt
Pointer to reaction derivatives.
AdvectionDiffusionReactionReactionFctPt & reaction_fct_pt()
Access function: Pointer to reaction function.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
Vector< double > * Tau_pt
Pointer to global timescales.
AdvectionDiffusionReactionReactionDerivFctPt & reaction_deriv_fct_pt()
Access function: Pointer to reaction derivatives function.
AdvectionDiffusionReactionWindFctPt Wind_fct_pt
Pointer to wind function:
Vector< double > *& tau_pt()
Pointer to vector of dimensionless timescales.
virtual unsigned c_index_adv_diff_react(const unsigned &i) const
Return the index at which the unknown i-th reagent is stored. The default value, i,...
AdvectionDiffusionReactionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
AdvectionDiffusionReactionReactionFctPt Reaction_fct_pt
Pointer to reaction function.
Vector< double > * Diff_pt
Pointer to global diffusion coefficients.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2491
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 Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2500
A general mesh class.
Definition: mesh.h:67
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
A version of the Advection Diffusion Reaction equations that can be used with non-uniform mesh refine...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
RefineableAdvectionDiffusionReactionEquations(const RefineableAdvectionDiffusionReactionEquations< NREAGENT, DIM > &dummy)=delete
Broken copy constructor.
void operator=(const RefineableAdvectionDiffusionReactionEquations< NREAGENT, DIM > &)=delete
Broken assignment operator.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function values c in Vector. Note: Given the generality of the interface (this function is us...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Standard flux.from AdvectionDiffusionReaction equations.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function values c in Vector. Note: Given the generality of the interface (this function is us...
void fill_in_generic_residual_contribution_adv_diff_react(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add the element's contribution to the elemental residual vector and/or Jacobian matrix flag=1: comput...
void further_build()
Further build: Copy all pointers from the father element.
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Refineable version of QAdvectionDiffusionReactionElement. Inherit from the standard QAdvectionDiffusi...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: NREAGENT.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void operator=(const RefineableQAdvectionDiffusionReactionElement< NREAGENT, DIM, NNODE_1D > &)=delete
Broken assignment operator.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
RefineableQAdvectionDiffusionReactionElement(const RefineableQAdvectionDiffusionReactionElement< NREAGENT, DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...