refineable_time_harmonic_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_TIME_HARMONIC_LINEAR_ELASTICITY_ELEMENTS_HEADER
30#define OOMPH_REFINEABLE_TIME_HARMONIC_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 TimeHarmonicLinearElasticity equations
42 //========================================================================
43 template<unsigned DIM>
45 : public virtual TimeHarmonicLinearElasticityEquations<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(2 * 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 std::complex<unsigned> u_nodal_index =
83 for (unsigned l = 0; l < n_node; l++)
84 {
85 values[i] += this->nodal_value(t, l, u_nodal_index.real()) * psi(l);
86 values[i + DIM] +=
87 this->nodal_value(t, l, u_nodal_index.imag()) * psi(l);
88 }
89 }
90 }
91
92 /// 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 this->get_interpolated_values(0, s, values);
100 }
101
102 /// Number of 'flux' terms for Z2 error estimation
104 {
105 // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
106 return 2 * (DIM + DIM * (DIM - 1) / 2);
107 }
108
109 /// Get 'flux' for Z2 error recovery: Upper triangular entries
110 /// in strain tensor.
112 {
113#ifdef PARANOID
114 unsigned num_entries = 2 * (DIM + ((DIM * DIM) - DIM) / 2);
115 if (flux.size() != num_entries)
116 {
117 std::ostringstream error_message;
118 error_message << "The flux vector has the wrong number of entries, "
119 << flux.size() << ", whereas it should be " << num_entries
120 << std::endl;
121 throw OomphLibError(error_message.str(),
122 OOMPH_CURRENT_FUNCTION,
123 OOMPH_EXCEPTION_LOCATION);
124 }
125#endif
126
127 // Get strain matrix
129 this->get_strain(s, strain);
130
131 // Pack into flux Vector
132 unsigned icount = 0;
133
134 // Start with diagonal terms
135 for (unsigned i = 0; i < DIM; i++)
136 {
137 flux[icount] = strain(i, i).real();
138 icount++;
139 flux[icount] = strain(i, i).imag();
140 icount++;
141 }
142
143 // Off diagonals row by row
144 for (unsigned i = 0; i < DIM; i++)
145 {
146 for (unsigned j = i + 1; j < DIM; j++)
147 {
148 flux[icount] = strain(i, j).real();
149 icount++;
150 flux[icount] = strain(i, j).imag();
151 icount++;
152 }
153 }
154 }
155
156 /// Number of continuously interpolated values: 2*DIM
158 {
159 return 2 * DIM;
160 }
161
162 /// Further build function, pass the pointers down to the sons
164 {
166 cast_father_element_pt =
168 this->father_element_pt());
169
170 // Set pointer to body force function
171 this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
172
173 // Set pointer to the contitutive law
175 cast_father_element_pt->elasticity_tensor_pt();
176
177 // Set the frequency
178 this->Omega_sq_pt = cast_father_element_pt->omega_sq_pt();
179 }
180
181
182 private:
183 /// Overloaded helper function to take hanging nodes into account
185 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
186 };
187
188
189 /// ///////////////////////////////////////////////////////////////////
190 /// ///////////////////////////////////////////////////////////////////
191 /// ///////////////////////////////////////////////////////////////////
192
193 //========================================================================
194 /// Class for refineable QTimeHarmonicLinearElasticityElement elements
195 //========================================================================
196 template<unsigned DIM, unsigned NNODE_1D>
198 : public virtual QTimeHarmonicLinearElasticityElement<DIM, NNODE_1D>,
200 public virtual RefineableQElement<DIM>
201 {
202 public:
203 /// Constructor:
205 : QTimeHarmonicLinearElasticityElement<DIM, NNODE_1D>(),
208 RefineableQElement<DIM>()
209 {
210 }
211
212 /// Empty rebuild from sons, no need to reconstruct anything here
213 void rebuild_from_sons(Mesh*& mesh_pt) {}
214
215 /// Number of vertex nodes in the element
216 unsigned nvertex_node() const
217 {
219 NNODE_1D>::nvertex_node();
220 }
221
222 /// Pointer to the j-th vertex node in the element
223 Node* vertex_node_pt(const unsigned& j) const
224 {
226 NNODE_1D>::vertex_node_pt(j);
227 }
228
229 /// Order of recovery shape functions for Z2 error estimation:
230 /// Same order as shape functions.
232 {
233 return NNODE_1D - 1;
234 }
235
236 /// No additional hanging node procedures are required
238 };
239
240
241 //==============================================================
242 /// FaceGeometry of the 2D
243 /// RefineableQTimeHarmonicLinearElasticityElement elements
244 //==============================================================
245 template<unsigned NNODE_1D>
248 : public virtual QElement<1, NNODE_1D>
249 {
250 public:
251 // Make sure that we call the constructor of the QElement
252 // Only the Intel compiler seems to need this!
253 FaceGeometry() : QElement<1, NNODE_1D>() {}
254 };
255
256 //==============================================================
257 /// FaceGeometry of the FaceGeometry of the 2D
258 /// RefineableQTimeHarmonicLinearElasticityElement
259 //==============================================================
260 template<unsigned NNODE_1D>
263 : public virtual PointElement
264 {
265 public:
266 // Make sure that we call the constructor of the QElement
267 // Only the Intel compiler seems to need this!
269 };
270
271
272 //==============================================================
273 /// FaceGeometry of the 3D RefineableQTimeHarmonicLinearElasticityElement
274 /// elements
275 //==============================================================
276 template<unsigned NNODE_1D>
279 : public virtual QElement<2, NNODE_1D>
280 {
281 public:
282 // Make sure that we call the constructor of the QElement
283 // Only the Intel compiler seems to need this!
284 FaceGeometry() : QElement<2, NNODE_1D>() {}
285 };
286
287 //==============================================================
288 /// FaceGeometry of the FaceGeometry of the 3D
289 /// RefineableQTimeHarmonicLinearElasticityElement
290 //==============================================================
291 template<unsigned NNODE_1D>
294 : public virtual QElement<1, NNODE_1D>
295 {
296 public:
297 // Make sure that we call the constructor of the QElement
298 // Only the Intel compiler seems to need this!
299 FaceGeometry() : QElement<1, NNODE_1D>() {}
300 };
301
302
303} // namespace oomph
304
305#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
Definition: matrices.h:386
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
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 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
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....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
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.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void further_setup_hanging_nodes()
No additional hanging node procedures are required.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void rebuild_from_sons(Mesh *&mesh_pt)
Empty rebuild from sons, no need to reconstruct anything here.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 2*DIM.
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 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_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_time_harmonic_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overloaded helper function to take hanging nodes into account.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
TimeHarmonicElasticityTensor *& elasticity_tensor_pt()
Return the pointer to the elasticity_tensor.
double *& omega_sq_pt()
Access function for square of non-dim frequency.
TimeHarmonicElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
virtual std::complex< unsigned > u_index_time_harmonic_linear_elasticity(const unsigned i) const
Return the index at which the i-th real or imag unknown displacement component is stored....
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double > > &strain) const
Return the strain tensor.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...