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-2022 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
28namespace 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
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
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...