pml_time_harmonic_elasticity_tensor.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 objects representing the fourth-rank elasticity tensor
27// for linear elasticity problems
28
29// Include guards to prevent multiple inclusion of the header
30#ifndef OOMPH_PML_TIME_HARMONIC_ELASTICITY_TENSOR_HEADER
31#define OOMPH_PML_TIME_HARMONIC_ELASTICITY_TENSOR_HEADER
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <oomph-lib-config.h>
36#endif
37
38#include "../generic/oomph_utilities.h"
39#include "math.h"
40#include <complex>
41
42namespace oomph
43{
44 //=====================================================================
45 /// A base class that represents the fourth-rank elasticity tensor
46 /// \f$E_{ijkl}\f$ defined such that
47 /// \f[\tau_{ij} = E_{ijkl} e_{kl},\f]
48 /// where \f$e_{ij}\f$ is the infinitessimal (Cauchy) strain tensor
49 /// and \f$\tau_{ij}\f$ is the stress tensor. The symmetries of the
50 /// tensor are such that
51 /// \f[E_{ijkl} = E_{jikl} = E_{ijlk} = E_{klij}\f]
52 /// and thus there are relatively few independent components. These
53 /// symmetries are included in the definition of the object so that
54 /// non-physical symmetries cannot be accidentally imposed.
55 //=====================================================================
57 {
58 /// Translation table from the four indices to the corresponding
59 /// independent component
60 static const unsigned Index[3][3][3][3];
61
62 protected:
63 /// Member function that returns the i-th independent component of the
64 /// elasticity tensor
65 virtual inline std::complex<double> independent_component(
66 const unsigned& i) const
67 {
68 return std::complex<double>(0.0, 0.0);
69 }
70
71
72 /// Helper range checking function
73 /// (Note that this only captures over-runs in 3D but
74 /// errors are likely to be caught in evaluation of the
75 /// stress and strain tensors anyway...)
76 void range_check(const unsigned& i,
77 const unsigned& j,
78 const unsigned& k,
79 const unsigned& l) const
80 {
81 if ((i > 2) || (j > 2) || (k > 2) || (l > 2))
82 {
83 std::ostringstream error_message;
84 if (i > 2)
85 {
86 error_message << "Range Error : Index 1 " << i
87 << " is not in the range (0,2)";
88 }
89 if (j > 2)
90 {
91 error_message << "Range Error : Index 2 " << j
92 << " is not in the range (0,2)";
93 }
94
95 if (k > 2)
96 {
97 error_message << "Range Error : Index 2 " << k
98 << " is not in the range (0,2)";
99 }
100
101 if (l > 2)
102 {
103 error_message << "Range Error : Index 4 " << l
104 << " is not in the range (0,2)";
105 }
106
107 // Throw the error
108 throw OomphLibError(error_message.str(),
109 OOMPH_CURRENT_FUNCTION,
110 OOMPH_EXCEPTION_LOCATION);
111 }
112 }
113
114
115 /// Empty Constructor
117
118 public:
119 /// Empty virtual Destructor
121
122 public:
123 /// Return the appropriate independent component
124 /// via the index translation scheme (const version).
125 std::complex<double> operator()(const unsigned& i,
126 const unsigned& j,
127 const unsigned& k,
128 const unsigned& l) const
129 {
130 // Range check
131#ifdef PARANOID
132 range_check(i, j, k, l);
133#endif
134 return independent_component(Index[i][j][k][l]);
135 }
136 };
137
138
139 //===================================================================
140 /// An isotropic elasticity tensor defined in terms of Young's modulus
141 /// and Poisson's ratio. The elasticity tensor is assumed to be
142 /// non-dimensionalised on some reference value for Young's modulus
143 /// so the value provided to the constructor (if any) is to be
144 /// interpreted as the ratio of the actual Young's modulus to the
145 /// Young's modulus used to non-dimensionalise the stresses/tractions
146 /// in the governing equations.
147 //===================================================================
150 {
151 /// Storage for the independent components of the elasticity tensor
152 std::complex<double> C[4];
153
154 /// Translation scheme between the 21 independent components of the
155 /// general elasticity tensor and the isotropic case
156 static const unsigned StaticIndex[21];
157
158 /// Poisson's ratio
159 double Nu;
160
161 public:
162 /// Constructor. Passing in the values of the Poisson's ratio
163 /// and Young's modulus (interpreted as the ratio of the actual
164 /// Young's modulus to the Young's modulus (or other reference stiffness)
165 /// used to non-dimensionalise stresses and tractions in the governing
166 /// equations).
167 PMLTimeHarmonicIsotropicElasticityTensor(const double& nu, const double& E)
169 {
170 // Set the three indepdent components
171 C[0] = std::complex<double>(0.0, 0.0);
172 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
173 double mu = E / (2.0 * (1.0 + nu));
174 this->set_lame_coefficients(lambda, mu);
175 }
176
177 /// Constructor. Passing in the value of the Poisson's ratio.
178 /// Stresses and tractions in the governing equations are assumed
179 /// to have been non-dimensionalised on Young's modulus.
182 {
183 // Set the three indepdent components
184 C[0] = std::complex<double>(0.0, 0.0);
185
186 // reference value
187 double E = 1.0;
188 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
189 double mu = E / (2.0 * (1.0 + nu));
190 this->set_lame_coefficients(lambda, mu);
191 }
192
193 /// Poisson's ratio
194 double nu() const
195 {
196 return Nu;
197 }
198
199 /// Update parameters: Specify values of the Poisson's ratio
200 /// and (optionally) Young's modulus (interpreted as the ratio of the actual
201 /// Young's modulus to the Young's modulus (or other reference stiffness)
202 /// used to non-dimensionalise stresses and tractions in the governing
203 /// equations).
204 void update_constitutive_parameters(const double& nu, const double& E = 1.0)
205 {
206 // Set the three indepdent components
207 C[0] = std::complex<double>(0.0, 0.0);
208 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
209 double mu = E / (2.0 * (1.0 + nu));
210 this->set_lame_coefficients(lambda, mu);
211 }
212
213
214 /// Overload the independent coefficient function
215 inline std::complex<double> independent_component(const unsigned& i) const
216 {
217 return C[StaticIndex[i]];
218 }
219
220
221 private:
222 // Set the values of the lame coefficients
223 void set_lame_coefficients(const double& lambda, const double& mu)
224 {
225 C[1] = lambda + 2.0 * mu;
226 C[2] = lambda;
227 C[3] = mu;
228 }
229 };
230
231} // namespace oomph
232#endif
cstr elem_len * i
Definition: cfortran.h:603
An OomphLibError object which should be thrown when an run-time error is encountered....
A base class that represents the fourth-rank elasticity tensor defined such that.
static const unsigned Index[3][3][3][3]
Translation table from the four indices to the corresponding independent component.
virtual ~PMLTimeHarmonicElasticityTensor()
Empty virtual Destructor.
virtual std::complex< double > independent_component(const unsigned &i) const
Member function that returns the i-th independent component of the elasticity tensor.
std::complex< double > operator()(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l) const
Return the appropriate independent component via the index translation scheme (const version).
void range_check(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l) const
Helper range checking function (Note that this only captures over-runs in 3D but errors are likely to...
An isotropic elasticity tensor defined in terms of Young's modulus and Poisson's ratio....
PMLTimeHarmonicIsotropicElasticityTensor(const double &nu, const double &E)
Constructor. Passing in the values of the Poisson's ratio and Young's modulus (interpreted as the rat...
PMLTimeHarmonicIsotropicElasticityTensor(const double &nu)
Constructor. Passing in the value of the Poisson's ratio. Stresses and tractions in the governing equ...
void update_constitutive_parameters(const double &nu, const double &E=1.0)
Update parameters: Specify values of the Poisson's ratio and (optionally) Young's modulus (interprete...
std::complex< double > C[4]
Storage for the independent components of the elasticity tensor.
static const unsigned StaticIndex[21]
Translation scheme between the 21 independent components of the general elasticity tensor and the iso...
void set_lame_coefficients(const double &lambda, const double &mu)
std::complex< double > independent_component(const unsigned &i) const
Overload the independent coefficient function.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...