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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for refineable linear elasticity elements
27
28// Include guard to prevent multiple inclusions of this header
29#ifndef OOMPH_REFINEABLE_LINEAR_ELASTICITY_ELEMENTS_HEADER
30#define OOMPH_REFINEABLE_LINEAR_ELASTICITY_ELEMENTS_HEADER
31
32// oomph-lib headers
34#include "../generic/refineable_quad_element.h"
35#include "../generic/refineable_brick_element.h"
36#include "../generic/error_estimator.h"
37
38namespace oomph
39{
40 //========================================================================
41 /// Class for Refineable LinearElasticity equations
42 //========================================================================
43 template<unsigned DIM>
45 : public virtual LinearElasticityEquations<DIM>,
46 public virtual RefineableElement,
47 public virtual ElementWithZ2ErrorEstimator
48 {
49 public:
50 /// Constructor
55 {
56 }
57
58
59 /// Get the function value u in Vector.
60 /// Note: Given the generality of the interface (this function
61 /// is usually called from black-box documentation or interpolation
62 /// routines), the values Vector sets its own size in here.
63 void get_interpolated_values(const unsigned& t,
64 const Vector<double>& s,
65 Vector<double>& values)
66 {
67 // Create enough initialised storage
68 values.resize(DIM, 0.0);
69
70 // Find out how many nodes there are
71 unsigned n_node = this->nnode();
72
73 // Shape functions
74 Shape psi(n_node);
75 this->shape(s, psi);
76
77 // Calculate displacements
78 for (unsigned i = 0; i < DIM; i++)
79 {
80 // Get the index at which the i-th velocity is stored
81 unsigned u_nodal_index = this->u_index_linear_elasticity(i);
82 for (unsigned l = 0; l < n_node; l++)
83 {
84 values[i] += this->nodal_value(t, l, u_nodal_index) * psi(l);
85 }
86 }
87 }
88
89 /// Get the current interpolated values (displacements).
90 /// Note: Given the generality of the interface (this function
91 /// is usually called from black-box documentation or interpolation
92 /// routines),the values Vector sets its own size in here.
94 Vector<double>& values)
95 {
96 values.resize(DIM);
97 this->interpolated_u_linear_elasticity(s, values);
98 }
99
100 /// Number of 'flux' terms for Z2 error estimation
102 {
103 // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
104 return DIM + DIM * (DIM - 1) / 2;
105 }
106
107 /// Get 'flux' for Z2 error recovery: Upper triangular entries
108 /// in strain tensor.
110 {
111#ifdef PARANOID
112 unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
113 if (flux.size() != num_entries)
114 {
115 std::ostringstream error_message;
116 error_message << "The flux vector has the wrong number of entries, "
117 << flux.size() << ", whereas it should be " << num_entries
118 << std::endl;
119 throw OomphLibError(error_message.str(),
120 OOMPH_CURRENT_FUNCTION,
121 OOMPH_EXCEPTION_LOCATION);
122 }
123#endif
124
125 // Get strain matrix
126 DenseMatrix<double> strain(DIM);
127 this->get_strain(s, strain);
128
129 // Pack into flux Vector
130 unsigned icount = 0;
131
132 // Start with diagonal terms
133 for (unsigned i = 0; i < DIM; i++)
134 {
135 flux[icount] = strain(i, i);
136 icount++;
137 }
138
139 // Off diagonals row by row
140 for (unsigned i = 0; i < DIM; i++)
141 {
142 for (unsigned j = i + 1; j < DIM; j++)
143 {
144 flux[icount] = strain(i, j);
145 icount++;
146 }
147 }
148 }
149
150 /// Number of continuously interpolated values: DIM
152 {
153 return DIM;
154 }
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
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 /// Overloaded helper function to take hanging nodes into account
181 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
182 };
183
184
185 /// ///////////////////////////////////////////////////////////////////
186 /// ///////////////////////////////////////////////////////////////////
187 /// ///////////////////////////////////////////////////////////////////
188
189 //========================================================================
190 /// Class for refineable QLinearElasticityElement elements
191 //========================================================================
192 template<unsigned DIM, unsigned NNODE_1D>
194 : public virtual QLinearElasticityElement<DIM, NNODE_1D>,
195 public virtual RefineableLinearElasticityEquations<DIM>,
196 public virtual RefineableQElement<DIM>
197 {
198 public:
199 /// Constructor:
201 : QLinearElasticityElement<DIM, NNODE_1D>(),
204 RefineableQElement<DIM>()
205 {
206 }
207
208 /// Empty rebuild from sons, no need to reconstruct anything here
209 void rebuild_from_sons(Mesh*& mesh_pt) {}
210
211 /// Number of vertex nodes in the element
212 unsigned nvertex_node() const
213 {
215 }
216
217 /// Pointer to the j-th vertex node in the element
218 Node* vertex_node_pt(const unsigned& j) const
219 {
221 }
222
223 /// Order of recovery shape functions for Z2 error estimation:
224 /// Same order as shape functions.
226 {
227 return NNODE_1D - 1;
228 }
229
230 /// No additional hanging node procedures are required
232 };
233
234
235 //======================================================================
236 /// p-refineable version of 2D QLinearElasticityElement elements
237 //======================================================================
238 template<unsigned DIM>
240 : public QLinearElasticityElement<DIM, 2>,
241 public virtual RefineableLinearElasticityEquations<DIM>,
242 public virtual PRefineableQElement<DIM>
243 {
244 public:
245 /// Constructor, simply call the other constructors
249 PRefineableQElement<DIM>(),
251 {
252 // Set integration scheme
253 // (To avoid memory leaks in pre-build and p-refine where new
254 // integration schemes are created)
256 }
257
258 /// Destructor (to avoid memory leaks)
260 {
261 delete this->integral_pt();
262 }
263
264
265 /// Broken copy constructor
267 const PRefineableQLinearElasticityElement<DIM>& dummy) = delete;
268
269 /// Broken assignment operator
270 /* void operator=(const PRefineableQLinearElasticityElement<DIM>&) =
271 * delete;*/
272
273 void further_build();
274
275 /// Number of continuously interpolated values: 1
277 {
278 return 1;
279 }
280
281 /// Number of vertex nodes in the element
282 unsigned nvertex_node() const
283 {
285 }
286
287 /// Pointer to the j-th vertex node in the element
288 Node* vertex_node_pt(const unsigned& j) const
289 {
291 }
292
293 /// Order of recovery shape functions for Z2 error estimation:
294 /// - Same order as shape functions.
295 // unsigned nrecovery_order()
296 // {
297 // if(this->nnode_1d() < 4) {return (this->nnode_1d()-1);}
298 // else {return 3;}
299 // }
300 /// - Constant recovery order, since recovery order of the first element
301 /// is used for the whole mesh.
303 {
304 return 3;
305 }
306
308 std::ostream& outfile,
310 double& error,
311 double& norm);
312 };
313
314
315 //==============================================================
316 /// FaceGeometry of the 2D RefineableQLinearElasticityElement elements
317 //==============================================================
318 template<unsigned NNODE_1D>
320 : public virtual QElement<1, NNODE_1D>
321 {
322 public:
323 // Make sure that we call the constructor of the QElement
324 // Only the Intel compiler seems to need this!
325 FaceGeometry() : QElement<1, NNODE_1D>() {}
326 };
327
328 //==============================================================
329 /// FaceGeometry of the FaceGeometry of the 2D
330 /// RefineableQLinearElasticityElement
331 //==============================================================
332 template<unsigned NNODE_1D>
335 : public virtual PointElement
336 {
337 public:
338 // Make sure that we call the constructor of the QElement
339 // Only the Intel compiler seems to need this!
341 };
342
343
344 //==============================================================
345 /// FaceGeometry of the 3D RefineableQLinearElasticityElement elements
346 //==============================================================
347 template<unsigned NNODE_1D>
349 : public virtual QElement<2, NNODE_1D>
350 {
351 public:
352 // Make sure that we call the constructor of the QElement
353 // Only the Intel compiler seems to need this!
354 FaceGeometry() : QElement<2, NNODE_1D>() {}
355 };
356
357 //==============================================================
358 /// FaceGeometry of the FaceGeometry of the 3D
359 /// RefineableQLinearElasticityElement
360 //==============================================================
361 template<unsigned NNODE_1D>
364 : public virtual QElement<1, NNODE_1D>
365 {
366 public:
367 // Make sure that we call the constructor of the QElement
368 // Only the Intel compiler seems to need this!
369 FaceGeometry() : QElement<1, NNODE_1D>() {}
370 };
371
372
373} // namespace oomph
374
375#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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:2593
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2491
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...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3210
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:2500
Class for multidimensional Gauss Lobatto Legendre integration rules empty - just establishes template...
Definition: integral.h:1281
ElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
ElasticityTensor *& elasticity_tensor_pt()
Return the pointer to the elasticity_tensor.
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
double *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
BodyForceFctPt Body_force_fct_pt
Pointer to body force function.
bool is_inertia_enabled() const
Access function to flag that switches inertia on/off (const version)
void interpolated_u_linear_elasticity(const Vector< double > &s, Vector< double > &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
bool Unsteady
Flag that switches inertia on/off.
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,...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
A general mesh class.
Definition: mesh.h:67
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
A class that is used to template the p-refineable Q elements by dimension. It's really nothing more t...
Definition: Qelements.h:2274
p-refineable version of 2D QLinearElasticityElement elements
PRefineableQLinearElasticityElement()
Constructor, simply call the other constructors.
PRefineableQLinearElasticityElement(const PRefineableQLinearElasticityElement< DIM > &dummy)=delete
Broken copy constructor.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void compute_energy_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_grad_pt, double &error, double &norm)
Get error against and norm of exact solution.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
~PRefineableQLinearElasticityElement()
Destructor (to avoid memory leaks)
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation:
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Class for Refineable LinearElasticity equations.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: DIM.
void further_build()
Further build function, pass the pointers down to the sons.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in 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...
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...
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.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
void further_setup_hanging_nodes()
No additional hanging node procedures are required.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void rebuild_from_sons(Mesh *&mesh_pt)
Empty rebuild from sons, no need to reconstruct anything here.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
unsigned nvertex_node() const
Number of vertex nodes in the element.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...