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