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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for objects representing the fourth-rank elasticity tensor
31 //for linear elasticity problems
32 
33 //Include guards to prevent multiple inclusion of the header
34 #ifndef OOMPH_POROELASTICITY_TENSOR_HEADER
35 #define OOMPH_POROELASTICITY_TENSOR_HEADER
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39  #include <oomph-lib-config.h>
40 #endif
41 
42 #include "../generic/oomph_utilities.h"
43 
44 namespace oomph
45 {
46  //=====================================================================
47  /// A base class that represents the fourth-rank elasticity tensor
48  /// \f$E_{ijkl}\f$ defined such that
49  /// \f[\tau_{ij} = E_{ijkl} e_{kl},\f]
50  /// where \f$e_{ij}\f$ is the infinitessimal (Cauchy) strain tensor
51  /// and \f$\tau_{ij}\f$ is the stress tensor. The symmetries of the
52  /// tensor are such that
53  /// \f[E_{ijkl} = E_{jikl} = E_{ijlk} = E_{klij}\f]
54  /// and thus there are relatively few independent components. These
55  /// symmetries are included in the definition of the object so that
56  /// non-physical symmetries cannot be accidentally imposed.
57  //=====================================================================
58 class ElasticityTensor
59  {
60 
61  ///\short Translation table from the four indices to the corresponding
62  ///independent component
63  static const unsigned Index[3][3][3][3];
64 
65  protected:
66 
67  ///Member function that returns the i-th independent component of the
68  ///elasticity tensor
69  virtual inline double independent_component(const unsigned &i) const
70  {return 0.0;}
71 
72 
73  ///\short Helper range checking function
74  /// (Note that this only captures over-runs in 3D but
75  /// errors are likely to be caught in evaluation of the
76  /// stress and strain tensors anyway...)
77  void range_check(const unsigned &i, const unsigned &j,
78  const unsigned &k, const unsigned &l) const
79  {
80  if((i > 2) || (j > 2) || (k> 2) || (l>2))
81  {
82  std::ostringstream error_message;
83  if(i > 2)
84  {
85  error_message << "Range Error : Index 1 " << i
86  << " is not in the range (0,2)";
87  }
88  if(j > 2)
89  {
90  error_message << "Range Error : Index 2 " << j
91  << " is not in the range (0,2)";
92  }
93 
94  if(k > 2)
95  {
96  error_message << "Range Error : Index 2 " << k
97  << " is not in the range (0,2)";
98  }
99 
100  if(l > 2)
101  {
102  error_message << "Range Error : Index 4 " << l
103  << " is not in the range (0,2)";
104  }
105 
106  //Throw the error
107  throw OomphLibError(error_message.str(),
108  OOMPH_CURRENT_FUNCTION,
109  OOMPH_EXCEPTION_LOCATION);
110  }
111  }
112 
113 
114  ///Empty Constructor
116 
117  public:
118 
119  ///Empty virtual Destructor
120  virtual ~ElasticityTensor() {}
121 
122  public:
123 
124  ///\short Return the appropriate independent component
125  ///via the index translation scheme (const version).
126  double operator()(const unsigned &i, const unsigned &j,
127  const unsigned &k, const unsigned &l) const
128  {
129  //Range check
130 #ifdef PARANOID
131  range_check(i,j,k,l);
132 #endif
133  return independent_component(Index[i][j][k][l]);
134  }
135  };
136 
137 
138 //===================================================================
139 /// An isotropic elasticity tensor defined in terms of Young's modulus
140 /// and Poisson's ratio. The elasticity tensor is assumed to be
141 /// non-dimensionalised on some reference value for Young's modulus
142 /// so the value provided to the constructor (if any) is to be
143 /// interpreted as the ratio of the actual Young's modulus to the
144 /// Young's modulus used to non-dimensionalise the stresses/tractions
145 /// in the governing equations.
146 //===================================================================
147  class IsotropicElasticityTensor : public ElasticityTensor
148  {
149  //Storage for the independent components of the elasticity tensor
150  double C[4];
151 
152  //Translation scheme between the 21 independent components of the general
153  //elasticity tensor and the isotropic case
154  static const unsigned StaticIndex[21];
155 
156  public:
157 
158  /// \short Constructor. Passing in the values of the Poisson's ratio
159  /// and Young's modulus (interpreted as the ratio of the actual
160  /// Young's modulus to the Young's modulus (or other reference stiffness)
161  /// used to non-dimensionalise stresses and tractions in the governing
162  /// equations).
163  IsotropicElasticityTensor(const double &nu,
164  const double &E) : ElasticityTensor()
165  {
166  //Set the three independent components
167  C[0] = 0.0;
168  double lambda=E*nu/((1.0+nu)*(1.0-2.0*nu));
169  double mu=E/(2.0*(1.0+nu));
170  this->set_lame_coefficients(lambda,mu);
171  }
172 
173  /// \short Constructor. Passing in the value of the Poisson's ratio.
174  /// Stresses and tractions in the governing equations are assumed
175  /// to have been non-dimensionalised on Young's modulus.
177  {
178  //Set the three independent components
179  C[0] = 0.0;
180 
181  // reference value
182  double E=1.0;
183  double lambda=E*nu/((1.0+nu)*(1.0-2.0*nu));
184  double mu=E/(2.0*(1.0+nu));
185  this->set_lame_coefficients(lambda,mu);
186  }
187 
188  /// \short Constructur. Passing in the values of the two lame
189  /// coefficients directly (interpreted as the ratios of these
190  /// quantities to a reference stiffness used to non-dimensionalised
192  {
193  //Set the three independent components
194  C[0] = 0.0;
195  this->set_lame_coefficients(lame[0],lame[1]);
196  }
197 
198 
199  ///Overload the independent coefficient function
200  inline double independent_component(const unsigned &i) const
201  {return C[StaticIndex[i]];}
202 
203 
204  private:
205 
206  //Set the values of the lame coefficients
207  void set_lame_coefficients(const double &lambda, const double &mu)
208  {
209  C[1] = lambda + 2.0*mu;
210  C[2] = lambda;
211  C[3] = mu;
212  }
213 
214  };
215 
216 
217 //===================================================================
218 /// An isotropic elasticity tensor defined in terms of Young's modulus
219 /// and Poisson's ratio. The elasticity tensor is assumed to be
220 /// non-dimensionalised on some reference value for Young's modulus
221 /// so the value provided to the constructor (if any) is to be
222 /// interpreted as the ratio of the actual Young's modulus to the
223 /// Young's modulus used to non-dimensionalise the stresses/tractions
224 /// in the governing equations.
225 //===================================================================
227  {
228  //Storage for the independent components of the elasticity tensor
229  double C[3];
230 
231  //Storage for the first Lame parameter
232  double Lambda;
233 
234  //Storage for the second Lame parameter
235  double Mu;
236 
237  //Translation scheme between the 21 independent components of the general
238  //elasticity tensor and the isotropic case
239  static const unsigned StaticIndex[21];
240 
241  public:
242 
243  /// \short Constructor. For use with incompressibility. Requires no
244  /// parameters since Poisson's ratio is fixed at 0.5 and lambda is set to a
245  /// dummy value of 0 (since it would be infinite)
247  {
248  C[0] = 0.0;
249  double E=1.0;
250  double nu=0.5;
251  double lambda=0.0;
252  double mu=E/(2.0*(1.0+nu));
253  this->set_lame_coefficients(lambda,mu);
254  }
255 
256  /// \short Constructor. Passing in the values of the Poisson's ratio
257  /// and Young's modulus (interpreted as the ratio of the actual
258  /// Young's modulus to the Young's modulus (or other reference stiffness)
259  /// used to non-dimensionalise stresses and tractions in the governing
260  /// equations).
262  const double &E) : ElasticityTensor()
263  {
264  //Set the three independent components
265  C[0] = 0.0;
266  double lambda=E*nu/((1.0+nu)*(1.0-2.0*nu));
267  double mu=E/(2.0*(1.0+nu));
268  this->set_lame_coefficients(lambda,mu);
269  }
270 
271  /// \short Constructor. Passing in the value of the Poisson's ratio.
272  /// Stresses and tractions in the governing equations are assumed
273  /// to have been non-dimensionalised on Young's modulus.
275  {
276  //Set the three independent components
277  C[0] = 0.0;
278 
279  // reference value
280  double E=1.0;
281  double lambda=E*nu/((1.0+nu)*(1.0-2.0*nu));
282  double mu=E/(2.0*(1.0+nu));
283  this->set_lame_coefficients(lambda,mu);
284  }
285 
286  /// \short Constructur. Passing in the values of the two lame
287  /// coefficients directly (interpreted as the ratios of these
288  /// quantities to a reference stiffness used to non-dimensionalised
290  {
291  //Set the three independent components
292  C[0] = 0.0;
293  this->set_lame_coefficients(lame[0],lame[1]);
294  }
295 
296 
297  ///Overload the independent coefficient function
298  inline double independent_component(const unsigned &i) const
299  {return C[StaticIndex[i]];}
300 
301  ///Accessor function for the first lame parameter
302  const double &lambda() const
303  {
304  return Lambda;
305  }
306 
307  ///Accessor function for the second lame parameter
308  const double &mu() const
309  {
310  return Mu;
311  }
312 
313  private:
314 
315  //Set the values of the lame coefficients
316  void set_lame_coefficients(const double &lambda, const double &mu)
317  {
318  C[1] = 2.0*mu;
319  C[2] = mu;
320 
321  Lambda=lambda;
322  Mu=mu;
323  }
324  };
325 
326 }
327 #endif
double independent_component(const unsigned &i) const
Overload the independent coefficient function.
cstr elem_len * i
Definition: cfortran.h:607
static const unsigned StaticIndex[21]
Translation scheme for the isotropic elasticity tensor.
static const unsigned StaticIndex[21]
Translation scheme for the deviatoric isotropic elasticity tensor.
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...
DeviatoricIsotropicElasticityTensor()
Constructor. For use with incompressibility. Requires no parameters since Poisson's ratio is fixed at...
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...
static const unsigned Index[3][3][3][3]
Translation table from the four indices to the corresponding independent component.
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)...
const double & mu() const
Accessor function for the second lame parameter.
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...
virtual double independent_component(const unsigned &i) const
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)
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)
IsotropicElasticityTensor(const double &nu)
Constructor. Passing in the value of the Poisson's ratio. Stresses and tractions in the governing equ...