refineable_spherical_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_SPHERICAL_ADVECTION_DIFFUSION_ELEMENTS_HEADER
30#define OOMPH_REFINEABLE_SPHERICAL_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 spherical
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_spherical_adv_diff()
50 /// function so that contributions
51 /// from hanging nodes (or alternatively in-compatible function values)
52 /// are taken into account.
53 //======================================================================
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
75 delete;
76
77 /// Number of 'flux' terms for Z2 error estimation
79 {
80 return 2;
81 }
82
83 /// Get 'flux' for Z2 error recovery:
84 /// Standard flux.from AdvectionDiffusion equations
86 {
87 this->get_flux(s, flux);
88 }
89
90
91 /// Get the function value u in Vector.
92 /// Note: Given the generality of the interface (this function
93 /// is usually called from black-box documentation or interpolation
94 /// routines), the values Vector sets its own size in here.
96 Vector<double>& values)
97 {
98 // Set size of Vector: u
99 values.resize(1);
100
101 // Find number of nodes
102 unsigned n_node = nnode();
103
104 // Find the index at which the unknown is stored
105 unsigned u_nodal_index = this->u_index_spherical_adv_diff();
106
107 // Local shape function
108 Shape psi(n_node);
109
110 // Find values of shape function
111 shape(s, psi);
112
113 // Initialise value of u
114 values[0] = 0.0;
115
116 // Loop over the local nodes and sum
117 for (unsigned l = 0; l < n_node; l++)
118 {
119 values[0] += this->nodal_value(l, u_nodal_index) * psi[l];
120 }
121 }
122
123 /// Get the function value u in Vector.
124 /// Note: Given the generality of the interface (this function
125 /// is usually called from black-box documentation or interpolation
126 /// routines), the values Vector sets its own size in here.
127 void get_interpolated_values(const unsigned& t,
128 const Vector<double>& s,
129 Vector<double>& values)
130 {
131 // Set size of Vector: u
132 values.resize(1);
133
134 // Find number of nodes
135 const unsigned n_node = nnode();
136
137 // Find the index at which the unknown is stored
138 const unsigned u_nodal_index = this->u_index_spherical_adv_diff();
139
140 // Local shape function
141 Shape psi(n_node);
142
143 // Find values of shape function
144 shape(s, psi);
145
146 // Initialise value of u
147 values[0] = 0.0;
148
149 // Loop over the local nodes and sum
150 for (unsigned l = 0; l < n_node; l++)
151 {
152 values[0] += this->nodal_value(t, l, u_nodal_index) * psi[l];
153 }
154 // }
155
156 if (t != 0)
157 {
158 std::string error_message =
159 "Time-dependent version of get_interpolated_values() ";
160 error_message += "not implemented for this element \n";
161 throw OomphLibError(error_message,
162 "RefineableSphericalAdvectionDiffusionEquations::"
163 "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)
171 s, values);
172 }
173 }
174
175 /// Fill in the geometric Jacobian, which in this case is r*r*sin(theta)
177 {
178 return x[0] * x[0] * sin(x[1]);
179 }
180
181 /// Further build: Copy source function pointer from father element
183 {
184 RefineableSphericalAdvectionDiffusionEquations* cast_father_element_pt =
186 this->father_element_pt());
187
188 // Set the values of the pointers from the father
189 this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
190 this->Wind_fct_pt = cast_father_element_pt->wind_fct_pt();
191 this->Pe_pt = cast_father_element_pt->pe_pt();
192 this->PeSt_pt = cast_father_element_pt->pe_st_pt();
193
194 // Set the ALE status
195 // this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
196 }
197
198 /// Compute the derivatives of the i-th component of
199 /// velocity at point s with respect
200 /// to all data that can affect its value. In addition, return the global
201 /// equation numbers corresponding to the data.
202 /// Overload the non-refineable version to take account of hanging node
203 /// information
205 Vector<double>& du_ddata,
206 Vector<unsigned>& global_eqn_number)
207 {
208 // Find number of nodes
209 unsigned n_node = this->nnode();
210 // Local shape function
211 Shape psi(n_node);
212 // Find values of shape function at the given local coordinate
213 this->shape(s, psi);
214
215 // Find the index at which the velocity component is stored
216 const unsigned u_nodal_index = this->u_index_spherical_adv_diff();
217
218 // Storage for hang info pointer
219 HangInfo* hang_info_pt = 0;
220 // Storage for global equation
221 int global_eqn = 0;
222
223 // Find the number of dofs associated with interpolated u
224 unsigned n_u_dof = 0;
225 for (unsigned l = 0; l < n_node; l++)
226 {
227 unsigned n_master = 1;
228
229 // Local bool (is the node hanging)
230 bool is_node_hanging = this->node_pt(l)->is_hanging();
231
232 // If the node is hanging, get the number of master nodes
233 if (is_node_hanging)
234 {
235 hang_info_pt = this->node_pt(l)->hanging_pt();
236 n_master = hang_info_pt->nmaster();
237 }
238 // Otherwise there is just one master node, the node itself
239 else
240 {
241 n_master = 1;
242 }
243
244 // Loop over the master nodes
245 for (unsigned m = 0; m < n_master; m++)
246 {
247 // Get the equation number
248 if (is_node_hanging)
249 {
250 // Get the equation number from the master node
251 global_eqn =
252 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
253 }
254 else
255 {
256 // Global equation number
257 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
258 }
259
260 // If it's positive add to the count
261 if (global_eqn >= 0)
262 {
263 ++n_u_dof;
264 }
265 }
266 }
267
268 // Now resize the storage schemes
269 du_ddata.resize(n_u_dof, 0.0);
270 global_eqn_number.resize(n_u_dof, 0);
271
272 // Loop over th nodes again and set the derivatives
273 unsigned count = 0;
274 // Loop over the local nodes and sum
275 for (unsigned l = 0; l < n_node; l++)
276 {
277 unsigned n_master = 1;
278 double hang_weight = 1.0;
279
280 // Local bool (is the node hanging)
281 bool is_node_hanging = this->node_pt(l)->is_hanging();
282
283 // If the node is hanging, get the number of master nodes
284 if (is_node_hanging)
285 {
286 hang_info_pt = this->node_pt(l)->hanging_pt();
287 n_master = hang_info_pt->nmaster();
288 }
289 // Otherwise there is just one master node, the node itself
290 else
291 {
292 n_master = 1;
293 }
294
295 // Loop over the master nodes
296 for (unsigned m = 0; m < n_master; m++)
297 {
298 // If the node is hanging get weight from master node
299 if (is_node_hanging)
300 {
301 // Get the hang weight from the master node
302 hang_weight = hang_info_pt->master_weight(m);
303 }
304 else
305 {
306 // Node contributes with full weight
307 hang_weight = 1.0;
308 }
309
310 // Get the equation number
311 if (is_node_hanging)
312 {
313 // Get the equation number from the master node
314 global_eqn =
315 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
316 }
317 else
318 {
319 // Global equation number
320 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
321 }
322
323 if (global_eqn >= 0)
324 {
325 // Set the global equation number
326 global_eqn_number[count] = global_eqn;
327 // Set the derivative with respect to the unknown
328 du_ddata[count] = psi[l] * hang_weight;
329 // Increase the counter
330 ++count;
331 }
332 }
333 }
334 }
335
336
337 protected:
338 /// Add the element's contribution to the elemental residual vector
339 /// and/or Jacobian matrix
340 /// flag=1: compute both
341 /// flag=0: compute only residual vector
343 Vector<double>& residuals,
344 DenseMatrix<double>& jacobian,
345 DenseMatrix<double>& mass_matrix,
346 unsigned flag);
347 };
348
349
350 //======================================================================
351 /// Refineable version of QSphericalAdvectionDiffusionElement.
352 /// Inherit from the standard QSphericalAdvectionDiffusionElement and the
353 /// appropriate refineable geometric element and the refineable equations.
354 //======================================================================
355 template<unsigned NNODE_1D>
357 : public QSphericalAdvectionDiffusionElement<NNODE_1D>,
359 public virtual RefineableQElement<2>
360 {
361 public:
362 /// Empty Constructor:
368 {
369 }
370
371
372 /// Broken copy constructor
375 delete;
376
377 /// Broken assignment operator
380
381 /// Number of continuously interpolated values: 1
383 {
384 return 1;
385 }
386
387 /// Number of vertex nodes in the element
388 unsigned nvertex_node() const
389 {
391 }
392
393 /// Pointer to the j-th vertex node in the element
394 Node* vertex_node_pt(const unsigned& j) const
395 {
397 }
398
399 /// Rebuild from sons: empty
400 void rebuild_from_sons(Mesh*& mesh_pt) {}
401
402 /// Order of recovery shape functions for Z2 error estimation:
403 /// Same order as shape functions.
405 {
406 return (NNODE_1D - 1);
407 }
408
409 /// Perform additional hanging node procedures for variables
410 /// that are not interpolated by all nodes. Empty.
412 };
413
414 /// /////////////////////////////////////////////////////////////////////
415 /// /////////////////////////////////////////////////////////////////////
416 /// /////////////////////////////////////////////////////////////////////
417
418
419 //=======================================================================
420 /// Face geometry for the RefineableQSphericalAdvectionDiffusionElement
421 /// elements: The spatial
422 /// dimension of the face elements is one lower than that of the
423 /// bulk element but they have the same number of points
424 /// along their 1D edges.
425 //=======================================================================
426 template<unsigned NNODE_1D>
428 : public virtual QElement<1, NNODE_1D>
429 {
430 public:
431 /// Constructor: Call the constructor for the
432 /// appropriate lower-dimensional QElement
433 FaceGeometry() : QElement<1, NNODE_1D>() {}
434 };
435
436} // namespace oomph
437
438#endif
static char t char * s
Definition: cfortran.h:568
char t
Definition: cfortran.h:568
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
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
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.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
Refineable version of QSphericalAdvectionDiffusionElement. Inherit from the standard QSphericalAdvect...
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.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
RefineableQSphericalAdvectionDiffusionElement(const RefineableQSphericalAdvectionDiffusionElement< NNODE_1D > &dummy)=delete
Broken copy constructor.
void operator=(const RefineableQSphericalAdvectionDiffusionElement< NNODE_1D > &)=delete
Broken assignment operator.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
A version of the Advection Diffusion in spherical coordinates equations that can be used with non-uni...
void fill_in_generic_residual_contribution_spherical_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 further_build()
Further build: Copy source function pointer from father element.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Standard flux.from AdvectionDiffusion equations.
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...
RefineableSphericalAdvectionDiffusionEquations(const RefineableSphericalAdvectionDiffusionEquations &dummy)=delete
Broken copy constructor.
double geometric_jacobian(const Vector< double > &x)
Fill in the geometric Jacobian, which in this case is r*r*sin(theta)
void operator=(const RefineableSphericalAdvectionDiffusionEquations &)=delete
Broken assignment operator.
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...
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 elements that solve the Advection Diffusion equations in a spherical polar coordinate...
SphericalAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
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:
SphericalAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
SphericalAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
SphericalAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
virtual unsigned u_index_spherical_adv_diff() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...