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