time_harmonic_fourier_decomposed_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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for elements that are used to apply surface loads to
31 //the equations of time-harmonic Fourier decomposed linear elasticity
32 
33 #ifndef OOMPH_FOURIER_DECOMPOSED_TIME_HARMONIC_LINEAR_ELASTICITY_TRACTION_ELEMENTS_HEADER
34 #define OOMPH_FOURIER_DECOMPOSED_TIME_HARMONIC_LINEAR_ELASTICITY_TRACTION_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 
42 //OOMPH-LIB headers
43 //#include "../generic/Qelements.h"
44 #include "src/generic/Qelements.h"
45 
46 namespace oomph
47 {
48 
49 
50 //=======================================================================
51 /// Namespace containing the zero traction function for time-harmonic
52 /// Fourier decomposed linear elasticity traction elements
53 //=======================================================================
54 namespace TimeHarmonicFourierDecomposedLinearElasticityTractionElementHelper
55  {
56 
57  //=======================================================================
58  /// Default load function (zero traction)
59  //=======================================================================
61  const Vector<double>& N,
62  Vector<std::complex<double> >& load)
63  {
64  unsigned n_dim=load.size();
65  for (unsigned i=0;i<n_dim;i++) {load[i]=std::complex<double>(0.0,0.0);}
66  }
67 
68  }
69 
70 
71 //======================================================================
72 /// A class for elements that allow the imposition of an applied traction
73 /// in the equations of time-harmonic Fourier decomposed linear elasticity.
74 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
75 /// class and and thus, we can be generic enough without the need to have
76 /// a separate equations class.
77 //======================================================================
78 template <class ELEMENT>
80 public virtual FaceGeometry<ELEMENT>,
81  public virtual FaceElement
82  {
83  protected:
84 
85  /// Index at which the i-th displacement component is stored
88 
89  /// \short Pointer to an imposed traction function. Arguments:
90  /// Eulerian coordinate; outer unit normal;
91  /// applied traction. (Not all of the input arguments will be
92  /// required for all specific load functions but the list should
93  /// cover all cases)
94  void (*Traction_fct_pt)(const Vector<double> &x,
95  const Vector<double> &n,
96  Vector<std::complex<double> > &result);
97 
98 
99  /// \short Get the traction vector: Pass number of integration point (dummy),
100  /// Eulerian coordinate and normal vector and return the load vector
101  /// (not all of the input arguments will be
102  /// required for all specific load functions but the list should
103  /// cover all cases). This function is virtual so it can be
104  /// overloaded for FSI.
105  virtual void get_traction(const unsigned& intpt,
106  const Vector<double>& x,
107  const Vector<double>& n,
108  Vector<std::complex<double> >& traction)
109  {
111  }
112 
113 
114  /// \short Helper function that actually calculates the residuals
115  // This small level of indirection is required to avoid calling
116  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
117  // which causes all kinds of pain if overloading later on
119  Vector<double> &residuals);
120 
121 
122  public:
123 
124  /// \short Constructor, which takes a "bulk" element and the
125  /// value of the index and its limit
127  (FiniteElement* const &element_pt, const int &face_index) :
129  {
130  //Attach the geometrical information to the element. N.B. This function
131  //also assigns nbulk_value from the required_nvalue of the bulk element
132  element_pt->build_face_element(face_index,this);
133 
134  //Find the dimension of the problem
135  unsigned n_dim = element_pt->nodal_dimension();
136 
137  //Find the index at which the displacement unknowns are stored
138  ELEMENT* cast_element_pt = dynamic_cast<ELEMENT*>(element_pt);
139  this->
141  for(unsigned i=0;i<n_dim+1;i++)
142  {
143  //this->
144  // U_index_time_harmonic_fourier_decomposed_linear_elasticity_traction[i].real() =
145  // cast_element_pt->
146  // u_index_time_harmonic_fourier_decomposed_linear_elasticity(i).real();
147  //
148  //this->
149  // U_index_time_harmonic_fourier_decomposed_linear_elasticity_traction[i].imag() =
150  // cast_element_pt->
151  // u_index_time_harmonic_fourier_decomposed_linear_elasticity(i).imag();
152 
153  this->
155  cast_element_pt->
156  u_index_time_harmonic_fourier_decomposed_linear_elasticity(i);
157  }
158 
159  // Zero traction
163  }
164 
165 
166  /// Reference to the traction function pointer
167  void (* &traction_fct_pt())(const Vector<double>& x,
168  const Vector<double>& n,
170  {return Traction_fct_pt;}
171 
172 
173  /// Return the residuals
175  {
177  (residuals);
178  }
179 
180 
181 
182  /// Fill in contribution from Jacobian
184  DenseMatrix<double> &jacobian)
185  {
186  //Call the residuals
188  (residuals);
189  }
190 
191  /// Specify the value of nodal zeta from the face geometry
192  /// \short The "global" intrinsic coordinate of the element when
193  /// viewed as part of a geometric object should be given by
194  /// the FaceElement representation, by default (needed to break
195  /// indeterminacy if bulk element is SolidElement)
196  double zeta_nodal(const unsigned &n, const unsigned &k,
197  const unsigned &i) const
198  {return FaceElement::zeta_nodal(n,k,i);}
199 
200  /// \short Output function
201  void output(std::ostream &outfile)
202  {FiniteElement::output(outfile);}
203 
204  /// \short Output function
205  void output(std::ostream &outfile, const unsigned &n_plot)
206  {FiniteElement::output(outfile,n_plot);}
207 
208  /// \short C_style output function
209  void output(FILE* file_pt)
210  {FiniteElement::output(file_pt);}
211 
212  /// \short C-style output function
213  void output(FILE* file_pt, const unsigned &n_plot)
214  {FiniteElement::output(file_pt,n_plot);}
215 
216 
217  /// \short Compute traction vector at specified local coordinate
218  /// Should only be used for post-processing; ignores dependence
219  /// on integration point!
220  void traction(const Vector<double>& s,
221  Vector<std::complex<double> >& traction);
222 
223  };
224 
225 ///////////////////////////////////////////////////////////////////////
226 ///////////////////////////////////////////////////////////////////////
227 ///////////////////////////////////////////////////////////////////////
228 
229 //=====================================================================
230 /// Compute traction vector at specified local coordinate
231 /// Should only be used for post-processing; ignores dependence
232 /// on integration point!
233 //=====================================================================
234 template<class ELEMENT>
236  traction(const Vector<double>& s, Vector<std::complex<double> >& traction)
237  {
238  unsigned n_dim = this->nodal_dimension();
239 
240  // Position vector
241  Vector<double> x(n_dim);
242  interpolated_x(s,x);
243 
244  // Outer unit normal (only in r and z direction!)
245  Vector<double> unit_normal(n_dim);
246  outer_unit_normal(s,unit_normal);
247 
248  // Dummy
249  unsigned ipt=0;
250 
251  // Traction vector
252  get_traction(ipt,x,unit_normal,traction);
253 
254  }
255 
256 
257 //=====================================================================
258 /// Return the residuals for the
259 /// TimeHarmonicFourierDecomposedLinearElasticityTractionElement equations
260 //=====================================================================
261 template<class ELEMENT>
264  Vector<double> &residuals)
265  {
266 
267  //Find out how many nodes there are
268  unsigned n_node = nnode();
269 
270  #ifdef PARANOID
271  //Find out how many positional dofs there are
272  unsigned n_position_type = this->nnodal_position_type();
273  if(n_position_type != 1)
274  {
275  throw OomphLibError(
276  "TimeHarmonicFourierDecomposedLinearElasticity is not yet implemented for more than one position type",
277  OOMPH_CURRENT_FUNCTION,
278  OOMPH_EXCEPTION_LOCATION);
279  }
280 #endif
281 
282  //Find out the dimension of the node
283  const unsigned n_dim = this->nodal_dimension();
284 
285  //Cache the nodal indices at which the displacement components are stored
286  std::vector<std::complex<unsigned> > u_nodal_index(n_dim+1);
287  for(unsigned i=0;i<n_dim+1;i++)
288  {
289  //u_nodal_index[i].real() =
290  // this->
291  // U_index_time_harmonic_fourier_decomposed_linear_elasticity_traction[i].real();
292  //
293  //u_nodal_index[i].imag() =
294  // this->
295  // U_index_time_harmonic_fourier_decomposed_linear_elasticity_traction[i].imag();
296 
297  u_nodal_index[i] =
298  this->
299  U_index_time_harmonic_fourier_decomposed_linear_elasticity_traction[i];
300  }
301 
302  //Integer to hold the local equation number
303  int local_eqn=0;
304 
305  //Set up memory for the shape functions
306  //Note that in this case, the number of lagrangian coordinates is always
307  //equal to the dimension of the nodes
308  Shape psi(n_node);
309  DShape dpsids(n_node,n_dim-1);
310 
311  //Set the value of n_intpt
312  unsigned n_intpt = integral_pt()->nweight();
313 
314  //Loop over the integration points
315  for(unsigned ipt=0;ipt<n_intpt;ipt++)
316  {
317  //Get the integral weight
318  double w = integral_pt()->weight(ipt);
319 
320  //Only need to call the local derivatives
321  dshape_local_at_knot(ipt,psi,dpsids);
322 
323  //Calculate the Eulerian and Lagrangian coordinates
324  Vector<double> interpolated_x(n_dim,0.0);
325 
326  //Also calculate the surface Vectors (derivatives wrt local coordinates)
327  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
328 
329  //Calculate displacements and derivatives
330  for(unsigned l=0;l<n_node;l++)
331  {
332  //Loop over directions
333  for(unsigned i=0;i<n_dim;i++)
334  {
335  //Calculate the Eulerian coords
336  const double x_local = nodal_position(l,i);
337  interpolated_x[i] += x_local*psi(l);
338 
339  //Loop over LOCAL derivative directions, to calculate the tangent(s)
340  for(unsigned j=0;j<n_dim-1;j++)
341  {
342  interpolated_A(j,i) += x_local*dpsids(l,j);
343  }
344  }
345  }
346 
347  //Now find the local metric tensor from the tangent Vectors
348  DenseMatrix<double> A(n_dim-1);
349  for(unsigned i=0;i<n_dim-1;i++)
350  {
351  for(unsigned j=0;j<n_dim-1;j++)
352  {
353  //Initialise surface metric tensor to zero
354  A(i,j) = 0.0;
355 
356  //Take the dot product
357  for(unsigned k=0;k<n_dim;k++)
358  {
359  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
360  }
361  }
362  }
363 
364  //Get the outer unit normal
365  Vector<double> interpolated_normal(n_dim);
366  outer_unit_normal(ipt,interpolated_normal);
367 
368  //Find the determinant of the metric tensor
369  double Adet =0.0;
370  switch(n_dim)
371  {
372  case 2:
373  Adet = A(0,0);
374  break;
375  case 3:
376  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
377  break;
378  default:
379  throw
381  "Wrong dimension in TimeHarmonicFourierDecomposedLinearElasticityTractionElement",
382  "TimeHarmonicFourierDecomposedLinearElasticityTractionElement::fill_in_contribution_to_residuals()",
383  OOMPH_EXCEPTION_LOCATION);
384  }
385 
386  //Premultiply the weights and the square-root of the determinant of
387  //the metric tensor
388  double W = w*sqrt(Adet);
389 
390  //Now calculate the load
391  Vector<std::complex<double> > traction(n_dim+1);
392  get_traction(ipt,
393  interpolated_x,
394  interpolated_normal,
395  traction);
396 
397  //Loop over the test functions, nodes of the element
398  for(unsigned l=0;l<n_node;l++)
399  {
400  //Loop over the displacement components
401  for(unsigned i=0;i<n_dim+1;i++)
402  {
403 
404  // Real eqn
405  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i].real());
406  /*IF it's not a boundary condition*/
407  if(local_eqn >= 0)
408  {
409  //Add the loading terms to the residuals
410  residuals[local_eqn] -= traction[i].real()*psi(l)*interpolated_x[0]*W;
411  }
412 
413  // Imag eqn
414  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i].imag());
415  /*IF it's not a boundary condition*/
416  if(local_eqn >= 0)
417  {
418  //Add the loading terms to the residuals
419  residuals[local_eqn] -= traction[i].imag()*psi(l)*interpolated_x[0]*W;
420  }
421 
422  }
423  } //End of loop over shape functions
424  } //End of loop over integration points
425 
426  }
427 
428 
429 
430 
431 
432 }
433 
434 #endif
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:2872
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
void Zero_traction_fct(const Vector< double > &x, const Vector< double > &N, Vector< std::complex< double > > &load)
Default load function (zero traction)
A general Finite Element class.
Definition: elements.h:1271
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
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_residuals_time_harmonic_fourier_decomposed_linear_elasticity_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
FaceElement()
Constructor: Initialise all appropriate member data.
Definition: elements.h:4201
TimeHarmonicFourierDecomposedLinearElasticityTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
static char t char * s
Definition: cfortran.h:572
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...
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, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
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_fourier_decomposed_linear_elasticity_traction
Index at which the i-th displacement component is stored.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
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:4922
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...