pml_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-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 // Header file for elements that are used to apply surface loads to
27 // the equations of linear elasticity
28 
29 #ifndef OOMPH_PML_TIME_HARMONIC_LINEAR_ELASTICITY_TRACTION_ELEMENTS_HEADER
30 #define OOMPH_PML_TIME_HARMONIC_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 PMLTimeHarmonicLinearElasticityTractionElementHelper
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 PMLTimeHarmonicLinearElasticityTractionElementHelper
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)
88  void (*Traction_fct_pt)(const Vector<double>& x,
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  {
103  Traction_fct_pt(x, n, traction);
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  /// PMLTimeHarmonicLinearElasticityTractionElement 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("PMLTimeHarmonicLinearElasticity is not yet "
329  "implemented 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  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 PMLTimeHarmonicLinearElasticityTractionElement",
431  "PMLTimeHarmonicLinearElasticityTractionElement::fill_in_"
432  "contribution_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....
A class for elements that allow the imposition of an applied traction in the equations of time-harmon...
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...
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 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...
PMLTimeHarmonicLinearElasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
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.
Vector< std::complex< unsigned > > U_index_time_harmonic_linear_elasticity_traction
Index at which the i-th displacement component is stored.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
void output(std::ostream &outfile, const unsigned &nplot)
Output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
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...
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 Vector< double > &x, const Vector< double > &N, Vector< std::complex< double >> &load)
Default load function (zero traction)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...