poroelasticity/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_POROELASTICITY_TENSOR_HEADER
31#define OOMPH_POROELASTICITY_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
40namespace oomph
41{
42 //=====================================================================
43 /// A base class that represents the fourth-rank elasticity tensor
44 /// \f$E_{ijkl}\f$ defined such that
45 /// \f[\tau_{ij} = E_{ijkl} e_{kl},\f]
46 /// where \f$e_{ij}\f$ is the infinitessimal (Cauchy) strain tensor
47 /// and \f$\tau_{ij}\f$ is the stress tensor. The symmetries of the
48 /// tensor are such that
49 /// \f[E_{ijkl} = E_{jikl} = E_{ijlk} = E_{klij}\f]
50 /// and thus there are relatively few independent components. These
51 /// symmetries are included in the definition of the object so that
52 /// non-physical symmetries cannot be accidentally imposed.
53 //=====================================================================
54 class ElasticityTensor
55 {
56 /// Translation table from the four indices to the corresponding
57 /// independent component
58 static const unsigned Index[3][3][3][3];
59
60 protected:
61 /// Member function that returns the i-th independent component of the
62 /// elasticity tensor
63 virtual inline double independent_component(const unsigned& i) const
64 {
65 return 0.0;
66 }
67
68
69 /// Helper range checking function
70 /// (Note that this only captures over-runs in 3D but
71 /// errors are likely to be caught in evaluation of the
72 /// stress and strain tensors anyway...)
73 void range_check(const unsigned& i,
74 const unsigned& j,
75 const unsigned& k,
76 const unsigned& l) const
77 {
78 if ((i > 2) || (j > 2) || (k > 2) || (l > 2))
79 {
80 std::ostringstream error_message;
81 if (i > 2)
82 {
83 error_message << "Range Error : Index 1 " << i
84 << " is not in the range (0,2)";
85 }
86 if (j > 2)
87 {
88 error_message << "Range Error : Index 2 " << j
89 << " is not in the range (0,2)";
90 }
91
92 if (k > 2)
93 {
94 error_message << "Range Error : Index 2 " << k
95 << " is not in the range (0,2)";
96 }
97
98 if (l > 2)
99 {
100 error_message << "Range Error : Index 4 " << l
101 << " is not in the range (0,2)";
102 }
103
104 // Throw the error
105 throw OomphLibError(error_message.str(),
106 OOMPH_CURRENT_FUNCTION,
107 OOMPH_EXCEPTION_LOCATION);
108 }
109 }
110
111
112 /// Empty Constructor
114
115 public:
116 /// Empty virtual Destructor
117 virtual ~ElasticityTensor() {}
118
119 public:
120 /// Return the appropriate independent component
121 /// via the index translation scheme (const version).
122 double operator()(const unsigned& i,
123 const unsigned& j,
124 const unsigned& k,
125 const unsigned& l) const
126 {
127 // Range check
128#ifdef PARANOID
129 range_check(i, j, k, l);
130#endif
131 return independent_component(Index[i][j][k][l]);
132 }
133 };
134
135
136 //===================================================================
137 /// An isotropic elasticity tensor defined in terms of Young's modulus
138 /// and Poisson's ratio. The elasticity tensor is assumed to be
139 /// non-dimensionalised on some reference value for Young's modulus
140 /// so the value provided to the constructor (if any) is to be
141 /// interpreted as the ratio of the actual Young's modulus to the
142 /// Young's modulus used to non-dimensionalise the stresses/tractions
143 /// in the governing equations.
144 //===================================================================
145 class IsotropicElasticityTensor : public ElasticityTensor
146 {
147 // Storage for the independent components of the elasticity tensor
148 double C[4];
149
150 // Translation scheme between the 21 independent components of the general
151 // elasticity tensor and the isotropic case
152 static const unsigned StaticIndex[21];
153
154 public:
155 /// Constructor. Passing in the values of the Poisson's ratio
156 /// and Young's modulus (interpreted as the ratio of the actual
157 /// Young's modulus to the Young's modulus (or other reference stiffness)
158 /// used to non-dimensionalise stresses and tractions in the governing
159 /// equations).
160 IsotropicElasticityTensor(const double& nu, const double& E)
162 {
163 // Set the three independent components
164 C[0] = 0.0;
165 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
166 double mu = E / (2.0 * (1.0 + nu));
167 this->set_lame_coefficients(lambda, mu);
168 }
169
170 /// Constructor. Passing in the value of the Poisson's ratio.
171 /// Stresses and tractions in the governing equations are assumed
172 /// to have been non-dimensionalised on Young's modulus.
174 {
175 // Set the three independent components
176 C[0] = 0.0;
177
178 // reference value
179 double E = 1.0;
180 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
181 double mu = E / (2.0 * (1.0 + nu));
182 this->set_lame_coefficients(lambda, mu);
183 }
184
185 /// Constructur. Passing in the values of the two lame
186 /// coefficients directly (interpreted as the ratios of these
187 /// quantities to a reference stiffness used to non-dimensionalised
189 {
190 // Set the three independent components
191 C[0] = 0.0;
192 this->set_lame_coefficients(lame[0], lame[1]);
193 }
194
195
196 /// Overload the independent coefficient function
197 inline double independent_component(const unsigned& i) const
198 {
199 return C[StaticIndex[i]];
200 }
201
202
203 private:
204 // Set the values of the lame coefficients
205 void set_lame_coefficients(const double& lambda, const double& mu)
206 {
207 C[1] = lambda + 2.0 * mu;
208 C[2] = lambda;
209 C[3] = mu;
210 }
211 };
212
213
214 //===================================================================
215 /// An isotropic elasticity tensor defined in terms of Young's modulus
216 /// and Poisson's ratio. The elasticity tensor is assumed to be
217 /// non-dimensionalised on some reference value for Young's modulus
218 /// so the value provided to the constructor (if any) is to be
219 /// interpreted as the ratio of the actual Young's modulus to the
220 /// Young's modulus used to non-dimensionalise the stresses/tractions
221 /// in the governing equations.
222 //===================================================================
224 {
225 // Storage for the independent components of the elasticity tensor
226 double C[3];
227
228 // Storage for the first Lame parameter
229 double Lambda;
230
231 // Storage for the second Lame parameter
232 double Mu;
233
234 // Translation scheme between the 21 independent components of the general
235 // elasticity tensor and the isotropic case
236 static const unsigned StaticIndex[21];
237
238 public:
239 /// Constructor. For use with incompressibility. Requires no
240 /// parameters since Poisson's ratio is fixed at 0.5 and lambda is set to a
241 /// dummy value of 0 (since it would be infinite)
243 {
244 C[0] = 0.0;
245 double E = 1.0;
246 double nu = 0.5;
247 double lambda = 0.0;
248 double mu = E / (2.0 * (1.0 + nu));
249 this->set_lame_coefficients(lambda, mu);
250 }
251
252 /// Constructor. Passing in the values of the Poisson's ratio
253 /// and Young's modulus (interpreted as the ratio of the actual
254 /// Young's modulus to the Young's modulus (or other reference stiffness)
255 /// used to non-dimensionalise stresses and tractions in the governing
256 /// equations).
257 DeviatoricIsotropicElasticityTensor(const double& nu, const double& E)
259 {
260 // Set the three independent components
261 C[0] = 0.0;
262 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
263 double mu = E / (2.0 * (1.0 + nu));
264 this->set_lame_coefficients(lambda, mu);
265 }
266
267 /// Constructor. Passing in the value of the Poisson's ratio.
268 /// Stresses and tractions in the governing equations are assumed
269 /// to have been non-dimensionalised on Young's modulus.
271 {
272 // Set the three independent components
273 C[0] = 0.0;
274
275 // reference value
276 double E = 1.0;
277 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
278 double mu = E / (2.0 * (1.0 + nu));
279 this->set_lame_coefficients(lambda, mu);
280 }
281
282 /// Constructur. Passing in the values of the two lame
283 /// coefficients directly (interpreted as the ratios of these
284 /// quantities to a reference stiffness used to non-dimensionalised
286 {
287 // Set the three independent components
288 C[0] = 0.0;
289 this->set_lame_coefficients(lame[0], lame[1]);
290 }
291
292
293 /// Overload the independent coefficient function
294 inline double independent_component(const unsigned& i) const
295 {
296 return C[StaticIndex[i]];
297 }
298
299 /// Accessor function for the first lame parameter
300 const double& lambda() const
301 {
302 return Lambda;
303 }
304
305 /// Accessor function for the second lame parameter
306 const double& mu() const
307 {
308 return Mu;
309 }
310
311 private:
312 // Set the values of the lame coefficients
313 void set_lame_coefficients(const double& lambda, const double& mu)
314 {
315 C[1] = 2.0 * mu;
316 C[2] = mu;
317
318 Lambda = lambda;
319 Mu = mu;
320 }
321 };
322
323} // namespace oomph
324#endif
cstr elem_len * i
Definition: cfortran.h:603
An isotropic elasticity tensor defined in terms of Young's modulus and Poisson's ratio....
static const unsigned StaticIndex[21]
Translation scheme for the deviatoric isotropic elasticity tensor.
DeviatoricIsotropicElasticityTensor(const double &nu)
Constructor. Passing in the value of the Poisson's ratio. Stresses and tractions in the governing equ...
void set_lame_coefficients(const double &lambda, const double &mu)
DeviatoricIsotropicElasticityTensor()
Constructor. For use with incompressibility. Requires no parameters since Poisson's ratio is fixed at...
const double & mu() const
Accessor function for the second lame parameter.
DeviatoricIsotropicElasticityTensor(const double &nu, const double &E)
Constructor. Passing in the values of the Poisson's ratio and Young's modulus (interpreted as the rat...
const double & lambda() const
Accessor function for the first lame parameter.
double independent_component(const unsigned &i) const
Overload the independent coefficient function.
DeviatoricIsotropicElasticityTensor(const Vector< double > &lame)
Constructur. Passing in the values of the two lame coefficients directly (interpreted as the ratios o...
A base class that represents the fourth-rank elasticity tensor defined such that.
virtual ~ElasticityTensor()
Empty virtual Destructor.
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).
static const unsigned Index[3][3][3][3]
Translation table from the four indices to the corresponding independent component.
virtual double independent_component(const unsigned &i) const
Member function that returns the i-th independent component of the elasticity tensor.
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...
IsotropicElasticityTensor(const double &nu, const double &E)
Constructor. Passing in the values of the Poisson's ratio and Young's modulus (interpreted as the rat...
IsotropicElasticityTensor(const Vector< double > &lame)
Constructur. Passing in the values of the two lame coefficients directly (interpreted as the ratios o...
void set_lame_coefficients(const double &lambda, const double &mu)
IsotropicElasticityTensor(const double &nu)
Constructor. Passing in the value of the Poisson's ratio. Stresses and tractions in the governing equ...
double independent_component(const unsigned &i) const
Overload the independent coefficient function.
static const unsigned StaticIndex[21]
Translation scheme for the isotropic elasticity tensor.
An OomphLibError object which should be thrown when an run-time error is encountered....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...