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// 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 /// Namespace containing the zero traction function for linear elasticity
51 /// traction elements
52 //=======================================================================
53 namespace TimeHarmonicLinearElasticityTractionElementHelper
54  {
55 
56  //=======================================================================
57  /// Default load function (zero traction)
58  //=======================================================================
60  const Vector<double>& N,
61  Vector<std::complex<double> >& load)
62  {
63  unsigned n_dim=load.size();
64  for (unsigned i=0;i<n_dim;i++) {load[i]=std::complex<double>(0.0,0.0);}
65  }
66 
67  }
68 
69 
70 //======================================================================
71 /// A class for elements that allow the imposition of an applied traction
72 /// in the equations of time-harmonic linear elasticity.
73 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
74 /// class and and thus, we can be generic enough without the need to have
75 /// a separate equations class.
76 //======================================================================
77 template <class ELEMENT>
79 public virtual FaceGeometry<ELEMENT>,
80  public virtual FaceElement
81  {
82  protected:
83 
84  /// 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 Vector<double> &x,
94  const Vector<double> &n,
95  Vector<std::complex<double> > &result);
96 
97 
98  /// \short Get the traction vector: Pass number of integration point (dummy),
99  /// Eulerian coordinate and normal vector and return the load vector
100  /// (not all of the input arguments will be
101  /// required for all specific load functions but the list should
102  /// cover all cases). This function is virtual so it can be
103  /// overloaded for FSI.
104  virtual void get_traction(const unsigned& intpt,
105  const Vector<double>& x,
106  const Vector<double>& n,
107  Vector<std::complex<double> >& traction)
108  {
110  }
111 
112 
113  /// \short Helper function that actually calculates the residuals
114  // This small level of indirection is required to avoid calling
115  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
116  // which causes all kinds of pain if overloading later on
118  Vector<double> &residuals);
119 
120 
121  public:
122 
123  /// \short Constructor, which takes a "bulk" element and the
124  /// value of the index and its limit
126  const int &face_index) :
127  FaceGeometry<ELEMENT>(), FaceElement()
128  {
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);
140  for(unsigned i=0;i<n_dim;i++)
141  {
143  cast_element_pt->u_index_time_harmonic_linear_elasticity(i);
144  }
145 
146  // Zero traction
149 
150 #ifdef PARANOID
151  {
152  //Check that the element is not a refineable 3d element
153  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
154  //If it's three-d
155  if(elem_pt->dim()==3)
156  {
157  //Is it refineable
158  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
159  if(ref_el_pt!=0)
160  {
161  if (this->has_hanging_nodes())
162  {
163  throw OomphLibError(
164  "This flux element will not work correctly if nodes are hanging\n",
165  OOMPH_CURRENT_FUNCTION,
166  OOMPH_EXCEPTION_LOCATION);
167  }
168  }
169  }
170  }
171 #endif
172  }
173 
174 
175  /// Reference to the traction function pointer
176  void (* &traction_fct_pt())(const Vector<double>& x,
177  const Vector<double>& n,
179  {return Traction_fct_pt;}
180 
181 
182  /// Return the residuals
184  {
186  (residuals);
187  }
188 
189 
190 
191  /// Fill in contribution from Jacobian
193  DenseMatrix<double> &jacobian)
194  {
195  //Call the residuals
197  (residuals);
198  }
199 
200  /// Specify the value of nodal zeta from the face geometry
201  /// \short The "global" intrinsic coordinate of the element when
202  /// viewed as part of a geometric object should be given by
203  /// the FaceElement representation, by default (needed to break
204  /// indeterminacy if bulk element is SolidElement)
205  double zeta_nodal(const unsigned &n, const unsigned &k,
206  const unsigned &i) const
207  {return FaceElement::zeta_nodal(n,k,i);}
208 
209  /// \short Output function
210  void output(std::ostream &outfile)
211  {
212  unsigned nplot=5;
213  output(outfile,nplot);
214  }
215 
216  /// \short Output function
217  void output(std::ostream &outfile, const unsigned &nplot)
218  {
219  unsigned ndim=dim();
220  Vector<double> s(ndim);
221  Vector<double> x(ndim+1);
223 
224  // Tecplot header info
225  outfile << this->tecplot_zone_string(nplot);
226 
227  // Loop over plot points
228  unsigned num_plot_points=this->nplot_points(nplot);
229  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
230  {
231  // Get local coordinates of plot point
232  this->get_s_plot(iplot,nplot,s);
233 
234  // Get Eulerian coordinates and displacements
235  this->interpolated_x(s,x);
236  this->traction(s,traction);
237 
238  //Output the x,y,..
239  for(unsigned i=0;i<ndim+1;i++)
240  {outfile << x[i] << " ";}
241 
242  // Output u,v,..
243  for(unsigned i=0;i<ndim+1;i++)
244  {outfile << traction[i].real() << " ";}
245 
246  // Output u,v,..
247  for(unsigned i=0;i<ndim+1;i++)
248  {outfile << traction[i].imag() << " ";}
249 
250  outfile << std::endl;
251  }
252 
253  // Write tecplot footer (e.g. FE connectivity lists)
254  this->write_tecplot_zone_footer(outfile,nplot);
255  }
256 
257  /// \short C_style output function
258  void output(FILE* file_pt)
259  {FiniteElement::output(file_pt);}
260 
261  /// \short C-style output function
262  void output(FILE* file_pt, const unsigned &n_plot)
263  {FiniteElement::output(file_pt,n_plot);}
264 
265 
266  /// \short Compute traction vector at specified local coordinate
267  /// Should only be used for post-processing; ignores dependence
268  /// on integration point!
269  void traction(const Vector<double>& s,
270  Vector<std::complex<double> >& traction);
271 
272  };
273 
274 ///////////////////////////////////////////////////////////////////////
275 ///////////////////////////////////////////////////////////////////////
276 ///////////////////////////////////////////////////////////////////////
277 
278 //=====================================================================
279 /// Compute traction vector at specified local coordinate
280 /// Should only be used for post-processing; ignores dependence
281 /// on integration point!
282 //=====================================================================
283 template<class ELEMENT>
285  const Vector<double>& s, Vector<std::complex<double> >& traction)
286  {
287  unsigned n_dim = this->nodal_dimension();
288 
289  // Position vector
290  Vector<double> x(n_dim);
291  interpolated_x(s,x);
292 
293  // Outer unit normal
294  Vector<double> unit_normal(n_dim);
295  outer_unit_normal(s,unit_normal);
296 
297  // Dummy
298  unsigned ipt=0;
299 
300  // Traction vector
301  get_traction(ipt,x,unit_normal,traction);
302 
303  }
304 
305 
306 //=====================================================================
307 /// Return the residuals for the
308 /// TimeHarmonicLinearElasticityTractionElement equations
309 //=====================================================================
310 template<class ELEMENT>
313  Vector<double> &residuals)
314  {
315 
316  //Find out how many nodes there are
317  unsigned n_node = nnode();
318 
319  #ifdef PARANOID
320  //Find out how many positional dofs there are
321  unsigned n_position_type = this->nnodal_position_type();
322  if(n_position_type != 1)
323  {
324  throw OomphLibError(
325  "TimeHarmonicLinearElasticity is not yet implemented for more than one position type",
326  OOMPH_CURRENT_FUNCTION,
327  OOMPH_EXCEPTION_LOCATION);
328  }
329 #endif
330 
331  //Find out the dimension of the node
332  const unsigned n_dim = this->nodal_dimension();
333 
334  //Cache the nodal indices at which the displacement components are stored
335  std::vector<std::complex<unsigned> > u_nodal_index(n_dim);
336  for(unsigned i=0;i<n_dim;i++)
337  {
338  //u_nodal_index[i].real() =
339  // this->U_index_time_harmonic_linear_elasticity_traction[i].real();
340  //
341  //u_nodal_index[i].imag() =
342  // this->U_index_time_harmonic_linear_elasticity_traction[i].imag();
343 
344  u_nodal_index[i] =
345  this->U_index_time_harmonic_linear_elasticity_traction[i];
346  }
347 
348  //Integer to hold the local equation number
349  int local_eqn=0;
350 
351  //Set up memory for the shape functions
352  //Note that in this case, the number of lagrangian coordinates is always
353  //equal to the dimension of the nodes
354  Shape psi(n_node);
355  DShape dpsids(n_node,n_dim-1);
356 
357  //Set the value of n_intpt
358  unsigned n_intpt = integral_pt()->nweight();
359 
360  //Loop over the integration points
361  for(unsigned ipt=0;ipt<n_intpt;ipt++)
362  {
363  //Get the integral weight
364  double w = integral_pt()->weight(ipt);
365 
366  //Only need to call the local derivatives
367  dshape_local_at_knot(ipt,psi,dpsids);
368 
369  //Calculate the Eulerian and Lagrangian coordinates
370  Vector<double> interpolated_x(n_dim,0.0);
371 
372  //Also calculate the surface Vectors (derivatives wrt local coordinates)
373  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
374 
375  //Calculate displacements and derivatives
376  for(unsigned l=0;l<n_node;l++)
377  {
378  //Loop over directions
379  for(unsigned i=0;i<n_dim;i++)
380  {
381  //Calculate the Eulerian coords
382  const double x_local = nodal_position(l,i);
383  interpolated_x[i] += x_local*psi(l);
384 
385  //Loop over LOCAL derivative directions, to calculate the tangent(s)
386  for(unsigned j=0;j<n_dim-1;j++)
387  {
388  interpolated_A(j,i) += x_local*dpsids(l,j);
389  }
390  }
391  }
392 
393  //Now find the local metric tensor from the tangent Vectors
394  DenseMatrix<double> A(n_dim-1);
395  for(unsigned i=0;i<n_dim-1;i++)
396  {
397  for(unsigned j=0;j<n_dim-1;j++)
398  {
399  //Initialise surface metric tensor to zero
400  A(i,j) = 0.0;
401 
402  //Take the dot product
403  for(unsigned k=0;k<n_dim;k++)
404  {
405  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
406  }
407  }
408  }
409 
410  //Get the outer unit normal
411  Vector<double> interpolated_normal(n_dim);
412  outer_unit_normal(ipt,interpolated_normal);
413 
414  //Find the determinant of the metric tensor
415  double Adet =0.0;
416  switch(n_dim)
417  {
418  case 2:
419  Adet = A(0,0);
420  break;
421  case 3:
422  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
423  break;
424  default:
425  throw
427  "Wrong dimension in TimeHarmonicLinearElasticityTractionElement",
428  "TimeHarmonicLinearElasticityTractionElement::fill_in_contribution_to_residuals()",
429  OOMPH_EXCEPTION_LOCATION);
430  }
431 
432  //Premultiply the weights and the square-root of the determinant of
433  //the metric tensor
434  double W = w*sqrt(Adet);
435 
436  //Now calculate the load
437  Vector<std::complex<double> > traction(n_dim);
438  get_traction(ipt,
439  interpolated_x,
440  interpolated_normal,
441  traction);
442 
443  //Loop over the test functions, nodes of the element
444  for(unsigned l=0;l<n_node;l++)
445  {
446  //Loop over the displacement components
447  for(unsigned i=0;i<n_dim;i++)
448  {
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  }
469  } //End of loop over shape functions
470  } //End of loop over integration points
471 
472  }
473 
474 
475 
476 
477 
478 }
479 
480 #endif
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...
Vector< std::complex< unsigned > > U_index_time_harmonic_linear_elasticity_traction
Index at which the i-th displacement component is stored.
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
unsigned ndim() const
Access function to # of Eulerian coordinates.
Definition: geom_objects.h:185
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
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:2962
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
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...
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:2976
A general Finite Element class.
Definition: elements.h:1271
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...
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
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(*&)(const Vector< double > &x, const Vector< double > &n, Vector< std::complex< double > > &traction) traction_fct_pt()
Reference to the traction function pointer.
void output(std::ostream &outfile, const unsigned &nplot)
Output function.
static char t char * s
Definition: cfortran.h:572
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
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
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:4277
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
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:2949
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 fill_in_contribution_to_residuals_time_harmonic_linear_elasticity_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
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 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)