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// Version 1.0; svn revision $LastChangedRevision$
7 //LIC//
8 //LIC// $LastChangedDate$
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for elements that solve the advection diffusion equation
31 //and that can be refined.
32 
33 #ifndef OOMPH_REFINEABLE_SPHERICAL_ADVECTION_DIFFUSION_ELEMENTS_HEADER
34 #define OOMPH_REFINEABLE_SPHERICAL_ADVECTION_DIFFUSION_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 //oomph-lib headers
42 #include "../generic/refineable_quad_element.h"
43 #include "../generic/refineable_brick_element.h"
44 #include "../generic/error_estimator.h"
46 
47 namespace oomph
48 {
49 
50 //======================================================================
51 /// \short A version of the Advection Diffusion in spherical
52 /// coordinates equations that can be
53 /// used with non-uniform mesh refinement. In essence, the class overloads
54 /// the fill_in_generic_residual_contribution_spherical_adv_diff()
55 /// function so that contributions
56 /// from hanging nodes (or alternatively in-compatible function values)
57 /// are taken into account.
58 //======================================================================
61  public virtual RefineableElement,
62  public virtual ElementWithZ2ErrorEstimator
63 {
64  public:
65 
66  /// \short Empty Constructor
71  {}
72 
73 
74  /// Broken copy constructor
77  {
78  BrokenCopy::broken_copy("RefineableSphericalAdvectionDiffusionEquations");
79  }
80 
81  /// Broken assignment operator
83  {
84  BrokenCopy::broken_assign("RefineableSphericalAdvectionDiffusionEquations");
85  }
86 
87  /// Number of 'flux' terms for Z2 error estimation
88  unsigned num_Z2_flux_terms() {return 2;}
89 
90  /// \short Get 'flux' for Z2 error recovery:
91  /// Standard flux.from AdvectionDiffusion equations
93  {this->get_flux(s,flux);}
94 
95 
96  /// \short 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 routines),
99  /// the values Vector sets its own size in here.
101  {
102  // Set size of Vector: u
103  values.resize(1);
104 
105  //Find number of nodes
106  unsigned n_node = nnode();
107 
108  //Find the index at which the unknown is stored
109  unsigned u_nodal_index = this->u_index_spherical_adv_diff();
110 
111  //Local shape function
112  Shape psi(n_node);
113 
114  //Find values of shape function
115  shape(s,psi);
116 
117  //Initialise value of u
118  values[0] = 0.0;
119 
120  //Loop over the local nodes and sum
121  for(unsigned l=0;l<n_node;l++)
122  {
123  values[0] += this->nodal_value(l,u_nodal_index)*psi[l];
124  }
125  }
126 
127  /// \short Get the function value u in Vector.
128  /// Note: Given the generality of the interface (this function
129  /// is usually called from black-box documentation or interpolation routines),
130  /// the values Vector sets its own size in here.
131  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
132  Vector<double>& values)
133  {
134  // Set size of Vector: u
135  values.resize(1);
136 
137  //Find number of nodes
138  const unsigned n_node = nnode();
139 
140  //Find the index at which the unknown is stored
141  const unsigned u_nodal_index = this->u_index_spherical_adv_diff();
142 
143  //Local shape function
144  Shape psi(n_node);
145 
146  //Find values of shape function
147  shape(s,psi);
148 
149  //Initialise value of u
150  values[0] = 0.0;
151 
152  //Loop over the local nodes and sum
153  for(unsigned l=0;l<n_node;l++)
154  {
155  values[0] += this->nodal_value(t,l,u_nodal_index)*psi[l];
156  }
157  // }
158 
159  if (t!=0)
160  {
161  std::string error_message =
162  "Time-dependent version of get_interpolated_values() ";
163  error_message += "not implemented for this element \n";
164  throw
166  error_message,
167  "RefineableSphericalAdvectionDiffusionEquations::get_interpolated_values()",
168  OOMPH_EXCEPTION_LOCATION);
169  }
170  else
171  {
172  //Make sure that we call the appropriate steady version
173  //(the entire function might be overloaded lower down)
175  get_interpolated_values(s,values);
176  }
177  }
178 
179  /// Fill in the geometric Jacobian, which in this case is r*r*sin(theta)
181  {return x[0]*x[0]*sin(x[1]);}
182 
183  /// Further build: Copy source function pointer from father element
185  {
186  RefineableSphericalAdvectionDiffusionEquations* cast_father_element_pt
188  this->father_element_pt());
189 
190  //Set the values of the pointers from the father
191  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
192  this->Wind_fct_pt = cast_father_element_pt->wind_fct_pt();
193  this->Pe_pt = cast_father_element_pt->pe_pt();
194  this->PeSt_pt = cast_father_element_pt->pe_st_pt();
195 
196  //Set the ALE status
197  //this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
198  }
199 
200  /// \short Compute the derivatives of the i-th component of
201  /// velocity at point s with respect
202  /// to all data that can affect its value. In addition, return the global
203  /// equation numbers corresponding to the data.
204  /// Overload the non-refineable version to take account of hanging node
205  /// information
207  Vector<double> &du_ddata,
208  Vector<unsigned> &global_eqn_number)
209  {
210  //Find number of nodes
211  unsigned n_node = this->nnode();
212  //Local shape function
213  Shape psi(n_node);
214  //Find values of shape function at the given local coordinate
215  this->shape(s,psi);
216 
217  //Find the index at which the velocity component is stored
218  const unsigned u_nodal_index = this->u_index_spherical_adv_diff();
219 
220  //Storage for hang info pointer
221  HangInfo* hang_info_pt=0;
222  //Storage for global equation
223  int global_eqn = 0;
224 
225  //Find the number of dofs associated with interpolated u
226  unsigned n_u_dof=0;
227  for(unsigned l=0;l<n_node;l++)
228  {
229  unsigned n_master = 1;
230 
231  //Local bool (is the node hanging)
232  bool is_node_hanging = this->node_pt(l)->is_hanging();
233 
234  //If the node is hanging, get the number of master nodes
235  if(is_node_hanging)
236  {
237  hang_info_pt = this->node_pt(l)->hanging_pt();
238  n_master = hang_info_pt->nmaster();
239  }
240  //Otherwise there is just one master node, the node itself
241  else
242  {
243  n_master = 1;
244  }
245 
246  //Loop over the master nodes
247  for(unsigned m=0;m<n_master;m++)
248  {
249  //Get the equation number
250  if(is_node_hanging)
251  {
252  //Get the equation number from the master node
253  global_eqn = hang_info_pt->master_node_pt(m)->
254  eqn_number(u_nodal_index);
255  }
256  else
257  {
258  // Global equation number
259  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
260  }
261 
262  //If it's positive add to the count
263  if (global_eqn >= 0) {++n_u_dof;}
264  }
265  }
266 
267  //Now resize the storage schemes
268  du_ddata.resize(n_u_dof,0.0);
269  global_eqn_number.resize(n_u_dof,0);
270 
271  //Loop over th nodes again and set the derivatives
272  unsigned count=0;
273  //Loop over the local nodes and sum
274  for(unsigned l=0;l<n_node;l++)
275  {
276  unsigned n_master = 1;
277  double hang_weight = 1.0;
278 
279  //Local bool (is the node hanging)
280  bool is_node_hanging = this->node_pt(l)->is_hanging();
281 
282  //If the node is hanging, get the number of master nodes
283  if(is_node_hanging)
284  {
285  hang_info_pt = this->node_pt(l)->hanging_pt();
286  n_master = hang_info_pt->nmaster();
287  }
288  //Otherwise there is just one master node, the node itself
289  else
290  {
291  n_master = 1;
292  }
293 
294  //Loop over the master nodes
295  for(unsigned m=0;m<n_master;m++)
296  {
297  //If the node is hanging get weight from master node
298  if(is_node_hanging)
299  {
300  //Get the hang weight from the master node
301  hang_weight = hang_info_pt->master_weight(m);
302  }
303  else
304  {
305  // Node contributes with full weight
306  hang_weight = 1.0;
307  }
308 
309  //Get the equation number
310  if(is_node_hanging)
311  {
312  //Get the equation number from the master node
313  global_eqn = hang_info_pt->master_node_pt(m)->
314  eqn_number(u_nodal_index);
315  }
316  else
317  {
318  // Global equation number
319  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
320  }
321 
322  if (global_eqn >= 0)
323  {
324  //Set the global equation number
325  global_eqn_number[count] = global_eqn;
326  //Set the derivative with respect to the unknown
327  du_ddata[count] = psi[l]*hang_weight;
328  //Increase the counter
329  ++count;
330  }
331  }
332  }
333  }
334 
335 
336  protected:
337 
338 /// \short 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, DenseMatrix<double> &jacobian,
344  DenseMatrix<double> &mass_matrix, unsigned flag);
345 
346 };
347 
348 
349 //======================================================================
350 /// \short Refineable version of QSphericalAdvectionDiffusionElement.
351 /// Inherit from the standard QSphericalAdvectionDiffusionElement and the
352 /// appropriate refineable geometric element and the refineable equations.
353 //======================================================================
354 template <unsigned NNODE_1D>
356 public QSphericalAdvectionDiffusionElement<NNODE_1D>,
358  public virtual RefineableQElement<2>
359 {
360  public:
361 
362  /// \short Empty Constructor:
366  RefineableQElement<2>(),
368  {}
369 
370 
371  /// Broken copy constructor
374  dummy)
375  {
376  BrokenCopy::broken_copy("RefineableQSphericalAdvectionDiffusionElement");
377  }
378 
379  /// Broken assignment operator
380  void operator=(const
382  {
383  BrokenCopy::broken_assign("RefineableQSphericalAdvectionDiffusionElement");
384  }
385 
386  /// Number of continuously interpolated values: 1
387  unsigned ncont_interpolated_values() const {return 1;}
388 
389  /// \short Number of vertex nodes in the element
390  unsigned nvertex_node() const
392 
393  /// \short Pointer to the j-th vertex node in the element
394  Node* vertex_node_pt(const unsigned& j) const
396 
397  /// Rebuild from sons: empty
398  void rebuild_from_sons(Mesh* &mesh_pt) {}
399 
400  /// \short Order of recovery shape functions for Z2 error estimation:
401  /// Same order as shape functions.
402  unsigned nrecovery_order() {return (NNODE_1D-1);}
403 
404  /// \short 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 
416 //=======================================================================
417 /// Face geometry for the RefineableQSphericalAdvectionDiffusionElement
418 /// elements: The spatial
419 /// dimension of the face elements is one lower than that of the
420 /// bulk element but they have the same number of points
421 /// along their 1D edges.
422 //=======================================================================
423 template<unsigned NNODE_1D>
425  public virtual QElement<1,NNODE_1D>
426 {
427 
428  public:
429 
430  /// \short Constructor: Call the constructor for the
431  /// appropriate lower-dimensional QElement
432  FaceGeometry() : QElement<1,NNODE_1D>() {}
433 
434 };
435 
436 }
437 
438 #endif
439 
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Definition: Qelements.h:910
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
unsigned num_Z2_flux_terms()
Number of &#39;flux&#39; terms for Z2 error estimation.
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
SphericalAdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
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 get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
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...
RefineableQSphericalAdvectionDiffusionElement(const RefineableQSphericalAdvectionDiffusionElement< NNODE_1D > &dummy)
Broken copy constructor.
char t
Definition: cfortran.h:572
unsigned nvertex_node() const
Number of vertex nodes in the element.
Definition: Qelements.h:906
double geometric_jacobian(const Vector< double > &x)
Fill in the geometric Jacobian, which in this case is r*r*sin(theta)
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
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...
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
void operator=(const RefineableQSphericalAdvectionDiffusionElement< NNODE_1D > &)
Broken assignment operator.
A class for all elements that solve the Advection Diffusion equations in a spherical polar coordinate...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get &#39;flux&#39; for Z2 error recovery: Standard flux.from AdvectionDiffusion equations.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
static char t char * s
Definition: cfortran.h:572
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
void fill_in_generic_residual_contribution_spherical_adv_diff(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add the element&#39;s contribution to the elemental residual vector and/or Jacobian matrix flag=1: comput...
Refineable version of QSphericalAdvectionDiffusionElement. Inherit from the standard QSphericalAdvect...
SphericalAdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2109
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
A version of the Advection Diffusion in spherical coordinates equations that can be used with non-uni...
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:2470
QSphericalAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Axisymmetric Advec...
Class that contains data for hanging nodes.
Definition: nodes.h:684
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn&#39;t been defined.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
SphericalAdvectionDiffusionWindFctPt Wind_fct_pt
Pointer to wind function:
void operator=(const RefineableSphericalAdvectionDiffusionEquations &)
Broken assignment operator.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2146
void further_build()
Further build: Copy source function pointer from father element.
SphericalAdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
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...
RefineableSphericalAdvectionDiffusionEquations(const RefineableSphericalAdvectionDiffusionEquations &dummy)
Broken copy constructor.
A general mesh class.
Definition: mesh.h:74
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...