refineable_unsteady_heat_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
28
29namespace oomph
30{
31 //========================================================================
32 /// Add element's contribution to the elemental
33 /// residual vector and/or Jacobian matrix.
34 /// flag=1: compute both
35 /// flag=0: compute only residual vector
36 //========================================================================
37 template<unsigned DIM>
40 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
41 {
42 // Find out how many nodes there are in the element
43 unsigned n_node = nnode();
44
45 // Get continuous time from timestepper of first node
46 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
47
48 // Find the index at which the unknown is stored
49 unsigned u_nodal_index = this->u_index_ust_heat();
50
51 // Set up memory for the shape and test functions
52 Shape psi(n_node), test(n_node);
53 DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
54
55 // Set the value of n_intpt
56 unsigned n_intpt = integral_pt()->nweight();
57
58 // Set the Vector to hold local coordinates
59 Vector<double> s(DIM);
60
61 // Get Alpha and beta parameters number
62 double alpha_local = this->alpha();
63 double beta_local = this->beta();
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 < DIM; 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
85 double J = this->dshape_and_dtest_eulerian_at_knot_ust_heat(
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
92 double dudt = 0.0;
93 double interpolated_u = 0.0;
94
95 // This needs to be a Vector to be ANSI C++, initialise to zero
96 Vector<double> interpolated_x(DIM, 0.0);
97 Vector<double> interpolated_dudx(DIM, 0.0);
98 Vector<double> mesh_velocity(DIM, 0.0);
99
100
101 // Calculate function value and derivatives:
102 //-----------------------------------------
103
104 // Loop over nodes
105 for (unsigned l = 0; l < n_node; l++)
106 {
107 // Get the value of the unknown at the node
108 double u_value = this->nodal_value(l, u_nodal_index);
109 interpolated_u += u_value * psi(l);
110 dudt += this->du_dt_ust_heat(l) * psi(l);
111 // Loop over directions
112 for (unsigned j = 0; j < DIM; j++)
113 {
114 interpolated_x[j] += nodal_position(l, j) * psi(l);
115 interpolated_dudx[j] += u_value * dpsidx(l, j);
116 }
117 }
118
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 < DIM; j++)
125 {
126 mesh_velocity[j] += dnodal_position_dt(l, j) * psi(l);
127 }
128 }
129 }
130
131 // Get body force
132 double source;
133 this->get_source_ust_heat(time, ipt, interpolated_x, source);
134
135
136 // Assemble residuals and Jacobian
137 //================================
138
139 // Loop over the nodes for the test functions
140 //-------------------------------------------
141 for (unsigned l = 0; l < n_node; l++)
142 {
143 // Local variables to store the number of master nodes and
144 // the weight associated with the shape function if the node is hanging
145 unsigned n_master = 1;
146 double hang_weight = 1.0;
147 // Local bool (is the node hanging)
148 bool is_node_hanging = this->node_pt(l)->is_hanging();
149
150
151 // If the node is hanging, get the number of master nodes
152 if (is_node_hanging)
153 {
154 hang_info_pt = this->node_pt(l)->hanging_pt();
155 n_master = hang_info_pt->nmaster();
156 }
157 // Otherwise there is just one master node, the node itself
158 else
159 {
160 n_master = 1;
161 }
162
163 // Loop over the number of master nodes
164 for (unsigned m = 0; m < n_master; m++)
165 {
166 // Get the local equation number and hang_weight
167 // If the node is hanging
168 if (is_node_hanging)
169 {
170 // Read out the local equation from the master node
171 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
172 u_nodal_index);
173 // Read out the weight from the master node
174 hang_weight = hang_info_pt->master_weight(m);
175 }
176 // If the node is not hanging
177 else
178 {
179 // The local equation number comes from the node itself
180 local_eqn = this->nodal_local_eqn(l, u_nodal_index);
181 // The hang weight is one
182 hang_weight = 1.0;
183 }
184
185 // If the nodal equation is not a boundary condition
186 if (local_eqn >= 0)
187 {
188 // Add body force/source term and time derivative
189 residuals[local_eqn] +=
190 (alpha_local * dudt + source) * test(l) * W * hang_weight;
191
192 // Mesh velocity and Laplace operator itself
193 for (unsigned k = 0; k < DIM; k++)
194 {
195 double tmp = dtestdx(l, k) * beta_local;
196 if (!ALE_is_disabled_flag)
197 tmp -= alpha_local * mesh_velocity[k] * test(l);
198 residuals[local_eqn] +=
199 interpolated_dudx[k] * tmp * W * hang_weight;
200 }
201
202 // Calculate the Jacobian
203 if (flag)
204 {
205 // Local variables to store the number of master nodes
206 // and the weights associated with each hanging node
207 unsigned n_master2 = 1;
208 double hang_weight2 = 1.0;
209 // Loop over the nodes for the variables
210 for (unsigned l2 = 0; l2 < n_node; l2++)
211 {
212 // Local bool (is the node hanging)
213 bool is_node2_hanging = this->node_pt(l2)->is_hanging();
214 // If the node is hanging, get the number of master nodes
215 if (is_node2_hanging)
216 {
217 hang_info2_pt = this->node_pt(l2)->hanging_pt();
218 n_master2 = hang_info2_pt->nmaster();
219 }
220 // Otherwise there is one master node, the node itself
221 else
222 {
223 n_master2 = 1;
224 }
225
226 // Loop over the master nodes
227 for (unsigned m2 = 0; m2 < n_master2; m2++)
228 {
229 // Get the local unknown and weight
230 // If the node is hanging
231 if (is_node2_hanging)
232 {
233 // Read out the local unknown from the master node
234 local_unknown = this->local_hang_eqn(
235 hang_info2_pt->master_node_pt(m2), u_nodal_index);
236 // Read out the hanging weight from the master node
237 hang_weight2 = hang_info2_pt->master_weight(m2);
238 }
239 // If the node is not hanging
240 else
241 {
242 // The local unknown number comes from the node
243 local_unknown = this->nodal_local_eqn(l2, u_nodal_index);
244 // The hang weight is one
245 hang_weight2 = 1.0;
246 }
247
248 // If the unknown is not pinned
249 if (local_unknown >= 0)
250 {
251 // Add contribution to Elemental Matrix
252 // Mass matrix
253 jacobian(local_eqn, local_unknown) +=
254 alpha_local * test(l) * psi(l2) *
255 this->node_pt(l2)->time_stepper_pt()->weight(1, 0) * W *
256 hang_weight * hang_weight2;
257
258 // Laplace operator and mesh veloc bit
259 for (unsigned k = 0; k < DIM; k++)
260 {
261 double tmp = dtestdx(l, k) * beta_local;
262 if (!ALE_is_disabled_flag)
263 {
264 tmp -= alpha_local * mesh_velocity[k] * test(l);
265 }
266 jacobian(local_eqn, local_unknown) +=
267 dpsidx(l2, k) * tmp * W * hang_weight * hang_weight2;
268 }
269 }
270 } // End of loop over master nodes
271 }
272 } // End of Jacobian calculation
273 }
274 } // End of loop over master nodes for residuals
275 } // End of loop over nodes
276
277 } // End of loop over integration points
278 }
279
280 //====================================================================
281 // Force build of templates
282 //====================================================================
286
290
291} // 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
Refineable version of 2D QUnsteadyHeatElement elements.
void fill_in_generic_residual_contribution_ust_heat(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...