refineable_advection_diffusion_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_ELEMENTS_HEADER
30#define OOMPH_REFINEABLE_ADVECTION_DIFFUSION_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 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()
49 /// function so that contributions
50 /// from hanging nodes (or alternatively in-compatible function values)
51 /// are taken into account.
52 //======================================================================
53 template<unsigned DIM>
55 : public virtual AdvectionDiffusionEquations<DIM>,
56 public virtual RefineableElement,
57 public virtual ElementWithZ2ErrorEstimator
58 {
59 public:
60 /// Empty Constructor
65 {
66 }
67
68
69 /// Broken copy constructor
71 const RefineableAdvectionDiffusionEquations<DIM>& dummy) = delete;
72
73 /// Broken assignment operator
75
76 /// Number of 'flux' terms for Z2 error estimation
78 {
79 return DIM;
80 }
81
82 /// Get 'flux' for Z2 error recovery:
83 /// Standard flux.from AdvectionDiffusion equations
85 {
86 this->get_flux(s, flux);
87 }
88
89
90 /// Get the function value u in Vector.
91 /// Note: Given the generality of the interface (this function
92 /// is usually called from black-box documentation or interpolation
93 /// routines), the values Vector sets its own size in here.
95 Vector<double>& values)
96 {
97 // Set size of Vector: u
98 values.resize(1);
99
100 // Find number of nodes
101 unsigned n_node = nnode();
102
103 // Find the index at which the unknown is stored
104 unsigned u_nodal_index = this->u_index_adv_diff();
105
106 // Local shape function
107 Shape psi(n_node);
108
109 // Find values of shape function
110 shape(s, psi);
111
112 // Initialise value of u
113 values[0] = 0.0;
114
115 // Loop over the local nodes and sum
116 for (unsigned l = 0; l < n_node; l++)
117 {
118 values[0] += this->nodal_value(l, u_nodal_index) * psi[l];
119 }
120 }
121
122 /// Get the function value u in Vector.
123 /// Note: Given the generality of the interface (this function
124 /// is usually called from black-box documentation or interpolation
125 /// routines), the values Vector sets its own size in here.
126 void get_interpolated_values(const unsigned& t,
127 const Vector<double>& s,
128 Vector<double>& values)
129 {
130 // Set size of Vector: u
131 values.resize(1);
132
133 // Find number of nodes
134 const unsigned n_node = nnode();
135
136 // Find the index at which the unknown is stored
137 const unsigned u_nodal_index = this->u_index_adv_diff();
138
139 // Local shape function
140 Shape psi(n_node);
141
142 // Find values of shape function
143 shape(s, psi);
144
145 // Initialise value of u
146 values[0] = 0.0;
147
148 // Loop over the local nodes and sum
149 for (unsigned l = 0; l < n_node; l++)
150 {
151 values[0] += this->nodal_value(t, l, u_nodal_index) * psi[l];
152 }
153 }
154
155 /*if (t!=0)
156 {
157 std::string error_message =
158 "Time-dependent version of get_interpolated_values() ";
159 error_message += "not implemented for this element \n";
160 throw
161 OomphLibError(
162 error_message,
163 "RefineableAdvectionDiffusionEquations::get_interpolated_values()",
164 OOMPH_EXCEPTION_LOCATION);
165 }
166 else
167 {
168 //Make sure that we call the appropriate steady version
169 //(the entire function might be overloaded lower down)
170 RefineableAdvectionDiffusionEquations<DIM>::
171 get_interpolated_values(s,values);
172 }
173 }*/
174
175
176 /// Further build: Copy source function pointer from father element
178 {
179 RefineableAdvectionDiffusionEquations<DIM>* cast_father_element_pt =
181 this->father_element_pt());
182
183 // Set the values of the pointers from the father
184 this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
185 this->Wind_fct_pt = cast_father_element_pt->wind_fct_pt();
186 this->Pe_pt = cast_father_element_pt->pe_pt();
187 this->PeSt_pt = cast_father_element_pt->pe_st_pt();
188
189 // Set the ALE status
190 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
191 }
192
193 /// Compute the derivatives of the i-th component of
194 /// velocity at point s with respect
195 /// to all data that can affect its value. In addition, return the global
196 /// equation numbers corresponding to the data.
197 /// Overload the non-refineable version to take account of hanging node
198 /// information
200 Vector<double>& du_ddata,
201 Vector<unsigned>& global_eqn_number)
202 {
203 // Find number of nodes
204 unsigned n_node = this->nnode();
205 // Local shape function
206 Shape psi(n_node);
207 // Find values of shape function at the given local coordinate
208 this->shape(s, psi);
209
210 // Find the index at which the velocity component is stored
211 const unsigned u_nodal_index = this->u_index_adv_diff();
212
213 // Storage for hang info pointer
214 HangInfo* hang_info_pt = 0;
215 // Storage for global equation
216 int global_eqn = 0;
217
218 // Find the number of dofs associated with interpolated u
219 unsigned n_u_dof = 0;
220 for (unsigned l = 0; l < n_node; l++)
221 {
222 unsigned n_master = 1;
223
224 // Local bool (is the node hanging)
225 bool is_node_hanging = this->node_pt(l)->is_hanging();
226
227 // If the node is hanging, get the number of master nodes
228 if (is_node_hanging)
229 {
230 hang_info_pt = this->node_pt(l)->hanging_pt();
231 n_master = hang_info_pt->nmaster();
232 }
233 // Otherwise there is just one master node, the node itself
234 else
235 {
236 n_master = 1;
237 }
238
239 // Loop over the master nodes
240 for (unsigned m = 0; m < n_master; m++)
241 {
242 // Get the equation number
243 if (is_node_hanging)
244 {
245 // Get the equation number from the master node
246 global_eqn =
247 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
248 }
249 else
250 {
251 // Global equation number
252 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
253 }
254
255 // If it's positive add to the count
256 if (global_eqn >= 0)
257 {
258 ++n_u_dof;
259 }
260 }
261 }
262
263 // Now resize the storage schemes
264 du_ddata.resize(n_u_dof, 0.0);
265 global_eqn_number.resize(n_u_dof, 0);
266
267 // Loop over th nodes again and set the derivatives
268 unsigned count = 0;
269 // Loop over the local nodes and sum
270 for (unsigned l = 0; l < n_node; l++)
271 {
272 unsigned n_master = 1;
273 double hang_weight = 1.0;
274
275 // Local bool (is the node hanging)
276 bool is_node_hanging = this->node_pt(l)->is_hanging();
277
278 // If the node is hanging, get the number of master nodes
279 if (is_node_hanging)
280 {
281 hang_info_pt = this->node_pt(l)->hanging_pt();
282 n_master = hang_info_pt->nmaster();
283 }
284 // Otherwise there is just one master node, the node itself
285 else
286 {
287 n_master = 1;
288 }
289
290 // Loop over the master nodes
291 for (unsigned m = 0; m < n_master; m++)
292 {
293 // If the node is hanging get weight from master node
294 if (is_node_hanging)
295 {
296 // Get the hang weight from the master node
297 hang_weight = hang_info_pt->master_weight(m);
298 }
299 else
300 {
301 // Node contributes with full weight
302 hang_weight = 1.0;
303 }
304
305 // Get the equation number
306 if (is_node_hanging)
307 {
308 // Get the equation number from the master node
309 global_eqn =
310 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
311 }
312 else
313 {
314 // Global equation number
315 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
316 }
317
318 if (global_eqn >= 0)
319 {
320 // Set the global equation number
321 global_eqn_number[count] = global_eqn;
322 // Set the derivative with respect to the unknown
323 du_ddata[count] = psi[l] * hang_weight;
324 // Increase the counter
325 ++count;
326 }
327 }
328 }
329 }
330
331
332 protected:
333 /// Add the element's contribution to the elemental residual vector
334 /// and/or Jacobian matrix
335 /// flag=1: compute both
336 /// flag=0: compute only residual vector
338 Vector<double>& residuals,
339 DenseMatrix<double>& jacobian,
340 DenseMatrix<double>& mass_matrix,
341 unsigned flag);
342 };
343
344
345 //======================================================================
346 /// Refineable version of QAdvectionDiffusionElement.
347 /// Inherit from the standard QAdvectionDiffusionElement and the
348 /// appropriate refineable geometric element and the refineable equations.
349 //======================================================================
350 template<unsigned DIM, unsigned NNODE_1D>
352 : public QAdvectionDiffusionElement<DIM, NNODE_1D>,
353 public virtual RefineableAdvectionDiffusionEquations<DIM>,
354 public virtual RefineableQElement<DIM>
355 {
356 public:
357 /// Empty Constructor:
361 RefineableQElement<DIM>(),
362 QAdvectionDiffusionElement<DIM, NNODE_1D>()
363 {
364 }
365
366
367 /// Broken copy constructor
370 delete;
371
372 /// Broken assignment operator
374 delete;
375
376 /// Number of continuously interpolated values: 1
378 {
379 return 1;
380 }
381
382 /// Number of vertex nodes in the element
383 unsigned nvertex_node() const
384 {
386 }
387
388 /// Pointer to the j-th vertex node in the element
389 Node* vertex_node_pt(const unsigned& j) const
390 {
392 }
393
394 /// Rebuild from sons: empty
395 void rebuild_from_sons(Mesh*& mesh_pt) {}
396
397 /// Order of recovery shape functions for Z2 error estimation:
398 /// Same order as shape functions.
400 {
401 return (NNODE_1D - 1);
402 }
403
404 /// Perform additional hanging node procedures for variables
405 /// that are not interpolated by all nodes. Empty.
407 };
408
409 /// /////////////////////////////////////////////////////////////////////
410 /// /////////////////////////////////////////////////////////////////////
411 /// /////////////////////////////////////////////////////////////////////
412
413
414 //=======================================================================
415 /// Face geometry for the RefineableQuadAdvectionDiffusionElement elements:
416 /// The spatial dimension of the face elements is one lower than that of the
417 /// bulk element but they have the same number of points
418 /// along their 1D edges.
419 //=======================================================================
420 template<unsigned DIM, unsigned NNODE_1D>
422 : public virtual QElement<DIM - 1, NNODE_1D>
423 {
424 public:
425 /// Constructor: Call the constructor for the
426 /// appropriate lower-dimensional QElement
427 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
428 };
429
430} // namespace oomph
431
432#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 equations using isoparametric elements.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
double * Pe_pt
Pointer to global Peclet number.
double *& pe_pt()
Pointer to Peclet number.
AdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
AdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
AdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
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...
AdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
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
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
A version of the Advection Diffusion equations that can be used with non-uniform mesh refinement....
void further_build()
Further build: Copy source function pointer from father element.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
RefineableAdvectionDiffusionEquations(const RefineableAdvectionDiffusionEquations< DIM > &dummy)=delete
Broken copy constructor.
void operator=(const RefineableAdvectionDiffusionEquations< DIM > &)=delete
Broken assignment operator.
void fill_in_generic_residual_contribution_adv_diff(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 dinterpolated_u_adv_diff_ddata(const Vector< double > &s, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Compute the derivatives of the i-th component of velocity at point s with respect to all data that ca...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Standard flux.from AdvectionDiffusion equations.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
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 QAdvectionDiffusionElement. Inherit from the standard QAdvectionDiffusionElemen...
unsigned nvertex_node() const
Number of vertex nodes in the element.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
RefineableQAdvectionDiffusionElement(const RefineableQAdvectionDiffusionElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
void operator=(const RefineableQAdvectionDiffusionElement< DIM, NNODE_1D > &)=delete
Broken assignment operator.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...