time_harmonic_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 TimeHarmonicLinearElasticityTractionElementHelper
48 {
49 //=======================================================================
50 /// Default load function (zero traction)
51 //=======================================================================
53 const Vector<double>& N,
54 Vector<std::complex<double>>& load)
55 {
56 unsigned n_dim = load.size();
57 for (unsigned i = 0; i < n_dim; i++)
58 {
59 load[i] = std::complex<double>(0.0, 0.0);
60 }
61 }
62
63 } // namespace TimeHarmonicLinearElasticityTractionElementHelper
64
65
66 //======================================================================
67 /// A class for elements that allow the imposition of an applied traction
68 /// in the equations of time-harmonic linear elasticity.
69 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
70 /// class and and thus, we can be generic enough without the need to have
71 /// a separate equations class.
72 //======================================================================
73 template<class ELEMENT>
75 : public virtual FaceGeometry<ELEMENT>,
76 public virtual FaceElement
77 {
78 protected:
79 /// Index at which the i-th displacement component is stored
82
83 /// Pointer to an imposed traction function. Arguments:
84 /// Eulerian coordinate; outer unit normal;
85 /// applied traction. (Not all of the input arguments will be
86 /// required for all specific load functions but the list should
87 /// cover all cases)
89 const Vector<double>& n,
91
92
93 /// Get the traction vector: Pass number of integration point
94 /// (dummy), Eulerian 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 unsigned& intpt,
99 const Vector<double>& x,
100 const Vector<double>& n,
101 Vector<std::complex<double>>& traction)
102 {
104 }
105
106
107 /// Helper function that actually calculates the residuals
108 // This small level of indirection is required to avoid calling
109 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
110 // which causes all kinds of pain if overloading later on
112 Vector<double>& residuals);
113
114
115 public:
116 /// Constructor, which takes a "bulk" element and the
117 /// value of the index and its limit
119 FiniteElement* const& element_pt, const int& face_index)
120 : FaceGeometry<ELEMENT>(), FaceElement()
121 {
122 // Attach the geometrical information to the element. N.B. This function
123 // also assigns nbulk_value from the required_nvalue of the bulk element
124 element_pt->build_face_element(face_index, this);
125
126 // Find the dimension of the problem
127 unsigned n_dim = element_pt->nodal_dimension();
128
129 // Find the index at which the displacement unknowns are stored
130 ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
132 for (unsigned i = 0; i < n_dim; i++)
133 {
135 cast_element_pt->u_index_time_harmonic_linear_elasticity(i);
136 }
137
138 // Zero traction
141
142#ifdef PARANOID
143 {
144 // Check that the element is not a refineable 3d element
145 ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
146 // If it's three-d
147 if (elem_pt->dim() == 3)
148 {
149 // Is it refineable
150 RefineableElement* ref_el_pt =
151 dynamic_cast<RefineableElement*>(elem_pt);
152 if (ref_el_pt != 0)
153 {
154 if (this->has_hanging_nodes())
155 {
156 throw OomphLibError("This flux element will not work correctly "
157 "if nodes are hanging\n",
158 OOMPH_CURRENT_FUNCTION,
159 OOMPH_EXCEPTION_LOCATION);
160 }
161 }
162 }
163 }
164#endif
165 }
166
167
168 /// Reference to the traction function pointer
169 void (*&traction_fct_pt())(const Vector<double>& x,
170 const Vector<double>& n,
172 {
173 return Traction_fct_pt;
174 }
175
176
177 /// Return the residuals
179 {
181 residuals);
182 }
183
184
185 /// Fill in contribution from Jacobian
187 DenseMatrix<double>& jacobian)
188 {
189 // Call the residuals
191 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 unsigned nplot = 5;
210 output(outfile, nplot);
211 }
212
213 /// Output function
214 void output(std::ostream& outfile, const unsigned& nplot)
215 {
216 unsigned ndim = dim();
218 Vector<double> x(ndim + 1);
220
221 // Tecplot header info
222 outfile << this->tecplot_zone_string(nplot);
223
224 // Loop over plot points
225 unsigned num_plot_points = this->nplot_points(nplot);
226 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
227 {
228 // Get local coordinates of plot point
229 this->get_s_plot(iplot, nplot, s);
230
231 // Get Eulerian coordinates and displacements
232 this->interpolated_x(s, x);
233 this->traction(s, traction);
234
235 // Output the x,y,..
236 for (unsigned i = 0; i < ndim + 1; i++)
237 {
238 outfile << x[i] << " ";
239 }
240
241 // Output u,v,..
242 for (unsigned i = 0; i < ndim + 1; i++)
243 {
244 outfile << traction[i].real() << " ";
245 }
246
247 // Output u,v,..
248 for (unsigned i = 0; i < ndim + 1; i++)
249 {
250 outfile << traction[i].imag() << " ";
251 }
252
253 outfile << std::endl;
254 }
255
256 // Write tecplot footer (e.g. FE connectivity lists)
257 this->write_tecplot_zone_footer(outfile, nplot);
258 }
259
260 /// C_style output function
261 void output(FILE* file_pt)
262 {
263 FiniteElement::output(file_pt);
264 }
265
266 /// C-style output function
267 void output(FILE* file_pt, const unsigned& n_plot)
268 {
269 FiniteElement::output(file_pt, n_plot);
270 }
271
272
273 /// Compute traction vector at specified local coordinate
274 /// Should only be used for post-processing; ignores dependence
275 /// on integration point!
276 void traction(const Vector<double>& s,
277 Vector<std::complex<double>>& traction);
278 };
279
280 /// ////////////////////////////////////////////////////////////////////
281 /// ////////////////////////////////////////////////////////////////////
282 /// ////////////////////////////////////////////////////////////////////
283
284 //=====================================================================
285 /// Compute traction vector at specified local coordinate
286 /// Should only be used for post-processing; ignores dependence
287 /// on integration point!
288 //=====================================================================
289 template<class ELEMENT>
291 const Vector<double>& s, Vector<std::complex<double>>& traction)
292 {
293 unsigned n_dim = this->nodal_dimension();
294
295 // Position vector
296 Vector<double> x(n_dim);
297 interpolated_x(s, x);
298
299 // Outer unit normal
300 Vector<double> unit_normal(n_dim);
301 outer_unit_normal(s, unit_normal);
302
303 // Dummy
304 unsigned ipt = 0;
305
306 // Traction vector
307 get_traction(ipt, x, unit_normal, traction);
308 }
309
310
311 //=====================================================================
312 /// Return the residuals for the
313 /// TimeHarmonicLinearElasticityTractionElement equations
314 //=====================================================================
315 template<class ELEMENT>
318 Vector<double>& residuals)
319 {
320 // Find out how many nodes there are
321 unsigned n_node = nnode();
322
323#ifdef PARANOID
324 // Find out how many positional dofs there are
325 unsigned n_position_type = this->nnodal_position_type();
326 if (n_position_type != 1)
327 {
328 throw OomphLibError("TimeHarmonicLinearElasticity is not yet implemented "
329 "for more than one position type",
330 OOMPH_CURRENT_FUNCTION,
331 OOMPH_EXCEPTION_LOCATION);
332 }
333#endif
334
335 // Find out the dimension of the node
336 const unsigned n_dim = this->nodal_dimension();
337
338 // Cache the nodal indices at which the displacement components are stored
339 std::vector<std::complex<unsigned>> u_nodal_index(n_dim);
340 for (unsigned i = 0; i < n_dim; i++)
341 {
342 // u_nodal_index[i].real() =
343 // this->U_index_time_harmonic_linear_elasticity_traction[i].real();
344 //
345 // u_nodal_index[i].imag() =
346 // this->U_index_time_harmonic_linear_elasticity_traction[i].imag();
347
348 u_nodal_index[i] =
349 this->U_index_time_harmonic_linear_elasticity_traction[i];
350 }
351
352 // Integer to hold the local equation number
353 int local_eqn = 0;
354
355 // Set up memory for the shape functions
356 // Note that in this case, the number of lagrangian coordinates is always
357 // equal to the dimension of the nodes
358 Shape psi(n_node);
359 DShape dpsids(n_node, n_dim - 1);
360
361 // Set the value of n_intpt
362 unsigned n_intpt = integral_pt()->nweight();
363
364 // Loop over the integration points
365 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
366 {
367 // Get the integral weight
368 double w = integral_pt()->weight(ipt);
369
370 // Only need to call the local derivatives
371 dshape_local_at_knot(ipt, psi, dpsids);
372
373 // Calculate the Eulerian and Lagrangian coordinates
374 Vector<double> interpolated_x(n_dim, 0.0);
375
376 // Also calculate the surface Vectors (derivatives wrt local coordinates)
377 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
378
379 // Calculate displacements and derivatives
380 for (unsigned l = 0; l < n_node; l++)
381 {
382 // Loop over directions
383 for (unsigned i = 0; i < n_dim; i++)
384 {
385 // Calculate the Eulerian coords
386 const double x_local = nodal_position(l, i);
387 interpolated_x[i] += x_local * psi(l);
388
389 // Loop over LOCAL derivative directions, to calculate the tangent(s)
390 for (unsigned j = 0; j < n_dim - 1; j++)
391 {
392 interpolated_A(j, i) += x_local * dpsids(l, j);
393 }
394 }
395 }
396
397 // Now find the local metric tensor from the tangent Vectors
398 DenseMatrix<double> A(n_dim - 1);
399 for (unsigned i = 0; i < n_dim - 1; i++)
400 {
401 for (unsigned j = 0; j < n_dim - 1; j++)
402 {
403 // Initialise surface metric tensor to zero
404 A(i, j) = 0.0;
405
406 // Take the dot product
407 for (unsigned k = 0; k < n_dim; k++)
408 {
409 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
410 }
411 }
412 }
413
414 // Get the outer unit normal
415 Vector<double> interpolated_normal(n_dim);
416 outer_unit_normal(ipt, interpolated_normal);
417
418 // Find the determinant of the metric tensor
419 double Adet = 0.0;
420 switch (n_dim)
421 {
422 case 2:
423 Adet = A(0, 0);
424 break;
425 case 3:
426 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
427 break;
428 default:
429 throw OomphLibError(
430 "Wrong dimension in TimeHarmonicLinearElasticityTractionElement",
431 "TimeHarmonicLinearElasticityTractionElement::fill_in_contribution_"
432 "to_residuals()",
433 OOMPH_EXCEPTION_LOCATION);
434 }
435
436 // Premultiply the weights and the square-root of the determinant of
437 // the metric tensor
438 double W = w * sqrt(Adet);
439
440 // Now calculate the load
441 Vector<std::complex<double>> traction(n_dim);
442 get_traction(ipt, interpolated_x, interpolated_normal, traction);
443
444 // Loop over the test functions, nodes of the element
445 for (unsigned l = 0; l < n_node; l++)
446 {
447 // Loop over the displacement components
448 for (unsigned i = 0; i < n_dim; i++)
449 {
450 // Real eqn
451 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i].real());
452 /*IF it's not a boundary condition*/
453 if (local_eqn >= 0)
454 {
455 // Add the loading terms to the residuals
456 residuals[local_eqn] -= traction[i].real() * psi(l) * W;
457 }
458
459 // Imag eqn
460 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i].imag());
461 /*IF it's not a boundary condition*/
462 if (local_eqn >= 0)
463 {
464 // Add the loading terms to the residuals
465 residuals[local_eqn] -= traction[i].imag() * psi(l) * W;
466 }
467 }
468 } // End of loop over shape functions
469 } // End of loop over integration points
470 }
471
472
473} // namespace oomph
474
475#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
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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 std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
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
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
unsigned ndim() const
Access function to # of Eulerian coordinates.
Definition: geom_objects.h:177
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
A class for elements that allow the imposition of an applied traction in the equations of time-harmon...
void fill_in_contribution_to_residuals_time_harmonic_linear_elasticity_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void(*&)(const Vector< double > &x, const Vector< double > &n, Vector< std::complex< double > > &traction) traction_fct_pt()
Reference to the traction function pointer.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void traction(const Vector< double > &s, Vector< std::complex< double > > &traction)
Compute traction vector at specified local coordinate Should only be used for post-processing; ignore...
void(* Traction_fct_pt)(const Vector< double > &x, const Vector< double > &n, Vector< std::complex< double > > &result)
Pointer to an imposed traction function. Arguments: Eulerian coordinate; outer unit normal; applied t...
TimeHarmonicLinearElasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void output(std::ostream &outfile, const unsigned &nplot)
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...
Vector< std::complex< unsigned > > U_index_time_harmonic_linear_elasticity_traction
Index at which the i-th displacement component is stored.
virtual void get_traction(const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, Vector< std::complex< double > > &traction)
Get the traction vector: Pass number of integration point (dummy), Eulerian coordinate and normal vec...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void Zero_traction_fct(const Vector< double > &x, const Vector< double > &N, Vector< std::complex< double > > &load)
Default load function (zero traction)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...