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