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