linear_elasticity/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_ELASTICITY_TENSOR_HEADER
31#define OOMPH_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
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 //=====================================================================
55 {
56 protected:
57 /// Translation table from the four indices to the corresponding
58 /// independent component
59 static const unsigned Index[3][3][3][3];
60
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 /// Helper range checking function
69 /// (Note that this only captures over-runs in 3D but
70 /// errors are likely to be caught in evaluation of the
71 /// stress and strain tensors anyway...)
72 void range_check(const unsigned& i,
73 const unsigned& j,
74 const unsigned& k,
75 const unsigned& l) const
76 {
77 if ((i > 2) || (j > 2) || (k > 2) || (l > 2))
78 {
79 std::ostringstream error_message;
80 if (i > 2)
81 {
82 error_message << "Range Error : Index 1 " << i
83 << " is not in the range (0,2)";
84 }
85 if (j > 2)
86 {
87 error_message << "Range Error : Index 2 " << j
88 << " is not in the range (0,2)";
89 }
90
91 if (k > 2)
92 {
93 error_message << "Range Error : Index 2 " << k
94 << " is not in the range (0,2)";
95 }
96
97 if (l > 2)
98 {
99 error_message << "Range Error : Index 4 " << l
100 << " is not in the range (0,2)";
101 }
102
103 // Throw the error
104 throw OomphLibError(error_message.str(),
105 OOMPH_CURRENT_FUNCTION,
106 OOMPH_EXCEPTION_LOCATION);
107 }
108 }
109
110 /// Empty Constructor
112
113 public:
114 /// Empty virtual Destructor
115 virtual ~ElasticityTensor() {}
116
117 public:
118 /// Return the appropriate independent component
119 /// via the index translation scheme (const version).
120 double operator()(const unsigned& i,
121 const unsigned& j,
122 const unsigned& k,
123 const unsigned& l) const
124 {
125 // Range check
126#ifdef PARANOID
127 range_check(i, j, k, l);
128#endif
129 return independent_component(Index[i][j][k][l]);
130 }
131
132 /// Allow the values to be set (virtual function that must be
133 /// overloaded if values can be set directly
134 virtual void set_value(const unsigned& i,
135 const unsigned& j,
136 const unsigned& k,
137 const unsigned& l,
138 const double& value)
139 {
140 std::stringstream error_stream;
141 error_stream << "Broken base implementation.\n"
142 << "Must be overloaded for specific storage schemes.\n";
143 throw OomphLibError(error_stream.str().c_str(),
144 OOMPH_CURRENT_FUNCTION,
145 OOMPH_EXCEPTION_LOCATION);
146 }
147 };
148
149
150 //===================================================================
151 /// An isotropic elasticity tensor defined in terms of Young's modulus
152 /// and Poisson's ratio. The elasticity tensor is assumed to be
153 /// non-dimensionalised on some reference value for Young's modulus
154 /// so the value provided to the constructor (if any) is to be
155 /// interpreted as the ratio of the actual Young's modulus to the
156 /// Young's modulus used to non-dimensionalise the stresses/tractions
157 /// in the governing equations.
158 //===================================================================
160 {
161 // Storage for the independent components of the elasticity tensor
162 double C[4];
163
164 // Translation scheme between the 21 independent components of the general
165 // elasticity tensor and the isotropic case
166 static const unsigned StaticIndex[21];
167
168 public:
169 /// Constructor. Passing in the values of the Poisson's ratio
170 /// and Young's modulus (interpreted as the ratio of the actual
171 /// Young's modulus to the Young's modulus (or other reference stiffness)
172 /// used to non-dimensionalise stresses and tractions in the governing
173 /// equations).
174 IsotropicElasticityTensor(const double& nu, const double& E)
176 {
177 // Set the three indepdent components
178 C[0] = 0.0;
179 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
180 double mu = E / (2.0 * (1.0 + nu));
181 this->set_lame_coefficients(lambda, mu);
182 }
183
184 /// Constructor. Passing in the value of the Poisson's ratio.
185 /// Stresses and tractions in the governing equations are assumed
186 /// to have been non-dimensionalised on Young's modulus.
188 {
189 // Set the three indepdent components
190 C[0] = 0.0;
191
192 // reference value
193 double E = 1.0;
194 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
195 double mu = E / (2.0 * (1.0 + nu));
196 this->set_lame_coefficients(lambda, mu);
197 }
198
199 /// Constructur. Passing in the values of the two lame
200 /// coefficients directly (interpreted as the ratios of these
201 /// quantities to a reference stiffness used to non-dimensionalised
203 {
204 // Set the three independent componens
205 C[0] = 0.0;
206 this->set_lame_coefficients(lame[0], lame[1]);
207 }
208
209 /// Update parameters: Specify values of the Poisson's ratio
210 /// and (optionally) Young's modulus (interpreted as the ratio of the actual
211 /// Young's modulus to the Young's modulus (or other reference stiffness)
212 /// used to non-dimensionalise stresses and tractions in the governing
213 /// equations).
214 void update_constitutive_parameters(const double& nu, const double& E = 1.0)
215 {
216 // Set the three indepdent components
217 C[0] = 0.0;
218 double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
219 double mu = E / (2.0 * (1.0 + nu));
220 this->set_lame_coefficients(lambda, mu);
221 }
222
223
224 /// Overload the independent coefficient function
225 inline double independent_component(const unsigned& i) const
226 {
227 return C[StaticIndex[i]];
228 }
229
230
231 private:
232 // Set the values of the lame coefficients
233 void set_lame_coefficients(const double& lambda, const double& mu)
234 {
235 C[1] = lambda + 2.0 * mu;
236 C[2] = lambda;
237 C[3] = mu;
238 }
239 };
240
241
242 //========================================================================
243 /// A general elasticity tensor that provides storage for all 21 independent
244 /// components
245 //===================================================================
247 {
248 // Storage for the independent components of the elasticity tensor
249 double C[21];
250
251 public:
252 /// Empy Constructor.
254 {
255 // Initialise all independent components to zero
256 for (unsigned i = 0; i < 21; i++)
257 {
258 C[i] = 0.0;
259 }
260 }
261
262 /// Overload the independent coefficient function
263 inline double independent_component(const unsigned& i) const
264 {
265 return C[i];
266 }
267
268 /// Allow the values to be set
269 void set_value(const unsigned& i,
270 const unsigned& j,
271 const unsigned& k,
272 const unsigned& l,
273 const double& value)
274 {
275 C[this->Index[i][j][k][l]] = value;
276 }
277 };
278
279
280} // namespace oomph
281#endif
cstr elem_len * i
Definition: cfortran.h:603
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.
virtual void set_value(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l, const double &value)
Allow the values to be set (virtual function that must be overloaded if values can be set directly.
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...
A general elasticity tensor that provides storage for all 21 independent components.
double independent_component(const unsigned &i) const
Overload the independent coefficient function.
void set_value(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l, const double &value)
Allow the values to be set.
An isotropic elasticity tensor defined in terms of Young's modulus and Poisson's ratio....
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.
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...
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....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...