refineable_linear_elasticity_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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 2016) $
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 refineable linear elasticity elements
31 
32 //Include guard to prevent multiple inclusions of this header
33 #ifndef OOMPH_REFINEABLE_LINEAR_ELASTICITY_ELEMENTS_HEADER
34 #define OOMPH_REFINEABLE_LINEAR_ELASTICITY_ELEMENTS_HEADER
35 
36 //oomph-lib headers
38 #include "../generic/refineable_quad_element.h"
39 #include "../generic/refineable_brick_element.h"
40 #include "../generic/error_estimator.h"
41 
42 namespace oomph
43 {
44 
45 //========================================================================
46 /// Class for Refineable LinearElasticity equations
47 //========================================================================
48  template<unsigned DIM>
50  public virtual LinearElasticityEquations<DIM>,
51  public virtual RefineableElement,
52  public virtual ElementWithZ2ErrorEstimator
53  {
54  public:
55 
56  /// \short Constructor
60  {}
61 
62 
63  /// \short Get the function value u in Vector.
64  /// Note: Given the generality of the interface (this function
65  /// is usually called from black-box documentation or interpolation
66  /// routines), the values Vector sets its own size in here.
67  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
68  Vector<double>& values)
69  {
70  // Create enough initialised storage
71  values.resize(DIM,0.0);
72 
73  //Find out how many nodes there are
74  unsigned n_node = this->nnode();
75 
76  // Shape functions
77  Shape psi(n_node);
78  this->shape(s,psi);
79 
80  //Calculate displacements
81  for(unsigned i=0;i<DIM;i++)
82  {
83  //Get the index at which the i-th velocity is stored
84  unsigned u_nodal_index = this->u_index_linear_elasticity(i);
85  for(unsigned l=0;l<n_node;l++)
86  {
87  values[i]+=this->nodal_value(t,l,u_nodal_index)*psi(l);
88  }
89  }
90  }
91 
92  /// \short Get the current interpolated values (displacements).
93  /// Note: Given the generality of the interface (this function
94  /// is usually called from black-box documentation or interpolation
95  /// routines),the values Vector sets its own size in here.
97  Vector<double> & values)
98  {
99  values.resize(DIM);
100  this->interpolated_u_linear_elasticity(s,values);
101  }
102 
103  /// Number of 'flux' terms for Z2 error estimation
104  unsigned num_Z2_flux_terms()
105  {
106  // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
107  return DIM + DIM*(DIM-1)/2;
108  }
109 
110  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
111  /// in strain tensor.
113  {
114 #ifdef PARANOID
115  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
116  if (flux.size()!=num_entries)
117  {
118  std::ostringstream error_message;
119  error_message << "The flux vector has the wrong number of entries, "
120  << flux.size() << ", whereas it should be "
121  << num_entries << std::endl;
122  throw OomphLibError(error_message.str(),
123  OOMPH_CURRENT_FUNCTION,
124  OOMPH_EXCEPTION_LOCATION);
125  }
126 #endif
127 
128  // Get strain matrix
129  DenseMatrix<double> strain(DIM);
130  this->get_strain(s,strain);
131 
132  // Pack into flux Vector
133  unsigned icount=0;
134 
135  // Start with diagonal terms
136  for(unsigned i=0;i<DIM;i++)
137  {
138  flux[icount]=strain(i,i);
139  icount++;
140  }
141 
142  //Off diagonals row by row
143  for(unsigned i=0;i<DIM;i++)
144  {
145  for(unsigned j=i+1;j<DIM;j++)
146  {
147  flux[icount]=strain(i,j);
148  icount++;
149  }
150  }
151  }
152 
153  /// Number of continuously interpolated values: DIM
154  unsigned ncont_interpolated_values() const {return DIM;}
155 
156  /// Further build function, pass the pointers down to the sons
158  {
159  RefineableLinearElasticityEquations<DIM>* cast_father_element_pt
161  (this->father_element_pt());
162 
163  // Set pointer to body force function
164  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
165 
166  // Set pointer to the contitutive law
167  this->Elasticity_tensor_pt =
168  cast_father_element_pt->elasticity_tensor_pt();
169 
170  // Set the timescale ratio (non-dim. density)
171  this->Lambda_sq_pt = cast_father_element_pt->lambda_sq_pt();
172 
173  /// Set the flag that switches inertia on/off
174  this->Unsteady = cast_father_element_pt->is_inertia_enabled();
175  }
176 
177 
178  private:
179 
180  /// Overloaded helper function to take hanging nodes into account
182  Vector<double> &residuals,DenseMatrix<double> &jacobian,unsigned flag);
183 
184  };
185 
186 
187 //////////////////////////////////////////////////////////////////////
188 //////////////////////////////////////////////////////////////////////
189 //////////////////////////////////////////////////////////////////////
190 
191 //========================================================================
192 /// Class for refineable QLinearElasticityElement elements
193 //========================================================================
194  template<unsigned DIM, unsigned NNODE_1D>
196  public virtual QLinearElasticityElement<DIM,NNODE_1D>,
197  public virtual RefineableLinearElasticityEquations<DIM>,
198  public virtual RefineableQElement<DIM>
199  {
200  public:
201 
202  /// Constructor:
204  QLinearElasticityElement<DIM,NNODE_1D>(),
207  RefineableQElement<DIM>() {}
208 
209  /// Empty rebuild from sons, no need to reconstruct anything here
210  void rebuild_from_sons(Mesh* &mesh_pt) {}
211 
212  /// \short Number of vertex nodes in the element
213  unsigned nvertex_node() const
215 
216  /// \short Pointer to the j-th vertex node in the element
217  Node* vertex_node_pt(const unsigned& j) const
219 
220  /// \short Order of recovery shape functions for Z2 error estimation:
221  /// Same order as shape functions.
222  unsigned nrecovery_order() {return NNODE_1D-1;}
223 
224  /// \short No additional hanging node procedures are required
226 
227  };
228 
229 
230 //======================================================================
231 /// p-refineable version of 2D QLinearElasticityElement elements
232 //======================================================================
233 template<unsigned DIM>
235  public virtual RefineableLinearElasticityEquations<DIM>,
236  public virtual PRefineableQElement<DIM>
237 {
238  public:
239 
240  /// \short Constructor, simply call the other constructors
243  PRefineableQElement<DIM>(),
244  QLinearElasticityElement<DIM,2>()
245  {
246  // Set integration scheme
247  // (To avoid memory leaks in pre-build and p-refine where new
248  // integration schemes are created)
250  }
251 
252  /// Destructor (to avoid memory leaks)
254  {
255  delete this->integral_pt();
256  }
257 
258 
259  /// Broken copy constructor
261  {
262  BrokenCopy::broken_copy("PRefineableQLinearElasticityElement");
263  }
264 
265  /// Broken assignment operator
266 /* void operator=(const PRefineableQLinearElasticityElement<DIM>&)
267  {
268  BrokenCopy::broken_assign("PRefineableQLinearElasticityElement");
269  }*/
270 
271  void further_build();
272 
273  /// Number of continuously interpolated values: 1
274  unsigned ncont_interpolated_values() const {return 1;}
275 
276  /// \short Number of vertex nodes in the element
277  unsigned nvertex_node() const
279 
280  /// \short Pointer to the j-th vertex node in the element
281  Node* vertex_node_pt(const unsigned& j) const
283 
284  /// \short Order of recovery shape functions for Z2 error estimation:
285  /// - Same order as shape functions.
286  //unsigned nrecovery_order()
287  // {
288  // if(this->nnode_1d() < 4) {return (this->nnode_1d()-1);}
289  // else {return 3;}
290  // }
291  /// - Constant recovery order, since recovery order of the first element
292  /// is used for the whole mesh.
293  unsigned nrecovery_order() {return 3;}
294 
295  void compute_energy_error(std::ostream &outfile,
297  double& error, double& norm);
298 };
299 
300 
301 
302 //==============================================================
303 /// FaceGeometry of the 2D RefineableQLinearElasticityElement elements
304 //==============================================================
305  template<unsigned NNODE_1D>
307  public virtual QElement<1,NNODE_1D>
308  {
309  public:
310  //Make sure that we call the constructor of the QElement
311  //Only the Intel compiler seems to need this!
312  FaceGeometry() : QElement<1,NNODE_1D>() {}
313  };
314 
315 //==============================================================
316 /// FaceGeometry of the FaceGeometry of the 2D
317 ///RefineableQLinearElasticityElement
318 //==============================================================
319  template<unsigned NNODE_1D>
321  RefineableQLinearElasticityElement<2,NNODE_1D> > >:
322  public virtual PointElement
323  {
324  public:
325  //Make sure that we call the constructor of the QElement
326  //Only the Intel compiler seems to need this!
328  };
329 
330 
331 //==============================================================
332 /// FaceGeometry of the 3D RefineableQLinearElasticityElement elements
333 //==============================================================
334  template<unsigned NNODE_1D>
336  public virtual QElement<2,NNODE_1D>
337  {
338  public:
339  //Make sure that we call the constructor of the QElement
340  //Only the Intel compiler seems to need this!
341  FaceGeometry() : QElement<2,NNODE_1D>() {}
342  };
343 
344 //==============================================================
345 /// FaceGeometry of the FaceGeometry of the 3D
346 /// RefineableQLinearElasticityElement
347 //==============================================================
348  template<unsigned NNODE_1D>
350  RefineableQLinearElasticityElement<3,NNODE_1D> > >:
351  public virtual QElement<1,NNODE_1D>
352  {
353  public:
354  //Make sure that we call the constructor of the QElement
355  //Only the Intel compiler seems to need this!
356  FaceGeometry() : QElement<1,NNODE_1D>() {}
357  };
358 
359 
360 }
361 
362 #endif
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
~PRefineableQLinearElasticityElement()
Destructor (to avoid memory leaks)
Class for Refineable LinearElasticity equations.
void further_setup_hanging_nodes()
No additional hanging node procedures are required.
virtual unsigned u_index_linear_elasticity(const unsigned i) const
Return the index at which the i-th unknown displacement component is stored. The default value...
unsigned nvertex_node() const
Number of vertex nodes in the element.
cstr elem_len * i
Definition: cfortran.h:607
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain tensor.
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3148
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
char t
Definition: cfortran.h:572
unsigned nvertex_node() const
Number of vertex nodes in the element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation:
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: DIM.
ElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the current interpolated values (displacements). Note: Given the generality of the interface (thi...
p-refineable version of 2D QLinearElasticityElement elements
void rebuild_from_sons(Mesh *&mesh_pt)
Empty rebuild from sons, no need to reconstruct anything here.
bool Unsteady
Flag that switches inertia on/off.
void interpolated_u_linear_elasticity(const Vector< double > &s, Vector< double > &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
void compute_energy_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_grad_pt, double &error, double &norm)
Get error against and norm of exact solution.
static char t char * s
Definition: cfortran.h:572
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements...
Definition: elements.h:2371
PRefineableQLinearElasticityElement()
Constructor, simply call the other constructors.
bool is_inertia_enabled() const
Access function to flag that switches inertia on/off (const version)
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void fill_in_generic_contribution_to_residuals_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overloaded helper function to take hanging nodes into account.
virtual unsigned nvertex_node() const
Definition: elements.h:2362
double *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
Class for refineable QLinearElasticityElement elements.
PRefineableQLinearElasticityElement(const PRefineableQLinearElasticityElement< DIM > &dummy)
Broken copy constructor.
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
BodyForceFctPt Body_force_fct_pt
Pointer to body force function.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
A general mesh class.
Definition: mesh.h:74
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
void further_build()
Further build function, pass the pointers down to the sons.
ElasticityTensor *& elasticity_tensor_pt()
Return the pointer to the elasticity_tensor.