refineable_space_time_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-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//====================================================================
26 // Source file for refineable Van der Pol elements
28 
29 /// //////////////////////////////////////////////////////////////////////
30 /// //////////////////////////////////////////////////////////////////////
31 /// //////////////////////////////////////////////////////////////////////
32 
33 namespace oomph
34 {
35  //========================================================================
36  /// Add element's contribution to the elemental residual vector and/or
37  /// Jacobian matrix.
38  /// flag=0: compute residual vector only
39  /// flag=1: compute both
40  //========================================================================
41  template<unsigned SPATIAL_DIM>
44  Vector<double>& residuals,
45  DenseMatrix<double>& jacobian,
46  const unsigned& flag)
47  {
48  // Find out how many nodes there are in the element
49  unsigned n_node = nnode();
50 
51  // Find the index at which the unknown is stored
52  unsigned u_nodal_index = this->u_index_ust_heat();
53 
54  // Set up memory for the shape functions
55  Shape psi(n_node);
56 
57  // Set up memory for the test functions
58  Shape test(n_node);
59 
60  // Allocate space for the derivatives of the shape functions
61  DShape dpsidx(n_node, SPATIAL_DIM + 1);
62 
63  // Allocate space for the derivatives of the test functions
64  DShape dtestdx(n_node, SPATIAL_DIM + 1);
65 
66  // Set the value of n_intpt
67  unsigned n_intpt = integral_pt()->nweight();
68 
69  // Storage for the local coordinates
70  Vector<double> s(SPATIAL_DIM + 1, 0.0);
71 
72  // Get the Alpha parameter
73  double alpha_local = this->alpha();
74 
75  // Get the Beta parameter
76  double beta_local = this->beta();
77 
78  // Integer to hold the local equation
79  int local_eqn = 0;
80 
81  // Integer to hold the local unknowns
82  int local_unknown = 0;
83 
84  // Local storage for pointer to hang info object
85  HangInfo* hang_info_pt = 0;
86 
87  // Local storage for pointer to (another) hang info object
88  HangInfo* hang_info2_pt = 0;
89 
90  // Loop over the integration points
91  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
92  {
93  // Assign values of s
94  for (unsigned i = 0; i < SPATIAL_DIM + 1; i++)
95  {
96  // Calculate the i-th local coordinate
97  s[i] = integral_pt()->knot(ipt, i);
98  }
99 
100  // Get the integral weight
101  double w = integral_pt()->weight(ipt);
102 
103  // Call the derivatives of the shape and test functions
104  double J = this->dshape_and_dtest_eulerian_at_knot_ust_heat(
105  ipt, psi, dpsidx, test, dtestdx);
106 
107  // Premultiply the weights and the Jacobian
108  double W = w * J;
109 
110  // Storage for the interpolated time value
111  double interpolated_t = 0.0;
112 
113  // Storage for the interpolated solution value
114  double interpolated_u = 0.0;
115 
116  // Storage for the interpolated time-derivative of the solution
117  double interpolated_dudt = 0.0;
118 
119  // Storage for the spatial coordinates
120  Vector<double> interpolated_x(SPATIAL_DIM, 0.0);
121 
122  // Storage for the spatial derivatives of the solution
123  Vector<double> interpolated_dudx(SPATIAL_DIM + 1, 0.0);
124 
125  // Storage for the mesh velocity
126  Vector<double> mesh_velocity(SPATIAL_DIM, 0.0);
127 
128  //-------------------------------------------------
129  // Calculate derivatives and source function value:
130  //-------------------------------------------------
131  // Loop over the nodes
132  for (unsigned l = 0; l < n_node; l++)
133  {
134  // Get the nodal value at the l-th node
135  double u_value = nodal_value(l, u_nodal_index);
136 
137  // Update the interpolated u value
138  interpolated_u += u_value * psi(l);
139 
140  // Loop over the coordinate directions (both spatial AND time)
141  for (unsigned j = 0; j < SPATIAL_DIM; j++)
142  {
143  // Update the interpolated x value
144  interpolated_x[j] += nodal_position(l, j) * psi(l);
145 
146  // Update the interpolated du/dx_j value
147  interpolated_dudx[j] += u_value * dpsidx(l, j);
148  }
149 
150  // Update the interpolated time value
151  interpolated_t += nodal_position(l, SPATIAL_DIM) * psi(l);
152 
153  // Update the interpolated du/dt value
154  interpolated_dudt += u_value * dpsidx(l, SPATIAL_DIM);
155  } // for (unsigned l=0;l<n_node;l++)
156 
157  // If ALE is enabled
158  if (!(this->ALE_is_disabled))
159  {
160  // Loop over the coordinate directions
161  for (unsigned j = 0; j < SPATIAL_DIM; j++)
162  {
163  // Loop over the nodes
164  for (unsigned l = 0; l < n_node; l++)
165  {
166  // Update the mesh velocity
167  mesh_velocity[j] += nodal_position(l, j) * dpsidx(l, SPATIAL_DIM);
168  }
169  } // for (unsigned j=0;j<SPATIAL_DIM;j++)
170  } // if (!ALE_is_disabled)
171 
172  // Initialise the source term value
173  double source = 0.0;
174 
175  // Get the interpolated source term value
176  this->get_source_ust_heat(interpolated_t, ipt, interpolated_x, source);
177 
178  //--------------------------------
179  // Assemble residuals and Jacobian
180  //--------------------------------
181  // Loop over the nodes (or equivalently the test functions)
182  for (unsigned l = 0; l < n_node; l++)
183  {
184  // Storage for the number of master nodes
185  unsigned n_master = 1;
186 
187  // Storage for the hang weight associated with the shape function
188  double hang_weight = 1.0;
189 
190  // Check if the node is hanging
191  bool is_node_hanging = this->node_pt(l)->is_hanging();
192 
193  // If the node is hanging, get the number of master nodes
194  if (is_node_hanging)
195  {
196  // Get the hang info pointer associated with the l-th node in the
197  // element
198  hang_info_pt = this->node_pt(l)->hanging_pt();
199 
200  // Get the number of master nodes associated with the node
201  n_master = hang_info_pt->nmaster();
202  }
203  // If it's not hanging there is just one master node; the node itself
204  else
205  {
206  // Set the number of master nodes to one
207  n_master = 1;
208  }
209 
210  // Loop over the number of master nodes
211  for (unsigned m = 0; m < n_master; m++)
212  {
213  // Check if the node is hanging
214  if (is_node_hanging)
215  {
216  // Read out the local equation from the master node
217  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
218  u_nodal_index);
219 
220  // Read out the weight from the master node
221  hang_weight = hang_info_pt->master_weight(m);
222  }
223  // If the node is not hanging
224  else
225  {
226  // The local equation number comes from the node itself
227  local_eqn = this->nodal_local_eqn(l, u_nodal_index);
228 
229  // The hang weight is one
230  hang_weight = 1.0;
231  }
232 
233  // If the nodal equation is not a boundary condition
234  if (local_eqn >= 0)
235  {
236  // Add source term and time derivative
237  residuals[local_eqn] +=
238  ((source + alpha_local * interpolated_dudt) * test(l) * W *
239  hang_weight);
240 
241  // If ALE is enabled
242  if (!(this->ALE_is_disabled))
243  {
244  // Loop over the coordinate directions
245  for (unsigned k = 0; k < SPATIAL_DIM; k++)
246  {
247  // Add in the mesh velocity contribution
248  residuals[local_eqn] -=
249  (alpha_local * mesh_velocity[k] * interpolated_dudx[k] *
250  test(l) * W * hang_weight);
251  }
252  } // if (!ALE_is_disabled)
253 
254  // Loop over the coordinate directions
255  for (unsigned k = 0; k < SPATIAL_DIM; k++)
256  {
257  // Add in the contribution from the Laplace operator
258  residuals[local_eqn] += (beta_local * interpolated_dudx[k] *
259  dtestdx(l, k) * W * hang_weight);
260  }
261 
262  //-----------------------
263  // Calculate the Jacobian
264  //-----------------------
265  // If we also need to construct the Jacobian
266  if (flag)
267  {
268  // Storage for the number of master nodes
269  unsigned n_master2 = 1;
270 
271  // Storage for the hang weight associated with the shape function
272  double hang_weight2 = 1.0;
273 
274  // Loop over the nodes for the variables
275  for (unsigned l2 = 0; l2 < n_node; l2++)
276  {
277  // Check if the node is hanging
278  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
279 
280  // If the node is hanging, get the number of master nodes
281  if (is_node2_hanging)
282  {
283  // Get the hang info pointer associated with the l2-th node
284  hang_info2_pt = this->node_pt(l2)->hanging_pt();
285 
286  // Get the number of master nodes associated with the node
287  n_master2 = hang_info2_pt->nmaster();
288  }
289  // If it's not hanging there is just one master node; the node
290  // itself
291  else
292  {
293  // Set the number of master nodes to one
294  n_master2 = 1;
295  }
296 
297  // Loop over the master nodes
298  for (unsigned m2 = 0; m2 < n_master2; m2++)
299  {
300  // Check if the node is hanging
301  if (is_node2_hanging)
302  {
303  // Read out the local equation from the master node
304  local_unknown = this->local_hang_eqn(
305  hang_info2_pt->master_node_pt(m2), u_nodal_index);
306 
307  // Read out the weight from the master node
308  hang_weight2 = hang_info2_pt->master_weight(m2);
309  }
310  // If the node is not hanging
311  else
312  {
313  // The local equation number comes from the node itself
314  local_unknown = this->nodal_local_eqn(l2, u_nodal_index);
315 
316  // The hang weight is one
317  hang_weight2 = 1.0;
318  }
319 
320  // If the unknown is not pinned
321  if (local_unknown >= 0)
322  {
323  // Time-derivative contribution
324  jacobian(local_eqn, local_unknown) +=
325  (alpha_local * test(l) * dpsidx(l2, SPATIAL_DIM) * W *
326  hang_weight * hang_weight2);
327 
328  // Laplace operator contribution
329  for (unsigned k = 0; k < SPATIAL_DIM; k++)
330  {
331  // Add the Laplacian contribution to the elemental
332  // Jacobian
333  jacobian(local_eqn, local_unknown) +=
334  (beta_local * dpsidx(l2, k) * dtestdx(l, k) * W *
335  hang_weight * hang_weight2);
336  }
337 
338  // If ALE is enabled
339  if (!(this->ALE_is_disabled))
340  {
341  // Loop over the spatial coordinates
342  for (unsigned k = 0; k < SPATIAL_DIM; k++)
343  {
344  // Add the ALE contribution to the Jacobian
345  jacobian(local_eqn, local_unknown) -=
346  (alpha_local * mesh_velocity[k] * dpsidx(l2, k) *
347  test(l) * W * hang_weight * hang_weight2);
348  }
349  } // if (!(this->ALE_is_disabled))
350  } // if (local_unknown>=0)
351  } // for (unsigned m2=0;m2<n_master;m2++)
352  } // for (unsigned l2=0;l2<n_node;l2++)
353  } // if (flag)
354  } // if (local_eqn>=0)
355  } // for (unsigned m=0;m<n_master;m++)
356  } // for (unsigned l=0;l<n_node;l++)
357  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
358  } // End of fill_in_generic_residual_contribution_ust_heat
359 
360  //=========================================================================
361  // Force build of templates:
362  // Build 2D (in space) elements --> 3D (space-time) elements
363  //=========================================================================
367 } // 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
Refineable version of 2D QUnsteadyHeatSpaceTimeElement elements.
virtual void fill_in_generic_residual_contribution_ust_heat(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=0: compute residu...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...