refineable_gen_axisym_advection_diffusion_elements.cc
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//====================================================================
31 
32 namespace oomph
33 {
34 
35 //==========================================================================
36 /// Add the element's contribution to the elemental residual vector
37 /// and/or elemental jacobian matrix.
38 /// This function overloads the standard version so that the possible
39 /// presence of hanging nodes is taken into account.
40 //=========================================================================
43  Vector<double> &residuals,
44  DenseMatrix<double> &jacobian,
46  &mass_matrix,
47  unsigned flag)
48 {
49 
50 //Find out how many nodes there are in the element
51 unsigned n_node = nnode();
52 
53 //Get the nodal index at which the unknown is stored
54  unsigned u_nodal_index = this->u_index_cons_axisym_adv_diff();
55 
56 //Set up memory for the shape and test functions
57  Shape psi(n_node), test(n_node);
58  DShape dpsidx(n_node,2), dtestdx(n_node,2);
59 
60 //Set the value of n_intpt
61 unsigned n_intpt = integral_pt()->nweight();
62 
63 //Set the Vector to hold local coordinates
64 Vector<double> s(2);
65 
66 //Get Peclet number
67 double peclet = this->pe();
68 
69 //Get the Peclet multiplied by the Strouhal number
70 double peclet_st = this->pe_st();
71 
72 //Integers used to store the local equation number and local unknown
73 //indices for the residuals and jacobians
74 int local_eqn=0, local_unknown=0;
75 
76 // Local storage for pointers to hang_info objects
77  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
78 
79 //Local variable to determine the ALE stuff
80  bool ALE_is_disabled_flag = this->ALE_is_disabled;
81 
82 //Loop over the integration points
83 for(unsigned ipt=0;ipt<n_intpt;ipt++)
84 {
85 
86  //Assign values of s
87  for(unsigned i=0;i<2;i++) s[i] = integral_pt()->knot(ipt,i);
88 
89  //Get the integral weight
90  double w = integral_pt()->weight(ipt);
91 
92  //Call the derivatives of the shape and test functions
93  double J =
95  test,dtestdx);
96 
97  //Premultiply the weights and the Jacobian
98  double W = w*J;
99 
100  //Calculate local values of the function, initialise to zero
101  double dudt=0.0;
102  double interpolated_u=0.0;
103 
104  //These need to be a Vector to be ANSI C++, initialise to zero
106  Vector<double> interpolated_dudx(2,0.0);
107  Vector<double> mesh_velocity(2,0.0);
108 
109  //Calculate function value and derivatives:
110  //-----------------------------------------
111 
112  // Loop over nodes
113  for(unsigned l=0;l<n_node;l++)
114  {
115  //Get the value at the node
116  double u_value = this->nodal_value(l,u_nodal_index);
117  interpolated_u += u_value*psi(l);
118  dudt += this->du_dt_cons_axisym_adv_diff(l)*psi(l);
119  // Loop over directions
120  for(unsigned j=0;j<2;j++)
121  {
122  interpolated_x[j] += nodal_position(l,j)*psi(l);
123  interpolated_dudx[j] += u_value*dpsidx(l,j);
124  }
125  }
126 
127  //Get the mesh velocity, if required
128  if (!ALE_is_disabled_flag)
129  {
130  for(unsigned l=0;l<n_node;l++)
131  {
132  // Loop over directions
133  for(unsigned j=0;j<2;j++)
134  {
135  mesh_velocity[j] += dnodal_position_dt(l,j)*psi(l);
136  }
137  }
138  }
139 
140 
141  //Get body force
142  double source;
143  this->get_source_cons_axisym_adv_diff(ipt,interpolated_x,source);
144 
145 
146  //Get wind
147  //--------
148  Vector<double> wind(3);
149  this->get_wind_cons_axisym_adv_diff(ipt,s,interpolated_x,wind);
150 
151  //Get the conserved wind (non-divergence free)
152  Vector<double> conserved_wind(3);
153  this->get_conserved_wind_cons_axisym_adv_diff(ipt,s,interpolated_x,
154  conserved_wind);
155 
156  //Get diffusivity tensor
157  DenseMatrix<double> D(3,3);
158  this->get_diff_cons_axisym_adv_diff(ipt,s,interpolated_x,D);
159 
160  //r is the first component of position
161  double r = interpolated_x[0];
162 
163  // Assemble residuals and Jacobian
164  //================================
165 
166  // Loop over the nodes for the test functions
167  for(unsigned l=0;l<n_node;l++)
168  {
169  //Local variables to store the number of master nodes and
170  //the weight associated with the shape function if the node is hanging
171  unsigned n_master=1; double hang_weight=1.0;
172  //Local bool (is the node hanging)
173  bool is_node_hanging = this->node_pt(l)->is_hanging();
174 
175  //If the node is hanging, get the number of master nodes
176  if(is_node_hanging)
177  {
178  hang_info_pt = this->node_pt(l)->hanging_pt();
179  n_master = hang_info_pt->nmaster();
180  }
181  //Otherwise there is just one master node, the node itself
182  else
183  {
184  n_master = 1;
185  }
186 
187  //Loop over the number of master nodes
188  for(unsigned m=0;m<n_master;m++)
189  {
190  //Get the local equation number and hang_weight
191  //If the node is hanging
192  if(is_node_hanging)
193  {
194  //Read out the local equation from the master node
195  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
196  u_nodal_index);
197  //Read out the weight from the master node
198  hang_weight = hang_info_pt->master_weight(m);
199  }
200  //If the node is not hanging
201  else
202  {
203  //The local equation number comes from the node itself
204  local_eqn = this->nodal_local_eqn(l,u_nodal_index);
205  //The hang weight is one
206  hang_weight = 1.0;
207  }
208 
209  //If the nodal equation is not a boundary conditino
210  if(local_eqn >= 0)
211  {
212  //Add du/dt and body force/source term here
213  residuals[local_eqn] -=
214  (peclet_st*dudt + source)*test(l)*r*W*hang_weight;
215 
216  //The Mesh velocity and GeneralisedAxisymAdvection--Diffusion bit
217  for(unsigned k=0;k<2;k++)
218  {
219  //Terms that multiply the test function
220  double tmp = peclet*wind[k];
221  //If the mesh is moving need to subtract the mesh velocity
222  if(!ALE_is_disabled_flag)
223  {tmp -= peclet_st*mesh_velocity[k];}
224  tmp *= interpolated_dudx[k];
225 
226  //Terms that multiply the derivative of the test function
227  //Advective term
228  double tmp2 = - conserved_wind[k]*interpolated_u;
229  //Now the diffusive term
230  for(unsigned j=0;j<2;j++)
231  {
232  tmp2 += D(k,j)*interpolated_dudx[j];
233  }
234  //Now construct the contribution to the residuals
235  residuals[local_eqn] -= (tmp*test(l) + tmp2*dtestdx(l,k))*r*W*
236  hang_weight;
237  }
238 
239  // Calculate the Jacobian
240  if(flag)
241  {
242  //Local variables to store the number of master nodes
243  //and the weights associated with each hanging node
244  unsigned n_master2=1; double hang_weight2=1.0;
245  //Loop over the nodes for the variables
246  for(unsigned l2=0;l2<n_node;l2++)
247  {
248  //Local bool (is the node hanging)
249  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
250  //If the node is hanging, get the number of master nodes
251  if(is_node2_hanging)
252  {
253  hang_info2_pt = this->node_pt(l2)->hanging_pt();
254  n_master2 = hang_info2_pt->nmaster();
255  }
256  //Otherwise there is one master node, the node itself
257  else
258  {
259  n_master2 = 1;
260  }
261 
262  //Loop over the master nodes
263  for(unsigned m2=0;m2<n_master2;m2++)
264  {
265  //Get the local unknown and weight
266  //If the node is hanging
267  if(is_node2_hanging)
268  {
269  //Read out the local unknown from the master node
270  local_unknown =
271  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
272  u_nodal_index);
273  //Read out the hanging weight from the master node
274  hang_weight2 = hang_info2_pt->master_weight(m2);
275  }
276  //If the node is not hanging
277  else
278  {
279  //The local unknown number comes from the node
280  local_unknown = this->nodal_local_eqn(l2,u_nodal_index);
281  //The hang weight is one
282  hang_weight2 = 1.0;
283  }
284 
285  //If the unknown is not pinned
286  if(local_unknown >= 0)
287  {
288  //Add contribution to Elemental Matrix
289  // Mass matrix du/dt term
290  jacobian(local_eqn,local_unknown)
291  -= peclet_st*test(l)*psi(l2)*
292  this->node_pt(l2)->time_stepper_pt()->weight(1,0)
293  *r*W*hang_weight*hang_weight2;
294 
295  //Add contribution to mass matrix
296  if(flag==2)
297  {
298  mass_matrix(local_eqn,local_unknown) +=
299  peclet_st*test(l)*psi(l2)*r*W*hang_weight*hang_weight2;
300  }
301 
302  //Add contribution to Elemental Matrix
303  for(unsigned k=0;k<2;k++)
304  {
305  //Temporary term used in assembly
306  double tmp = peclet*wind[k];
307  if(!ALE_is_disabled_flag)
308  {tmp -= peclet_st*mesh_velocity[k];}
309  tmp *= dpsidx(l2,k);
310 
311  double tmp2 = - conserved_wind[k]*psi(l2);
312  //Now the diffusive term
313  for(unsigned j=0;j<2;j++)
314  {
315  tmp2 += D(k,j)*dpsidx(l2,j);
316  }
317 
318  //Now assemble Jacobian term
319  jacobian(local_eqn,local_unknown)
320  -= (tmp*test(l) + tmp2*dtestdx(l,k))*W*r*
321  hang_weight*hang_weight2;
322  }
323  }
324  } //End of loop over master nodes
325  } //End of loop over nodes
326  } //End of Jacobian calculation
327 
328  } //End of non-zero equation
329 
330  } //End of loop over the master nodes for residual
331  } //End of loop over nodes
332 
333  } // End of loop over integration points
334 }
335 
336 
337 
338 //====================================================================
339 // Force build of templates
340 //====================================================================
344 
345 }
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
virtual void get_wind_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get wind at (Eulerian) position x and/or local coordinate s. This function is virtual to allow overlo...
Refineable version of QGeneralisedAxisymAdvectionDiffusionElement. Inherit from the standard QGeneral...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
cstr elem_len * i
Definition: cfortran.h:607
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
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2226
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
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
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
const double & pe() const
Peclet number.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
virtual double dshape_and_dtest_eulerian_at_knot_cons_axisym_adv_diff(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
const double & pe_st() const
Peclet number multiplied by Strouhal number.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
virtual unsigned u_index_cons_axisym_adv_diff() const
A class for all elements that solve the Advection Diffusion equations in conservative form using isop...
static char t char * s
Definition: cfortran.h:572
virtual void get_source_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &x, double &source) const
Get source term at (Eulerian) position x. This function is virtual to allow overloading in multi-phys...
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void fill_in_generic_residual_contribution_cons_axisym_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...
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
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
Class that contains data for hanging nodes.
Definition: nodes.h:684
virtual void get_conserved_wind_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &wind) const
Get additional (conservative) wind at (Eulerian) position x and/or local coordinate s...
double du_dt_cons_axisym_adv_diff(const unsigned &n) const
du/dt at local node n.
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
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 void get_diff_cons_axisym_adv_diff(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &D) const
Get diffusivity tensor at (Eulerian) position x and/or local coordinate s. This function is virtual t...