axisym_solid_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 elasticity
32 
33 #ifndef OOMPH_AXISYMM_SOLID_TRACTION_ELEMENTS_HEADER
34 #define OOMPH_AXISYMM_SOLID_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 // OOMPH-LIB headers
42 #include "../generic/Qelements.h"
43 #include "../generic/hermite_elements.h"
44 
45 namespace oomph
46 {
47 
48 //======================================================================
49 /// A class for elements that allow the imposition of an applied traction
50 /// in the principle of virtual displacements.
51 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
52 /// class and and thus, we can be generic enough without the need to have
53 /// a separate equations class.
54 //======================================================================
55 template <class ELEMENT>
57 public virtual FaceGeometry<ELEMENT>, public virtual FaceElement
58 {
59 private:
60 
61  /// Pointer to an imposed traction function
62  void (*Traction_fct_pt)(const double& time, const Vector<double> &x,
63  const Vector<double> &n,
64  Vector<double> &result);
65 
66 protected:
67 
68  /// Return the surface traction force
69  void get_traction(const double& time, const Vector<double> &x,
70  const Vector<double> &n,Vector<double> &result) const
71  {
72  //If the function pointer is zero return zero
73  if(Traction_fct_pt == 0)
74  {
75  //Loop over dimensions and set body forces to zero
76  //It's axisymmetric, so we only have "two" dimensions
77  for(unsigned i=0;i<2;i++) {result[i] = 0.0;}
78  }
79  //Otherwise call the function
80  else
81  {
82  (*Traction_fct_pt)(time,x,n,result);
83  }
84  }
85 
86 public:
87 
88  /// \short Constructor, which takes a "bulk" element and
89  /// the value of the index and its limit
91  const int &face_index) :
92  FaceGeometry<ELEMENT>(), FaceElement()
93  {
94  //Attach the geometrical information to the element. N.B. This function
95  //also assigns nbulk_value from the required_nvalue of the bulk element
96  element_pt->build_face_element(face_index,this);
97 
98  //Set the body force function pointer to zero
99  Traction_fct_pt = 0;
100  }
101 
102  /// Return the imposed traction pointer
103  void (* &traction_fct_pt())(const double&, const Vector<double> &,
104  const Vector<double> &, Vector<double> &)
105  {return Traction_fct_pt;}
106 
107  /// Return the residuals
109 
110 /// Return the jacobian
112  DenseMatrix<double> &jacobian)
113  {
115  //Call the generic FD jacobian calculation
117  }
118 
119  /// Overload the output function
120  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
121 
122  /// Output function: x,y,[z],u,v,[w],p in tecplot format
123  void output(std::ostream &outfile, const unsigned &n_plot)
124  {FiniteElement::output(outfile,n_plot);}
125 
126  /// Overload the output function
127  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
128 
129  /// Output function: x,y,[z],u,v,[w],p in tecplot format
130  void output(FILE* file_pt, const unsigned &n_plot)
131  {FiniteElement::output(file_pt,n_plot);}
132 
133 };
134 
135 
136 
137 /////////////////////////////////////////////////////////////////////////
138 /////////////////////////////////////////////////////////////////////////
139 /////////////////////////////////////////////////////////////////////////
140 
141 
142 //=======================================================================
143 /// Return the residuals for the AxisymmetricSolidTractionElements
144 //=======================================================================
145 template<class ELEMENT>
148 {
149  //Find out how many nodes there are
150  unsigned n_node = nnode();
151  //Find out how many positional dofs there are
152  unsigned n_position_type = nnodal_position_type();
153 
154  //Integer to hold the local equation number
155  int local_eqn=0;
156 
157  //Set up memory for the shape functions
158  //The surface is 1D, so we only have one local derivative
159  Shape psi(n_node,n_position_type);
160  DShape dpsids(n_node,n_position_type,1);
161 
162  //Set the value of n_intpt
163  unsigned n_intpt = integral_pt()->nweight();
164 
165  //Loop over the integration points
166  for(unsigned ipt=0;ipt<n_intpt;ipt++)
167  {
168  //Get the integral weight
169  double w = integral_pt()->weight(ipt);
170 
171  //Only need to call the local derivatives
172  dshape_local_at_knot(ipt,psi,dpsids);
173 
174  //Calculate the global position and lagrangian coordinate
175  Vector<double> interpolated_x(2,0.0), interpolated_xi(2,0.0);
176  //Calculate the global and lagrangian derivtives wrt the local coordinates
177  Vector<double> interpolated_dxds(2,0.0), interpolated_dxids(2,0.0);
178 
179  //Calculate displacements and derivatives
180  for(unsigned l=0;l<n_node;l++)
181  {
182  //Loop over positional dofs
183  for(unsigned k=0;k<n_position_type;k++)
184  {
185  //Loop over the number of lagrangian coordinates (2)
186  for(unsigned i=0;i<2;i++)
187  {
188  //Calculate the global position
189  interpolated_x[i] +=
190  nodal_position_gen(l,bulk_position_type(k),i)*psi(l,k);
191  interpolated_xi[i] +=
192  this->lagrangian_position_gen(l,bulk_position_type(k),i)*psi(l,k);
193  //Calculate the derivatives of the global and lagrangian coordinates
194  interpolated_dxds[i] +=
195  nodal_position_gen(l,bulk_position_type(k),i)*dpsids(l,k,0);
196  interpolated_dxids[i] +=
197  this->lagrangian_position_gen(l,bulk_position_type(k),i)
198  *dpsids(l,k,0);
199  }
200  }
201  }
202 
203  //Now calculate the entries of the deformed surface metric tensor
204  //Now find the local deformed metric tensor from the tangent Vectors
205  DenseMatrix<double> A(2);
206  //The off-diagonal terms are Zero
207  A(0,1) = A(1,0) = 0.0;
208  //The diagonal terms are a little complicated
209  A(0,0) =
210  (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1])*
211  (interpolated_dxds[0] - interpolated_x[1]*interpolated_dxids[1]) +
212  (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1])*
213  (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);
214 
215 
216  A(1,1) = (interpolated_x[0]*sin(interpolated_xi[1]) +
217  interpolated_x[1]*cos(interpolated_xi[1]))*
218  (interpolated_x[0]*sin(interpolated_xi[1]) +
219  interpolated_x[1]*cos(interpolated_xi[1]));
220 
221  //Premultiply the weights and the square-root of the determinant of
222  //the metric tensor
223  double W = w*sqrt(A(0,0)*A(1,1));
224 
225  //Also find the normal -- just the cross product of the metric tensors
226  //but I want to express it in terms of e_r and e_theta components
227  //N.B. There is an issue at theta = 0,pi, where the normal is e_{r},
228  //but given that I never assemble it, should be OK!
229  //The minus sign is chosen to ensure that the normal is really outward
230  Vector<double> interpolated_normal(2);
231  //Component in the e_{r} direction
232  interpolated_normal[0] = -1.0*
233  (interpolated_x[0]*sin(interpolated_xi[1]) +
234  interpolated_x[1]*cos(interpolated_xi[1]))*
235  (interpolated_dxds[1] + interpolated_x[0]*interpolated_dxids[1]);
236  //Component in the e_{theta} direction
237  interpolated_normal[1] = -1.0*
238  (interpolated_x[0]*sin(interpolated_xi[1]) +
239  interpolated_x[1]*cos(interpolated_xi[1]))*
240  (interpolated_x[1]*interpolated_dxids[1] - interpolated_dxds[0]);
241 
242  //If we're on the north or south face need to flip normal
243  if(s_fixed_value()==-1)
244  {
245  interpolated_normal[0] *= -1.0;
246  interpolated_normal[1] *= -1.0;
247  }
248 
249  //Now adjust and scale the normal
250  double length = 0.0;
251  for(unsigned i=0;i<2;i++)
252  {
253  interpolated_normal[i] *= normal_sign();
254  length += interpolated_normal[i]*interpolated_normal[i];
255  }
256  for(unsigned i=0;i<2;i++)
257  {
258  interpolated_normal[i] /= sqrt(length);
259  }
260 
261  //Now calculate the load
262  Vector<double> traction(2);
263 
264  //Normal is outwards
265  get_traction(time(),interpolated_x,interpolated_normal,traction);
266 
267  //=====LOAD TERMS FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
268 
269  //Loop over the test functions, nodes of the element
270  for(unsigned l=0;l<n_node;l++)
271  {
272  //Loop of types of dofs
273  for(unsigned k=0;k<n_position_type;k++)
274  {
275  //Loop over the displacement components
276  for(unsigned i=0;i<2;i++)
277  {
278  local_eqn =
279  this->position_local_eqn(l,bulk_position_type(k),i);
280  /*IF it's not a boundary condition*/
281  if(local_eqn >= 0)
282  {
283  //Add the loading terms to the residuals
284  residuals[local_eqn] -= traction[i]*psi(l,k)*W;
285  }
286  }
287  } //End of if not boundary condition
288  } //End of loop over shape functions
289  } //End of loop over integration points
290 
291 }
292 
293 
294 }
295 
296 #endif
void(*&)(const double &, const Vector< double > &, const Vector< double > &, Vector< double > &) traction_fct_pt()
Return the imposed traction pointer.
AxisymmetricSolidTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
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
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
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
A general Finite Element class.
Definition: elements.h:1271
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Return the jacobian.
void get_traction(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result) const
Return the surface traction force.
void output(std::ostream &outfile)
Overload the output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
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 output(FILE *file_pt, const unsigned &n_plot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
void output(FILE *file_pt)
Overload the output function.