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-2023 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 
38 namespace 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
103  unsigned num_Z2_flux_terms()
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
157  unsigned ncont_interpolated_values() const
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
174  this->Elasticity_tensor_pt =
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.
231  unsigned nrecovery_order()
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
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
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.
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 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
TimeHarmonicElasticityTensor *& elasticity_tensor_pt()
Return the pointer to the elasticity_tensor.
TimeHarmonicElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
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....
double *& omega_sq_pt()
Access function for square of non-dim frequency.
void get_strain(const Vector< double > &s, DenseMatrix< std::complex< double >> &strain) const
Return the strain tensor.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...