refineable_advection_diffusion_reaction_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 NREAGENT,unsigned DIM>
45  Vector<double> &residuals,
46  DenseMatrix<double> &jacobian,
48  &mass_matrix, unsigned flag)
49 {
50 
51  //Find out how many nodes there are in the element
52  const unsigned n_node = nnode();
53 
54  //Get the nodal index at which the unknown is stored
55  unsigned c_nodal_index[NREAGENT];
56  for(unsigned r=0;r<NREAGENT;r++)
57  {c_nodal_index[r]= this->c_index_adv_diff_react(r);}
58 
59  //Set up memory for the shape and test functions
60  Shape psi(n_node), test(n_node);
61  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
62 
63  //Set the value of n_intpt
64  const unsigned n_intpt = integral_pt()->nweight();
65 
66  //Set the Vector to hold local coordinates
67  Vector<double> s(DIM);
68 
69  //Get diffusion coefficients
70  Vector<double> D = this->diff();
71 
72  //Get the timescales
73  Vector<double> T = this->tau();
74 
75  //Integers used to store the local equation number and local unknown
76  //indices for the residuals and jacobians
77  int local_eqn=0, local_unknown=0;
78 
79  // Local storage for pointers to hang_info objects
80  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
81 
82  //Local variable to determine the ALE stuff
83  bool ALE_is_disabled_flag = this->ALE_is_disabled;
84 
85  //Loop over the integration points
86  for(unsigned ipt=0;ipt<n_intpt;ipt++)
87  {
88 
89  //Assign values of s
90  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
91 
92  //Get the integral weight
93  double w = integral_pt()->weight(ipt);
94 
95  //Call the derivatives of the shape and test functions
96  double J =
97  this->dshape_and_dtest_eulerian_at_knot_adv_diff_react(ipt,psi,dpsidx,
98  test,dtestdx);
99 
100  //Premultiply the weights and the Jacobian
101  double W = w*J;
102 
103  //Calculate local values of the solution and its derivatives
104  //Allocate
105  Vector<double> interpolated_c(NREAGENT,0.0);
106  Vector<double> dcdt(NREAGENT,0.0);
107  Vector<double> interpolated_x(DIM,0.0);
108  DenseMatrix<double> interpolated_dcdx(NREAGENT,DIM,0.0);
109  Vector<double> mesh_velocity(DIM,0.0);
110 
111 
112  //Calculate function value and derivatives:
113  // Loop over nodes
114  for(unsigned l=0;l<n_node;l++)
115  {
116  // Loop over directions to calculate the position
117  for(unsigned j=0;j<DIM;j++)
118  {
119  interpolated_x[j] += nodal_position(l,j)*psi(l);
120  }
121 
122  //Loop over the unknown reagents
123  for(unsigned r=0;r<NREAGENT;r++)
124  {
125  //Get the value at the node
126  const double c_value = nodal_value(l,c_nodal_index[r]);
127 
128  //Calculate the interpolated value
129  interpolated_c[r] += c_value*psi(l);
130  dcdt[r] += this->dc_dt_adv_diff_react(l,r)*psi(l);
131 
132  // Loop over directions to calculate the derivatie
133  for(unsigned j=0;j<DIM;j++)
134  {
135  interpolated_dcdx(r,j) += c_value*dpsidx(l,j);
136  }
137  }
138  }
139 
140  // Mesh velocity?
141  if (!ALE_is_disabled_flag)
142  {
143  for(unsigned l=0;l<n_node;l++)
144  {
145  for(unsigned j=0;j<DIM;j++)
146  {
147  mesh_velocity[j] += dnodal_position_dt(l,j)*psi(l);
148  }
149  }
150  }
151 
152  //Get source function
153  Vector<double> source(NREAGENT);
154  this->get_source_adv_diff_react(ipt,interpolated_x,source);
155 
156 
157  //Get wind
158  Vector<double> wind(DIM);
159  this->get_wind_adv_diff_react(ipt,s,interpolated_x,wind);
160 
161  //Get reaction terms
162  Vector<double> R(NREAGENT);
163  this->get_reaction_adv_diff_react(ipt,interpolated_c,R);
164 
165  //If we are getting the jacobian the get the derivative terms
166  DenseMatrix<double> dRdC(NREAGENT);
167  if(flag)
168  {
169  this->get_reaction_deriv_adv_diff_react(ipt,interpolated_c,dRdC);
170  }
171 
172  // Assemble residuals and Jacobian
173  //================================
174 
175  // Loop over the nodes for the test functions
176  for(unsigned l=0;l<n_node;l++)
177  {
178  //Local variables to store the number of master nodes and
179  //the weight associated with the shape function if the node is hanging
180  unsigned n_master=1; double hang_weight=1.0;
181  //Local bool (is the node hanging)
182  bool is_node_hanging = this->node_pt(l)->is_hanging();
183 
184  //If the node is hanging, get the number of master nodes
185  if(is_node_hanging)
186  {
187  hang_info_pt = this->node_pt(l)->hanging_pt();
188  n_master = hang_info_pt->nmaster();
189  }
190  //Otherwise there is just one master node, the node itself
191  else
192  {
193  n_master = 1;
194  }
195 
196  //Loop over the number of master nodes
197  for(unsigned m=0;m<n_master;m++)
198  {
199  //Loop over the number of reagents
200  for(unsigned r=0;r<NREAGENT;r++)
201  {
202  //Get the local equation number and hang_weight
203  //If the node is hanging
204  if(is_node_hanging)
205  {
206  //Read out the local equation from the master node
207  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
208  c_nodal_index[r]);
209  //Read out the weight from the master node
210  hang_weight = hang_info_pt->master_weight(m);
211  }
212  //If the node is not hanging
213  else
214  {
215  //The local equation number comes from the node itself
216  local_eqn = this->nodal_local_eqn(l,c_nodal_index[r]);
217  //The hang weight is one
218  hang_weight = 1.0;
219  }
220 
221  //If the nodal equation is not a boundary conditino
222  if(local_eqn >= 0)
223  {
224  // Add body force/source/reaction term and time derivative
225  residuals[local_eqn] -=
226  (T[r]*dcdt[r] + source[r] + R[r])*test(l)*W*hang_weight;
227 
228  // The Advection Diffusion bit itself
229  for(unsigned k=0;k<DIM;k++)
230  {
231  //Terms that multiply the test function
232  double tmp = wind[k];
233  //If the mesh is moving need to subtract the mesh velocity
234  if(!ALE_is_disabled_flag) {tmp -= T[r]*mesh_velocity[k];}
235  //Now construct the contribution to the residuals
236  residuals[local_eqn] -=
237  interpolated_dcdx(r,k)*(tmp*test(l) +
238  D[r]*dtestdx(l,k))*W*hang_weight;
239  }
240 
241  // Calculate the Jacobian
242  if(flag)
243  {
244  //Local variables to store the number of master nodes
245  //and the weights associated with each hanging node
246  unsigned n_master2=1; double hang_weight2=1.0;
247  //Loop over the nodes for the variables
248  for(unsigned l2=0;l2<n_node;l2++)
249  {
250  //Local bool (is the node hanging)
251  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
252  //If the node is hanging, get the number of master nodes
253  if(is_node2_hanging)
254  {
255  hang_info2_pt = this->node_pt(l2)->hanging_pt();
256  n_master2 = hang_info2_pt->nmaster();
257  }
258  //Otherwise there is one master node, the node itself
259  else
260  {
261  n_master2 = 1;
262  }
263 
264  //Loop over the master nodes
265  for(unsigned m2=0;m2<n_master2;m2++)
266  {
267  //Loop over the reagents again
268  for(unsigned r2=0;r2<NREAGENT;r2++)
269  {
270  //Get the local unknown and weight
271  //If the node is hanging
272  if(is_node2_hanging)
273  {
274  //Read out the local unknown from the master node
275  local_unknown =
276  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
277  c_nodal_index[r2]);
278  //Read out the hanging weight from the master node
279  hang_weight2 = hang_info2_pt->master_weight(m2);
280  }
281  //If the node is not hanging
282  else
283  {
284  //The local unknown number comes from the node
285  local_unknown =
286  this->nodal_local_eqn(l2,c_nodal_index[r2]);
287  //The hang weight is one
288  hang_weight2 = 1.0;
289  }
290 
291  //If the unknown is not pinned
292  if(local_unknown >= 0)
293  {
294  //Diagonal terms (i.e. the basic equations)
295  if(r2==r)
296  {
297  //Mass matrix term
298  jacobian(local_eqn,local_unknown)
299  -= T[r]*test(l)*psi(l2)*
300  node_pt(l2)->time_stepper_pt()->weight(1,0)*
301  W*hang_weight*hang_weight2;
302 
303  //Add the mass matrix term
304  if(flag==2)
305  {
306  mass_matrix(local_eqn,local_unknown)
307  += T[r]*test(l)*psi(l2)*W*hang_weight*hang_weight2;
308  }
309 
310  //Add contribution to Elemental Matrix
311  for(unsigned i=0;i<DIM;i++)
312  {
313  //Temporary term used in assembly
314  double tmp = wind[i];
315  if(!ALE_is_disabled_flag)
316  {tmp -= T[r]*mesh_velocity[i];}
317  //Now assemble Jacobian term
318  jacobian(local_eqn,local_unknown)
319  -=
320  dpsidx(l2,i)*(tmp*test(l)+D[r]*dtestdx(l,i))*
321  W*hang_weight*hang_weight2;
322  }
323 
324  } //End of diagonal terms
325 
326  //Now add the cross-reaction terms
327  jacobian(local_eqn,local_unknown) -=
328  dRdC(r,r2)*psi(l2)*test(l)*W*hang_weight*hang_weight2;
329  }
330  } //End of loop over reagents
331  } //End of loop over master nodes
332  } //End of loop over nodes
333  } //End of Jacobian calculation
334 
335  } //End of non-zero equation
336 
337  } //End of loop over reagents
338  } //End of loop over the master nodes for residual
339  } //End of loop over nodes
340 
341  } // End of loop over integration points
342 }
343 
344 
345 
346 //====================================================================
347 // Force build of templates
348 //====================================================================
349 ///One reagent
353 
357 
358 //Two reagents
362 
366 
370 
371 }
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
void fill_in_generic_residual_contribution_adv_diff_react(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add the element's contribution to the elemental residual vector.
cstr elem_len * i
Definition: cfortran.h:607
Refineable version of QAdvectionDiffusionReactionElement. Inherit from the standard QAdvectionDiffusi...
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
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