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