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-2023 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 
43 namespace 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
75  void operator=(
77  delete;
78 
79  /// Number of 'flux' terms for Z2 error estimation
80  unsigned num_Z2_flux_terms()
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();
177  this->Reaction_deriv_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:
215  : RefineableElement(),
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
231  void operator=(
233  DIM,
234  NNODE_1D>&) = delete;
235 
236  /// Number of continuously interpolated values: NREAGENT
237  unsigned ncont_interpolated_values() const
238  {
239  return NREAGENT;
240  }
241 
242  /// Number of vertex nodes in the element
243  unsigned nvertex_node() const
244  {
246  nvertex_node();
247  }
248 
249  /// Pointer to the j-th vertex node in the element
250  Node* vertex_node_pt(const unsigned& j) const
251  {
253  vertex_node_pt(j);
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.
261  unsigned nrecovery_order()
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...
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....
AdvectionDiffusionReactionReactionDerivFctPt Reaction_deriv_fct_pt
Pointer to reaction derivatives.
AdvectionDiffusionReactionReactionDerivFctPt & reaction_deriv_fct_pt()
Access function: Pointer to reaction derivatives function.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
AdvectionDiffusionReactionReactionFctPt & reaction_fct_pt()
Access function: Pointer to reaction function.
Vector< double > * Tau_pt
Pointer to global timescales.
AdvectionDiffusionReactionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
Vector< double > *& diff_pt()
Pointer to vector of diffusion coefficients.
AdvectionDiffusionReactionWindFctPt Wind_fct_pt
Pointer to wind function:
AdvectionDiffusionReactionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
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,...
Vector< double > *& tau_pt()
Pointer to vector of dimensionless timescales.
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....
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.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...