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-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  //=========================================================================
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
57  Vector<double> s(2);
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  //--------
139  Vector<double> wind(3);
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
148  DenseMatrix<double> D(3, 3);
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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
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
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
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.
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
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...