linear_elasticity_traction_elements.h
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-2024 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 // Header file for elements that are used to apply surface loads to
27 // the equations of linear elasticity
28 
29 #ifndef OOMPH_LINEAR_ELASTICITY_TRACTION_ELEMENTS_HEADER
30 #define OOMPH_LINEAR_ELASTICITY_TRACTION_ELEMENTS_HEADER
31 
32 // Config header generated by autoconfig
33 #ifdef HAVE_CONFIG_H
34 #include <oomph-lib-config.h>
35 #endif
36 
37 
38 // OOMPH-LIB headers
39 #include "../generic/Qelements.h"
40 
41 namespace oomph
42 {
43  //=======================================================================
44  /// Namespace containing the zero traction function for linear elasticity
45  /// traction elements
46  //=======================================================================
47  namespace LinearElasticityTractionElementHelper
48  {
49  //=======================================================================
50  /// Default load function (zero traction)
51  //=======================================================================
52  void Zero_traction_fct(const double& time,
53  const Vector<double>& x,
54  const Vector<double>& N,
55  Vector<double>& load);
56 
57  } // namespace LinearElasticityTractionElementHelper
58 
59 
60  //======================================================================
61  /// A class for elements that allow the imposition of an applied traction
62  /// in the equations of linear elasticity.
63  /// The geometrical information can be read from the FaceGeometry<ELEMENT>
64  /// class and and thus, we can be generic enough without the need to have
65  /// a separate equations class.
66  //======================================================================
67  template<class ELEMENT>
68  class LinearElasticityTractionElement : public virtual FaceGeometry<ELEMENT>,
69  public virtual FaceElement
70  {
71  protected:
72  /// Index at which the i-th displacement component is stored
74 
75  /// Pointer to an imposed traction function. Arguments:
76  /// Eulerian coordinate; outer unit normal;
77  /// applied traction. (Not all of the input arguments will be
78  /// required for all specific load functions but the list should
79  /// cover all cases)
80  void (*Traction_fct_pt)(const double& time,
81  const Vector<double>& x,
82  const Vector<double>& n,
83  Vector<double>& result);
84 
85 
86  /// Get the traction vector: Pass number of integration point
87  /// (dummy), Eulerlian coordinate and normal vector and return the load
88  /// vector (not all of the input arguments will be required for all specific
89  /// load functions but the list should cover all cases). This function is
90  /// virtual so it can be overloaded for FSI.
91  virtual void get_traction(const double& time,
92  const unsigned& intpt,
93  const Vector<double>& x,
94  const Vector<double>& n,
96  {
97  Traction_fct_pt(time, x, n, traction);
98  }
99 
100 
101  /// Helper function that actually calculates the residuals
102  // This small level of indirection is required to avoid calling
103  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
104  // which causes all kinds of pain if overloading later on
106  Vector<double>& residuals);
107 
108 
109  public:
110  /// Constructor, which takes a "bulk" element and the
111  /// value of the index and its limit
113  const int& face_index)
114  : FaceGeometry<ELEMENT>(), FaceElement()
115  {
116  // Attach the geometrical information to the element. N.B. This function
117  // also assigns nbulk_value from the required_nvalue of the bulk element
118  element_pt->build_face_element(face_index, this);
119 
120 #ifdef PARANOID
121  {
122  // Check that the element is not a refineable 3d element
123  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
124  // If it's three-d
125  if (elem_pt->dim() == 3)
126  {
127  // Is it refineable
128  RefineableElement* ref_el_pt =
129  dynamic_cast<RefineableElement*>(elem_pt);
130  if (ref_el_pt != 0)
131  {
132  if (this->has_hanging_nodes())
133  {
134  throw OomphLibError("This flux element will not work correctly "
135  "if nodes are hanging\n",
136  OOMPH_CURRENT_FUNCTION,
137  OOMPH_EXCEPTION_LOCATION);
138  }
139  }
140  }
141  }
142 #endif
143 
144  // Find the dimension of the problem
145  unsigned n_dim = element_pt->nodal_dimension();
146 
147  // Find the index at which the displacemenet unknowns are stored
148  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
149  this->U_index_linear_elasticity_traction.resize(n_dim);
150  for (unsigned i = 0; i < n_dim; i++)
151  {
152  this->U_index_linear_elasticity_traction[i] =
153  cast_element_pt->u_index_linear_elasticity(i);
154  }
155 
156  // Zero traction
159  }
160 
161 
162  /// Reference to the traction function pointer
163  void (*&traction_fct_pt())(const double& time,
164  const Vector<double>& x,
165  const Vector<double>& n,
167  {
168  return Traction_fct_pt;
169  }
170 
171 
172  /// Return the residuals
174  {
176  }
177 
178 
179  /// Fill in contribution from Jacobian
181  DenseMatrix<double>& jacobian)
182  {
183  // Call the residuals
185  }
186 
187  /// Specify the value of nodal zeta from the face geometry
188  /// The "global" intrinsic coordinate of the element when
189  /// viewed as part of a geometric object should be given by
190  /// the FaceElement representation, by default (needed to break
191  /// indeterminacy if bulk element is SolidElement)
192  double zeta_nodal(const unsigned& n,
193  const unsigned& k,
194  const unsigned& i) const
195  {
196  return FaceElement::zeta_nodal(n, k, i);
197  }
198 
199  /// Output function
200  void output(std::ostream& outfile)
201  {
202  FiniteElement::output(outfile);
203  }
204 
205  /// Output function
206  void output(std::ostream& outfile, const unsigned& n_plot)
207  {
208  FiniteElement::output(outfile, n_plot);
209  }
210 
211  /// C_style output function
212  void output(FILE* file_pt)
213  {
214  FiniteElement::output(file_pt);
215  }
216 
217  /// C-style output function
218  void output(FILE* file_pt, const unsigned& n_plot)
219  {
220  FiniteElement::output(file_pt, n_plot);
221  }
222 
223 
224  /// Compute traction vector at specified local coordinate
225  /// Should only be used for post-processing; ignores dependence
226  /// on integration point!
227  void traction(const double& time,
228  const Vector<double>& s,
230  };
231 
232  /// ////////////////////////////////////////////////////////////////////
233  /// ////////////////////////////////////////////////////////////////////
234  /// ////////////////////////////////////////////////////////////////////
235 
236  //=====================================================================
237  /// Compute traction vector at specified local coordinate
238  /// Should only be used for post-processing; ignores dependence
239  /// on integration point!
240  //=====================================================================
241  template<class ELEMENT>
243  const double& time, const Vector<double>& s, Vector<double>& traction)
244  {
245  unsigned n_dim = this->nodal_dimension();
246 
247  // Position vector
248  Vector<double> x(n_dim);
249  interpolated_x(s, x);
250 
251  // Outer unit normal
252  Vector<double> unit_normal(n_dim);
253  outer_unit_normal(s, unit_normal);
254 
255  // Dummy
256  unsigned ipt = 0;
257 
258  // Traction vector
259  get_traction(time, ipt, x, unit_normal, traction);
260  }
261 
262 
263  //=====================================================================
264  /// Return the residuals for the LinearElasticityTractionElement equations
265  //=====================================================================
266  template<class ELEMENT>
269  Vector<double>& residuals)
270  {
271  // Find out how many nodes there are
272  unsigned n_node = nnode();
273 
274  // Get continuous time from timestepper of first node
275  double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
276 
277 #ifdef PARANOID
278  // Find out how many positional dofs there are
279  unsigned n_position_type = this->nnodal_position_type();
280  if (n_position_type != 1)
281  {
282  throw OomphLibError("LinearElasticity is not yet implemented for more "
283  "than one position type",
284  OOMPH_CURRENT_FUNCTION,
285  OOMPH_EXCEPTION_LOCATION);
286  }
287 #endif
288 
289  // Find out the dimension of the node
290  unsigned n_dim = this->nodal_dimension();
291 
292  // Cache the nodal indices at which the displacement components are stored
293  unsigned u_nodal_index[n_dim];
294  for (unsigned i = 0; i < n_dim; i++)
295  {
296  u_nodal_index[i] = this->U_index_linear_elasticity_traction[i];
297  }
298 
299  // Integer to hold the local equation number
300  int local_eqn = 0;
301 
302  // Set up memory for the shape functions
303  // Note that in this case, the number of lagrangian coordinates is always
304  // equal to the dimension of the nodes
305  Shape psi(n_node);
306  DShape dpsids(n_node, n_dim - 1);
307 
308  // Set the value of n_intpt
309  unsigned n_intpt = integral_pt()->nweight();
310 
311  // Loop over the integration points
312  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
313  {
314  // Get the integral weight
315  double w = integral_pt()->weight(ipt);
316 
317  // Only need to call the local derivatives
318  dshape_local_at_knot(ipt, psi, dpsids);
319 
320  // Calculate the Eulerian and Lagrangian coordinates
321  Vector<double> interpolated_x(n_dim, 0.0);
322 
323  // Also calculate the surface Vectors (derivatives wrt local coordinates)
324  DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
325 
326  // Calculate displacements and derivatives
327  for (unsigned l = 0; l < n_node; l++)
328  {
329  // Loop over directions
330  for (unsigned i = 0; i < n_dim; i++)
331  {
332  // Calculate the Eulerian coords
333  const double x_local = nodal_position(l, i);
334  interpolated_x[i] += x_local * psi(l);
335 
336  // Loop over LOCAL derivative directions, to calculate the tangent(s)
337  for (unsigned j = 0; j < n_dim - 1; j++)
338  {
339  interpolated_A(j, i) += x_local * dpsids(l, j);
340  }
341  }
342  }
343 
344  // Now find the local metric tensor from the tangent Vectors
345  DenseMatrix<double> A(n_dim - 1);
346  for (unsigned i = 0; i < n_dim - 1; i++)
347  {
348  for (unsigned j = 0; j < n_dim - 1; j++)
349  {
350  // Initialise surface metric tensor to zero
351  A(i, j) = 0.0;
352 
353  // Take the dot product
354  for (unsigned k = 0; k < n_dim; k++)
355  {
356  A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
357  }
358  }
359  }
360 
361  // Get the outer unit normal
362  Vector<double> interpolated_normal(n_dim);
363  outer_unit_normal(ipt, interpolated_normal);
364 
365  // Find the determinant of the metric tensor
366  double Adet = 0.0;
367  switch (n_dim)
368  {
369  case 2:
370  Adet = A(0, 0);
371  break;
372  case 3:
373  Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
374  break;
375  default:
376  throw OomphLibError(
377  "Wrong dimension in LinearElasticityTractionElement",
378  "LinearElasticityTractionElement::fill_in_contribution_to_"
379  "residuals()",
380  OOMPH_EXCEPTION_LOCATION);
381  }
382 
383  // Premultiply the weights and the square-root of the determinant of
384  // the metric tensor
385  double W = w * sqrt(Adet);
386 
387  // Now calculate the load
388  Vector<double> traction(n_dim);
389  get_traction(time, ipt, interpolated_x, interpolated_normal, traction);
390 
391  // Loop over the test functions, nodes of the element
392  for (unsigned l = 0; l < n_node; l++)
393  {
394  // Loop over the displacement components
395  for (unsigned i = 0; i < n_dim; i++)
396  {
397  local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
398  /*IF it's not a boundary condition*/
399  if (local_eqn >= 0)
400  {
401  // Add the loading terms to the residuals
402  residuals[local_eqn] -= traction[i] * psi(l) * W;
403  }
404  } // End of if not boundary condition
405  } // End of loop over shape functions
406  } // End of loop over integration points
407  }
408 
409 
410 } // namespace oomph
411 
412 #endif
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
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4342
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4630
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary,...
Definition: elements.h:4501
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5002
A general Finite Element class.
Definition: elements.h:1317
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3054
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
Definition: elements.cc:5162
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2488
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2474
A class for elements that allow the imposition of an applied traction in the equations of linear elas...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void output(std::ostream &outfile)
Output function.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Specify the value of nodal zeta from the face geometry The "global" intrinsic coordinate of the eleme...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
void fill_in_contribution_to_residuals_linear_elasticity_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function. Arguments: Eulerian coordinate; outer unit normal; applied t...
void output(FILE *file_pt)
C_style output function.
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction function pointer.
LinearElasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
virtual void get_traction(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Get the traction vector: Pass number of integration point (dummy), Eulerlian coordinate and normal ve...
Vector< unsigned > U_index_linear_elasticity_traction
Index at which the i-th displacement component is stored.
void traction(const double &time, const Vector< double > &s, Vector< double > &traction)
Compute traction vector at specified local coordinate Should only be used for post-processing; ignore...
An OomphLibError object which should be thrown when an run-time error is encountered....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
void Zero_traction_fct(const double &time, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...