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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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_ADVECTION_DIFFUSION_ELEMENTS_HEADER
34 #define OOMPH_REFINEABLE_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 equations that can be
52 /// used with non-uniform mesh refinement. In essence, the class overloads
53 /// the fill_in_generic_residual_contribution_adv_diff()
54 /// function so that contributions
55 /// from hanging nodes (or alternatively in-compatible function values)
56 /// are taken into account.
57 //======================================================================
58 template <unsigned DIM>
60  public virtual AdvectionDiffusionEquations<DIM>,
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("RefineableAdvectionDiffusionEquations");
79  }
80 
81  /// Broken assignment operator
83  {
84  BrokenCopy::broken_assign("RefineableAdvectionDiffusionEquations");
85  }
86 
87  /// Number of 'flux' terms for Z2 error estimation
88  unsigned num_Z2_flux_terms() {return DIM;}
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_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_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
165  OomphLibError(
166  error_message,
167  "RefineableAdvectionDiffusionEquations::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)
174  RefineableAdvectionDiffusionEquations<DIM>::
175  get_interpolated_values(s,values);
176  }
177  }*/
178 
179 
180  /// Further build: Copy source function pointer from father element
182  {
183  RefineableAdvectionDiffusionEquations<DIM>* cast_father_element_pt
185  this->father_element_pt());
186 
187  //Set the values of the pointers from the father
188  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
189  this->Wind_fct_pt = cast_father_element_pt->wind_fct_pt();
190  this->Pe_pt = cast_father_element_pt->pe_pt();
191  this->PeSt_pt = cast_father_element_pt->pe_st_pt();
192 
193  //Set the ALE status
194  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
195  }
196 
197  /// \short Compute the derivatives of the i-th component of
198  /// velocity at point s with respect
199  /// to all data that can affect its value. In addition, return the global
200  /// equation numbers corresponding to the data.
201  /// Overload the non-refineable version to take account of hanging node
202  /// information
204  Vector<double> &du_ddata,
205  Vector<unsigned> &global_eqn_number)
206  {
207  //Find number of nodes
208  unsigned n_node = this->nnode();
209  //Local shape function
210  Shape psi(n_node);
211  //Find values of shape function at the given local coordinate
212  this->shape(s,psi);
213 
214  //Find the index at which the velocity component is stored
215  const unsigned u_nodal_index = this->u_index_adv_diff();
216 
217  //Storage for hang info pointer
218  HangInfo* hang_info_pt=0;
219  //Storage for global equation
220  int global_eqn = 0;
221 
222  //Find the number of dofs associated with interpolated u
223  unsigned n_u_dof=0;
224  for(unsigned l=0;l<n_node;l++)
225  {
226  unsigned n_master = 1;
227 
228  //Local bool (is the node hanging)
229  bool is_node_hanging = this->node_pt(l)->is_hanging();
230 
231  //If the node is hanging, get the number of master nodes
232  if(is_node_hanging)
233  {
234  hang_info_pt = this->node_pt(l)->hanging_pt();
235  n_master = hang_info_pt->nmaster();
236  }
237  //Otherwise there is just one master node, the node itself
238  else
239  {
240  n_master = 1;
241  }
242 
243  //Loop over the master nodes
244  for(unsigned m=0;m<n_master;m++)
245  {
246  //Get the equation number
247  if(is_node_hanging)
248  {
249  //Get the equation number from the master node
250  global_eqn = hang_info_pt->master_node_pt(m)->
251  eqn_number(u_nodal_index);
252  }
253  else
254  {
255  // Global equation number
256  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
257  }
258 
259  //If it's positive add to the count
260  if (global_eqn >= 0) {++n_u_dof;}
261  }
262  }
263 
264  //Now resize the storage schemes
265  du_ddata.resize(n_u_dof,0.0);
266  global_eqn_number.resize(n_u_dof,0);
267 
268  //Loop over th nodes again and set the derivatives
269  unsigned count=0;
270  //Loop over the local nodes and sum
271  for(unsigned l=0;l<n_node;l++)
272  {
273  unsigned n_master = 1;
274  double hang_weight = 1.0;
275 
276  //Local bool (is the node hanging)
277  bool is_node_hanging = this->node_pt(l)->is_hanging();
278 
279  //If the node is hanging, get the number of master nodes
280  if(is_node_hanging)
281  {
282  hang_info_pt = this->node_pt(l)->hanging_pt();
283  n_master = hang_info_pt->nmaster();
284  }
285  //Otherwise there is just one master node, the node itself
286  else
287  {
288  n_master = 1;
289  }
290 
291  //Loop over the master nodes
292  for(unsigned m=0;m<n_master;m++)
293  {
294  //If the node is hanging get weight from master node
295  if(is_node_hanging)
296  {
297  //Get the hang weight from the master node
298  hang_weight = hang_info_pt->master_weight(m);
299  }
300  else
301  {
302  // Node contributes with full weight
303  hang_weight = 1.0;
304  }
305 
306  //Get the equation number
307  if(is_node_hanging)
308  {
309  //Get the equation number from the master node
310  global_eqn = hang_info_pt->master_node_pt(m)->
311  eqn_number(u_nodal_index);
312  }
313  else
314  {
315  // Global equation number
316  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
317  }
318 
319  if (global_eqn >= 0)
320  {
321  //Set the global equation number
322  global_eqn_number[count] = global_eqn;
323  //Set the derivative with respect to the unknown
324  du_ddata[count] = psi[l]*hang_weight;
325  //Increase the counter
326  ++count;
327  }
328  }
329  }
330  }
331 
332 
333  protected:
334 
335 /// \short Add the element's contribution to the elemental residual vector
336 /// and/or Jacobian matrix
337 /// flag=1: compute both
338 /// flag=0: compute only residual vector
340  Vector<double> &residuals, DenseMatrix<double> &jacobian,
341  DenseMatrix<double> &mass_matrix, unsigned flag);
342 
343 };
344 
345 
346 //======================================================================
347 /// \short Refineable version of QAdvectionDiffusionElement.
348 /// Inherit from the standard QAdvectionDiffusionElement and the
349 /// appropriate refineable geometric element and the refineable equations.
350 //======================================================================
351 template <unsigned DIM, unsigned NNODE_1D>
353 public QAdvectionDiffusionElement<DIM,NNODE_1D>,
354  public virtual RefineableAdvectionDiffusionEquations<DIM>,
355  public virtual RefineableQElement<DIM>
356 {
357  public:
358 
359  /// \short Empty Constructor:
363  RefineableQElement<DIM>(),
364  QAdvectionDiffusionElement<DIM,NNODE_1D>()
365  {}
366 
367 
368  /// Broken copy constructor
371  dummy)
372  {
373  BrokenCopy::broken_copy("RefineableQuadAdvectionDiffusionElement");
374  }
375 
376  /// Broken assignment operator
378  {
379  BrokenCopy::broken_assign("RefineableQuadAdvectionDiffusionElement");
380  }
381 
382  /// Number of continuously interpolated values: 1
383  unsigned ncont_interpolated_values() const {return 1;}
384 
385  /// \short Number of vertex nodes in the element
386  unsigned nvertex_node() const
388 
389  /// \short Pointer to the j-th vertex node in the element
390  Node* vertex_node_pt(const unsigned& j) const
392 
393  /// Rebuild from sons: empty
394  void rebuild_from_sons(Mesh* &mesh_pt) {}
395 
396  /// \short Order of recovery shape functions for Z2 error estimation:
397  /// Same order as shape functions.
398  unsigned nrecovery_order() {return (NNODE_1D-1);}
399 
400  /// \short Perform additional hanging node procedures for variables
401  /// that are not interpolated by all nodes. Empty.
403 
404 };
405 
406 ////////////////////////////////////////////////////////////////////////
407 ////////////////////////////////////////////////////////////////////////
408 ////////////////////////////////////////////////////////////////////////
409 
410 
411 
412 //=======================================================================
413 /// Face geometry for the RefineableQuadAdvectionDiffusionElement elements: The spatial
414 /// dimension of the face elements is one lower than that of the
415 /// bulk element but they have the same number of points
416 /// along their 1D edges.
417 //=======================================================================
418 template<unsigned DIM, unsigned NNODE_1D>
420  public virtual QElement<DIM-1,NNODE_1D>
421 {
422 
423  public:
424 
425  /// \short Constructor: Call the constructor for the
426  /// appropriate lower-dimensional QElement
427  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
428 
429 };
430 
431 }
432 
433 #endif
434 
unsigned nvertex_node() const
Number of vertex nodes in the element.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
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...
AdvectionDiffusionSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
AdvectionDiffusionSourceFctPt Source_fct_pt
Pointer to source function:
void further_build()
Further build: Copy source function pointer from father element.
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
A version of the Advection Diffusion equations that can be used with non-uniform mesh refinement...
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:2458
QAdvectionDiffusionElement elements are linear/quadrilateral/brick-shaped Advection Diffusion element...
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...
char t
Definition: cfortran.h:572
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
double * Pe_pt
Pointer to global Peclet number.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the 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...
A class for all elements that solve the Advection Diffusion equations using isoparametric elements...
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
Pointer to wind function:
double *& pe_pt()
Pointer to Peclet number.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Standard flux.from AdvectionDiffusion equations.
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...
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
RefineableAdvectionDiffusionEquations(const RefineableAdvectionDiffusionEquations< DIM > &dummy)
Broken copy constructor.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
static char t char * s
Definition: cfortran.h:572
AdvectionDiffusionWindFctPt & wind_fct_pt()
Access function: Pointer to wind function.
void operator=(const RefineableQAdvectionDiffusionElement< DIM, NNODE_1D > &)
Broken assignment operator.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
Class that contains data for hanging nodes.
Definition: nodes.h:684
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
void operator=(const RefineableAdvectionDiffusionEquations< DIM > &)
Broken assignment operator.
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:2371
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 broken_assign(const std::string &class_name)
Issue error message and terminate execution.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
double *& pe_st_pt()
Pointer to Peclet number multipled by Strouha number.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
virtual unsigned nvertex_node() const
Definition: elements.h:2362
double * PeSt_pt
Pointer to global Peclet number multiplied by Strouhal number.
Refineable version of QAdvectionDiffusionElement. Inherit from the standard QAdvectionDiffusionElemen...
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...
A general mesh class.
Definition: mesh.h:74
RefineableQAdvectionDiffusionElement(const RefineableQAdvectionDiffusionElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.