refineable_time_harmonic_linear_elasticity_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 // Non-inline member functions and static member data for refineable
27 // linear elasticity elements
28 
29 
30 #include "refineable_time_harmonic_linear_elasticity_elements.h"
31 
32 namespace oomph
33 {
34  //====================================================================
35  /// Residuals for Refineable QTimeHarmonicLinearElasticityElements
36  //====================================================================
37  template<unsigned DIM>
38  void RefineableTimeHarmonicLinearElasticityEquations<DIM>::
39  fill_in_generic_contribution_to_residuals_time_harmonic_linear_elasticity(
40  Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
41  {
42 #ifdef PARANOID
43 
44  // Find out how many positional dofs there are
45  unsigned n_position_type = this->nnodal_position_type();
46  if (n_position_type != 1)
47  {
48  throw OomphLibError("TimeHarmonicLinearElasticity is not yet implemented "
49  "for more than one position type",
50  OOMPH_CURRENT_FUNCTION,
51  OOMPH_EXCEPTION_LOCATION);
52  }
53 
54  // Throw and error if an elasticity tensor has not been set
55  if (this->Elasticity_tensor_pt == 0)
56  {
57  throw OomphLibError("No elasticity tensor set",
58  OOMPH_CURRENT_FUNCTION,
59  OOMPH_EXCEPTION_LOCATION);
60  }
61 
62 #endif
63 
64  // Find out how many nodes there are
65  unsigned n_node = this->nnode();
66 
67  // Find the indices at which the local displacements are stored
68  std::complex<unsigned> u_nodal_index[DIM];
69  for (unsigned i = 0; i < DIM; i++)
70  {
71  u_nodal_index[i] = this->u_index_time_harmonic_linear_elasticity(i);
72  }
73 
74  // Square of non-dimensional frequency
75  const double omega_sq_local = this->omega_sq();
76 
77  // Set up memory for the shape functions
78  Shape psi(n_node);
79  DShape dpsidx(n_node, DIM);
80 
81  // Set the value of Nintpt -- the number of integration points
82  unsigned n_intpt = this->integral_pt()->nweight();
83 
84  // Set the vector to hold the local coordinates in the element
85  Vector<double> s(DIM);
86 
87  // Integer to store the local equation number
88  int local_eqn = 0, local_unknown = 0;
89 
90  // Local storage for pointers to hang_info objects
91  HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
92 
93  // Loop over the integration points
94  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
95  {
96  // Assign the values of s
97  for (unsigned i = 0; i < DIM; ++i)
98  {
99  s[i] = this->integral_pt()->knot(ipt, i);
100  }
101 
102  // Get the integral weight
103  double w = this->integral_pt()->weight(ipt);
104 
105  // Call the derivatives of the shape functions (and get Jacobian)
106  double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
107 
108  // Storage for Eulerian coordinates (initialised to zero)
109  Vector<double> interpolated_x(DIM, 0.0);
110 
111  // Displacement
112  Vector<std::complex<double>> interpolated_u(
113  DIM, std::complex<double>(0.0, 0.0));
114 
115  // Calculate interpolated values of the derivative of displacements
116  DenseMatrix<std::complex<double>> interpolated_dudx(
117  DIM, DIM, std::complex<double>(0.0, 0.0));
118 
119  // Calculate displacements and derivatives
120  for (unsigned l = 0; l < n_node; l++)
121  {
122  // Loop over displacement components
123  for (unsigned i = 0; i < DIM; i++)
124  {
125  // Calculate the coordinates and the accelerations
126  interpolated_x[i] += this->nodal_position(l, i) * psi(l);
127 
128  // Get the nodal displacements
129  const std::complex<double> u_value =
130  std::complex<double>(nodal_value(l, u_nodal_index[i].real()),
131  nodal_value(l, u_nodal_index[i].imag()));
132 
133  interpolated_u[i] += u_value * psi(l);
134 
135  // Loop over derivative directions
136  for (unsigned j = 0; j < DIM; j++)
137  {
138  interpolated_dudx(i, j) += u_value * dpsidx(l, j);
139  }
140  }
141  }
142 
143  // Get body force at current time
144  Vector<std::complex<double>> b(DIM);
145  this->body_force(interpolated_x, b);
146 
147  // Premultiply the weights and the Jacobian
148  double W = w * J;
149 
150  // Number of master nodes and storage for the weight of the shape function
151  unsigned n_master = 1;
152  double hang_weight = 1.0;
153 
154  // Loop over the test functions, nodes of the element
155  for (unsigned l = 0; l < n_node; l++)
156  {
157  // Local boolean to indicate whether the node is hanging
158  bool is_node_hanging = node_pt(l)->is_hanging();
159 
160  // If the node is hanging
161  if (is_node_hanging)
162  {
163  hang_info_pt = node_pt(l)->hanging_pt();
164  // Read out number of master nodes from hanging data
165  n_master = hang_info_pt->nmaster();
166  }
167  // Otherwise the node is its own master
168  else
169  {
170  n_master = 1;
171  }
172 
173  // Loop over the master nodes
174  for (unsigned m = 0; m < n_master; m++)
175  {
176  // Loop over the displacement components
177  for (unsigned a = 0; a < DIM; a++)
178  {
179  // Real
180  //-----
181 
182  // Get the equation number
183  if (is_node_hanging)
184  {
185  // Get the equation number from the master node
186  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
187  u_nodal_index[a].real());
188  // Get the hang weight from the master node
189  hang_weight = hang_info_pt->master_weight(m);
190  }
191  // Otherwise the node is not hanging
192  else
193  {
194  local_eqn = this->nodal_local_eqn(l, u_nodal_index[a].real());
195  hang_weight = 1.0;
196  }
197 
198  /*IF it's not a boundary condition*/
199  if (local_eqn >= 0)
200  {
201  // Acceleration and body force
202  residuals[local_eqn] +=
203  (-omega_sq_local * interpolated_u[a].real() - b[a].real()) *
204  psi(l) * W * hang_weight;
205 
206  // Stress term
207  for (unsigned b = 0; b < DIM; b++)
208  {
209  for (unsigned c = 0; c < DIM; c++)
210  {
211  for (unsigned d = 0; d < DIM; d++)
212  {
213  // Add the stress terms to the residuals
214  residuals[local_eqn] += this->E(a, b, c, d) *
215  interpolated_dudx(c, d).real() *
216  dpsidx(l, b) * W * hang_weight;
217  }
218  }
219  }
220 
221  // Jacobian entries
222  if (flag)
223  {
224  // Number of master nodes and weights
225  unsigned n_master2 = 1;
226  double hang_weight2 = 1.0;
227  // Loop over the displacement basis functions again
228  for (unsigned l2 = 0; l2 < n_node; l2++)
229  {
230  // Local boolean to indicate whether the node is hanging
231  bool is_node2_hanging = node_pt(l2)->is_hanging();
232 
233  // If the node is hanging
234  if (is_node2_hanging)
235  {
236  hang_info2_pt = node_pt(l2)->hanging_pt();
237  // Read out number of master nodes from hanging data
238  n_master2 = hang_info2_pt->nmaster();
239  }
240  // Otherwise the node is its own master
241  else
242  {
243  n_master2 = 1;
244  }
245 
246  // Loop over the master nodes
247  for (unsigned m2 = 0; m2 < n_master2; m2++)
248  {
249  // Loop over the displacement components again
250  for (unsigned c = 0; c < DIM; c++)
251  {
252  // Get the number of the unknown
253  // If the node is hanging
254  if (is_node2_hanging)
255  {
256  // Get the equation number from the master node
257  local_unknown = this->local_hang_eqn(
258  hang_info2_pt->master_node_pt(m2),
259  u_nodal_index[c].real());
260  // Get the hang weights
261  hang_weight2 = hang_info2_pt->master_weight(m2);
262  }
263  else
264  {
265  local_unknown =
266  this->nodal_local_eqn(l2, u_nodal_index[c].real());
267  hang_weight2 = 1.0;
268  }
269 
270  // If it's not pinned
271  if (local_unknown >= 0)
272  {
273  // Inertial term
274  if (a == c)
275  {
276  jacobian(local_eqn, local_unknown) -=
277  omega_sq_local * psi(l) * psi(l2) * W *
278  hang_weight * hang_weight2;
279  }
280 
281  // Stress term
282  for (unsigned b = 0; b < DIM; b++)
283  {
284  for (unsigned d = 0; d < DIM; d++)
285  {
286  // Add the contribution to the Jacobian matrix
287  jacobian(local_eqn, local_unknown) +=
288  this->E(a, b, c, d) * dpsidx(l2, d) *
289  dpsidx(l, b) * W * hang_weight * hang_weight2;
290  }
291  }
292  } // End of if not boundary condition
293  }
294  }
295  }
296  } // End of jacobian calculation
297 
298  } // End of if not boundary condition
299 
300 
301  // Imag
302  //-----
303 
304  // Get the equation number
305  if (is_node_hanging)
306  {
307  // Get the equation number from the master node
308  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
309  u_nodal_index[a].imag());
310  // Get the hang weight from the master node
311  hang_weight = hang_info_pt->master_weight(m);
312  }
313  // Otherwise the node is not hanging
314  else
315  {
316  local_eqn = this->nodal_local_eqn(l, u_nodal_index[a].imag());
317  hang_weight = 1.0;
318  }
319 
320  /*IF it's not a boundary condition*/
321  if (local_eqn >= 0)
322  {
323  // Acceleration and body force
324  residuals[local_eqn] +=
325  (-omega_sq_local * interpolated_u[a].imag() - b[a].imag()) *
326  psi(l) * W * hang_weight;
327 
328  // Stress term
329  for (unsigned b = 0; b < DIM; b++)
330  {
331  for (unsigned c = 0; c < DIM; c++)
332  {
333  for (unsigned d = 0; d < DIM; d++)
334  {
335  // Add the stress terms to the residuals
336  residuals[local_eqn] += this->E(a, b, c, d) *
337  interpolated_dudx(c, d).imag() *
338  dpsidx(l, b) * W * hang_weight;
339  }
340  }
341  }
342 
343  // Jacobian entries
344  if (flag)
345  {
346  // Number of master nodes and weights
347  unsigned n_master2 = 1;
348  double hang_weight2 = 1.0;
349  // Loop over the displacement basis functions again
350  for (unsigned l2 = 0; l2 < n_node; l2++)
351  {
352  // Local boolean to indicate whether the node is hanging
353  bool is_node2_hanging = node_pt(l2)->is_hanging();
354 
355  // If the node is hanging
356  if (is_node2_hanging)
357  {
358  hang_info2_pt = node_pt(l2)->hanging_pt();
359  // Read out number of master nodes from hanging data
360  n_master2 = hang_info2_pt->nmaster();
361  }
362  // Otherwise the node is its own master
363  else
364  {
365  n_master2 = 1;
366  }
367 
368  // Loop over the master nodes
369  for (unsigned m2 = 0; m2 < n_master2; m2++)
370  {
371  // Loop over the displacement components again
372  for (unsigned c = 0; c < DIM; c++)
373  {
374  // Get the number of the unknown
375  // If the node is hanging
376  if (is_node2_hanging)
377  {
378  // Get the equation number from the master node
379  local_unknown = this->local_hang_eqn(
380  hang_info2_pt->master_node_pt(m2),
381  u_nodal_index[c].imag());
382  // Get the hang weights
383  hang_weight2 = hang_info2_pt->master_weight(m2);
384  }
385  else
386  {
387  local_unknown =
388  this->nodal_local_eqn(l2, u_nodal_index[c].imag());
389  hang_weight2 = 1.0;
390  }
391 
392  // If it's not pinned
393  if (local_unknown >= 0)
394  {
395  // Inertial term
396  if (a == c)
397  {
398  jacobian(local_eqn, local_unknown) -=
399  omega_sq_local * psi(l) * psi(l2) * W *
400  hang_weight * hang_weight2;
401  }
402 
403  // Stress term
404  for (unsigned b = 0; b < DIM; b++)
405  {
406  for (unsigned d = 0; d < DIM; d++)
407  {
408  // Add the contribution to the Jacobian matrix
409  jacobian(local_eqn, local_unknown) +=
410  this->E(a, b, c, d) * dpsidx(l2, d) *
411  dpsidx(l, b) * W * hang_weight * hang_weight2;
412  }
413  }
414  } // End of if not boundary condition
415  }
416  }
417  }
418  } // End of jacobian calculation
419 
420  } // End of if not boundary condition
421 
422 
423  } // End of loop over coordinate directions
424  }
425  } // End of loop over shape functions
426  } // End of loop over integration points
427  }
428 
429 
430  //====================================================================
431  /// Force building of required templates
432  //====================================================================
433  template class RefineableTimeHarmonicLinearElasticityEquations<2>;
434  template class RefineableTimeHarmonicLinearElasticityEquations<3>;
435 
436 } // namespace oomph
s
static char t char * s
Definition: cfortran.h:568
i
cstr elem_len * i
Definition: cfortran.h:603
oomph::DShape
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
oomph::HangInfo
Class that contains data for hanging nodes.
Definition: nodes.h:742
oomph::HangInfo::master_weight
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
oomph::HangInfo::master_node_pt
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition: oomph_definitions.h:222
oomph::RefineableTimeHarmonicLinearElasticityEquations::fill_in_generic_contribution_to_residuals_time_harmonic_linear_elasticity
void fill_in_generic_contribution_to_residuals_time_harmonic_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overloaded helper function to take hanging nodes into account.
Definition: refineable_time_harmonic_linear_elasticity_elements.cc:39
oomph::Shape
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30
refineable_time_harmonic_linear_elasticity_elements.h