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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
27 
28 namespace oomph
29 {
30  //==========================================================================
31  /// Add the element's contribution to the elemental residual vector
32  /// and/or elemental jacobian matrix.
33  /// This function overloads the standard version so that the possible
34  /// presence of hanging nodes is taken into account.
35  //=========================================================================
36  template<unsigned NREAGENT, unsigned DIM>
39  Vector<double>& residuals,
40  DenseMatrix<double>& jacobian,
41  DenseMatrix<double>& mass_matrix,
42  unsigned flag)
43  {
44  // Find out how many nodes there are in the element
45  const unsigned n_node = nnode();
46 
47  // Get the nodal index at which the unknown is stored
48  unsigned c_nodal_index[NREAGENT];
49  for (unsigned r = 0; r < NREAGENT; r++)
50  {
51  c_nodal_index[r] = this->c_index_adv_diff_react(r);
52  }
53 
54  // Set up memory for the shape and test functions
55  Shape psi(n_node), test(n_node);
56  DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
57 
58  // Set the value of n_intpt
59  const unsigned n_intpt = integral_pt()->nweight();
60 
61  // Set the Vector to hold local coordinates
62  Vector<double> s(DIM);
63 
64  // Get diffusion coefficients
65  Vector<double> D = this->diff();
66 
67  // Get the timescales
68  Vector<double> T = this->tau();
69 
70  // Integers used to store the local equation number and local unknown
71  // indices for the residuals and jacobians
72  int local_eqn = 0, local_unknown = 0;
73 
74  // Local storage for pointers to hang_info objects
75  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
76 
77  // Local variable to determine the ALE stuff
78  bool ALE_is_disabled_flag = this->ALE_is_disabled;
79 
80  // Loop over the integration points
81  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
82  {
83  // Assign values of s
84  for (unsigned i = 0; i < DIM; i++) s[i] = integral_pt()->knot(ipt, i);
85 
86  // Get the integral weight
87  double w = integral_pt()->weight(ipt);
88 
89  // Call the derivatives of the shape and test functions
90  double J = this->dshape_and_dtest_eulerian_at_knot_adv_diff_react(
91  ipt, psi, dpsidx, test, dtestdx);
92 
93  // Premultiply the weights and the Jacobian
94  double W = w * J;
95 
96  // Calculate local values of the solution and its derivatives
97  // Allocate
98  Vector<double> interpolated_c(NREAGENT, 0.0);
99  Vector<double> dcdt(NREAGENT, 0.0);
100  Vector<double> interpolated_x(DIM, 0.0);
101  DenseMatrix<double> interpolated_dcdx(NREAGENT, DIM, 0.0);
102  Vector<double> mesh_velocity(DIM, 0.0);
103 
104 
105  // Calculate function value and derivatives:
106  // Loop over nodes
107  for (unsigned l = 0; l < n_node; l++)
108  {
109  // Loop over directions to calculate the position
110  for (unsigned j = 0; j < DIM; j++)
111  {
112  interpolated_x[j] += nodal_position(l, j) * psi(l);
113  }
114 
115  // Loop over the unknown reagents
116  for (unsigned r = 0; r < NREAGENT; r++)
117  {
118  // Get the value at the node
119  const double c_value = nodal_value(l, c_nodal_index[r]);
120 
121  // Calculate the interpolated value
122  interpolated_c[r] += c_value * psi(l);
123  dcdt[r] += this->dc_dt_adv_diff_react(l, r) * psi(l);
124 
125  // Loop over directions to calculate the derivatie
126  for (unsigned j = 0; j < DIM; j++)
127  {
128  interpolated_dcdx(r, j) += c_value * dpsidx(l, j);
129  }
130  }
131  }
132 
133  // Mesh velocity?
134  if (!ALE_is_disabled_flag)
135  {
136  for (unsigned l = 0; l < n_node; l++)
137  {
138  for (unsigned j = 0; j < DIM; j++)
139  {
140  mesh_velocity[j] += dnodal_position_dt(l, j) * psi(l);
141  }
142  }
143  }
144 
145  // Get source function
146  Vector<double> source(NREAGENT);
147  this->get_source_adv_diff_react(ipt, interpolated_x, source);
148 
149 
150  // Get wind
151  Vector<double> wind(DIM);
152  this->get_wind_adv_diff_react(ipt, s, interpolated_x, wind);
153 
154  // Get reaction terms
155  Vector<double> R(NREAGENT);
156  this->get_reaction_adv_diff_react(ipt, interpolated_c, R);
157 
158  // If we are getting the jacobian the get the derivative terms
159  DenseMatrix<double> dRdC(NREAGENT);
160  if (flag)
161  {
162  this->get_reaction_deriv_adv_diff_react(ipt, interpolated_c, dRdC);
163  }
164 
165  // Assemble residuals and Jacobian
166  //================================
167 
168  // Loop over the nodes for the test functions
169  for (unsigned l = 0; l < n_node; l++)
170  {
171  // Local variables to store the number of master nodes and
172  // the weight associated with the shape function if the node is hanging
173  unsigned n_master = 1;
174  double hang_weight = 1.0;
175  // Local bool (is the node hanging)
176  bool is_node_hanging = this->node_pt(l)->is_hanging();
177 
178  // If the node is hanging, get the number of master nodes
179  if (is_node_hanging)
180  {
181  hang_info_pt = this->node_pt(l)->hanging_pt();
182  n_master = hang_info_pt->nmaster();
183  }
184  // Otherwise there is just one master node, the node itself
185  else
186  {
187  n_master = 1;
188  }
189 
190  // Loop over the number of master nodes
191  for (unsigned m = 0; m < n_master; m++)
192  {
193  // Loop over the number of reagents
194  for (unsigned r = 0; r < NREAGENT; r++)
195  {
196  // Get the local equation number and hang_weight
197  // If the node is hanging
198  if (is_node_hanging)
199  {
200  // Read out the local equation from the master node
201  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
202  c_nodal_index[r]);
203  // Read out the weight from the master node
204  hang_weight = hang_info_pt->master_weight(m);
205  }
206  // If the node is not hanging
207  else
208  {
209  // The local equation number comes from the node itself
210  local_eqn = this->nodal_local_eqn(l, c_nodal_index[r]);
211  // The hang weight is one
212  hang_weight = 1.0;
213  }
214 
215  // If the nodal equation is not a boundary conditino
216  if (local_eqn >= 0)
217  {
218  // Add body force/source/reaction term and time derivative
219  residuals[local_eqn] -=
220  (T[r] * dcdt[r] + source[r] + R[r]) * test(l) * W * hang_weight;
221 
222  // The Advection Diffusion bit itself
223  for (unsigned k = 0; k < DIM; k++)
224  {
225  // Terms that multiply the test function
226  double tmp = wind[k];
227  // If the mesh is moving need to subtract the mesh velocity
228  if (!ALE_is_disabled_flag)
229  {
230  tmp -= T[r] * mesh_velocity[k];
231  }
232  // Now construct the contribution to the residuals
233  residuals[local_eqn] -= interpolated_dcdx(r, k) *
234  (tmp * test(l) + D[r] * dtestdx(l, k)) *
235  W * hang_weight;
236  }
237 
238  // Calculate the Jacobian
239  if (flag)
240  {
241  // Local variables to store the number of master nodes
242  // and the weights associated with each hanging node
243  unsigned n_master2 = 1;
244  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  // Loop over the reagents again
266  for (unsigned r2 = 0; r2 < NREAGENT; r2++)
267  {
268  // Get the local unknown and weight
269  // If the node is hanging
270  if (is_node2_hanging)
271  {
272  // Read out the local unknown from the master node
273  local_unknown = this->local_hang_eqn(
274  hang_info2_pt->master_node_pt(m2), c_nodal_index[r2]);
275  // Read out the hanging weight from the master node
276  hang_weight2 = hang_info2_pt->master_weight(m2);
277  }
278  // If the node is not hanging
279  else
280  {
281  // The local unknown number comes from the node
282  local_unknown =
283  this->nodal_local_eqn(l2, c_nodal_index[r2]);
284  // The hang weight is one
285  hang_weight2 = 1.0;
286  }
287 
288  // If the unknown is not pinned
289  if (local_unknown >= 0)
290  {
291  // Diagonal terms (i.e. the basic equations)
292  if (r2 == r)
293  {
294  // Mass matrix term
295  jacobian(local_eqn, local_unknown) -=
296  T[r] * test(l) * psi(l2) *
297  node_pt(l2)->time_stepper_pt()->weight(1, 0) * W *
298  hang_weight * hang_weight2;
299 
300  // Add the mass matrix term
301  if (flag == 2)
302  {
303  mass_matrix(local_eqn, local_unknown) +=
304  T[r] * test(l) * psi(l2) * W * hang_weight *
305  hang_weight2;
306  }
307 
308  // Add contribution to Elemental Matrix
309  for (unsigned i = 0; i < DIM; i++)
310  {
311  // Temporary term used in assembly
312  double tmp = wind[i];
313  if (!ALE_is_disabled_flag)
314  {
315  tmp -= T[r] * mesh_velocity[i];
316  }
317  // Now assemble Jacobian term
318  jacobian(local_eqn, local_unknown) -=
319  dpsidx(l2, i) *
320  (tmp * test(l) + D[r] * dtestdx(l, i)) * W *
321  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 *
329  hang_weight2;
330  }
331  } // End of loop over reagents
332  } // End of loop over master nodes
333  } // End of loop over nodes
334  } // End of Jacobian calculation
335 
336  } // End of non-zero equation
337 
338  } // End of loop over reagents
339  } // End of loop over the master nodes for residual
340  } // End of loop over nodes
341 
342  } // End of loop over integration points
343  }
344 
345 
346  //====================================================================
347  // Force build of templates
348  //====================================================================
349  /// One reagent
353 
357 
358  // Two reagents
362 
366 
370 
371 } // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
Class that contains data for hanging nodes.
Definition: nodes.h:742
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
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 and/or Jacobian matrix flag=1: comput...
Refineable version of QAdvectionDiffusionReactionElement. Inherit from the standard QAdvectionDiffusi...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...