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// 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_ELASTICITY_TENSOR_HEADER
35 #define OOMPH_ELASTICITY_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  //=====================================================================
59  {
60 
61  protected:
62 
63  ///\short Translation table from the four indices to the corresponding
64  ///independent component
65  static const unsigned Index[3][3][3][3];
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  ///\short 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, const unsigned &j,
77  const unsigned &k, const unsigned &l) const
78  {
79  if((i > 2) || (j > 2) || (k> 2) || (l>2))
80  {
81  std::ostringstream error_message;
82  if(i > 2)
83  {
84  error_message << "Range Error : Index 1 " << i
85  << " is not in the range (0,2)";
86  }
87  if(j > 2)
88  {
89  error_message << "Range Error : Index 2 " << j
90  << " is not in the range (0,2)";
91  }
92 
93  if(k > 2)
94  {
95  error_message << "Range Error : Index 2 " << k
96  << " is not in the range (0,2)";
97  }
98 
99  if(l > 2)
100  {
101  error_message << "Range Error : Index 4 " << l
102  << " is not in the range (0,2)";
103  }
104 
105  //Throw the error
106  throw OomphLibError(error_message.str(),
107  OOMPH_CURRENT_FUNCTION,
108  OOMPH_EXCEPTION_LOCATION);
109  }
110  }
111 
112  ///Empty Constructor
114 
115  public:
116 
117  ///Empty virtual Destructor
118  virtual ~ElasticityTensor() {}
119 
120  public:
121 
122  ///\short Return the appropriate independent component
123  ///via the index translation scheme (const version).
124  double operator()(const unsigned &i, const unsigned &j,
125  const unsigned &k, 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  /// \short Allow the values to be set (virtual function that must be
135  /// overloaded if values can be set directly
136  virtual void set_value(const unsigned &i, const unsigned &j,
137  const unsigned &k, 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 //===================================================================
152 /// An isotropic elasticity tensor defined in terms of Young's modulus
153 /// and Poisson's ratio. The elasticity tensor is assumed to be
154 /// non-dimensionalised on some reference value for Young's modulus
155 /// so the value provided to the constructor (if any) is to be
156 /// interpreted as the ratio of the actual Young's modulus to the
157 /// Young's modulus used to non-dimensionalise the stresses/tractions
158 /// in the governing equations.
159 //===================================================================
161  {
162  //Storage for the independent components of the elasticity tensor
163  double C[4];
164 
165  //Translation scheme between the 21 independent components of the general
166  //elasticity tensor and the isotropic case
167  static const unsigned StaticIndex[21];
168 
169  public:
170 
171  /// \short Constructor. Passing in the values of the Poisson's ratio
172  /// and Young's modulus (interpreted as the ratio of the actual
173  /// Young's modulus to the Young's modulus (or other reference stiffness)
174  /// used to non-dimensionalise stresses and tractions in the governing
175  /// equations).
176  IsotropicElasticityTensor(const double &nu,
177  const double &E) : ElasticityTensor()
178  {
179  //Set the three indepdent components
180  C[0] = 0.0;
181  double lambda=E*nu/((1.0+nu)*(1.0-2.0*nu));
182  double mu=E/(2.0*(1.0+nu));
183  this->set_lame_coefficients(lambda,mu);
184  }
185 
186  /// \short Constructor. Passing in the value of the Poisson's ratio.
187  /// Stresses and tractions in the governing equations are assumed
188  /// to have been non-dimensionalised on Young's modulus.
190  {
191  //Set the three indepdent components
192  C[0] = 0.0;
193 
194  // reference value
195  double E=1.0;
196  double lambda=E*nu/((1.0+nu)*(1.0-2.0*nu));
197  double mu=E/(2.0*(1.0+nu));
198  this->set_lame_coefficients(lambda,mu);
199  }
200 
201  /// \short Constructur. Passing in the values of the two lame
202  /// coefficients directly (interpreted as the ratios of these
203  /// quantities to a reference stiffness used to non-dimensionalised
205  {
206  //Set the three independent componens
207  C[0] = 0.0;
208  this->set_lame_coefficients(lame[0],lame[1]);
209  }
210 
211  /// \short Update parameters: Specify values of the Poisson's ratio
212  /// and (optionally) Young's modulus (interpreted as the ratio of the actual
213  /// Young's modulus to the Young's modulus (or other reference stiffness)
214  /// used to non-dimensionalise stresses and tractions in the governing
215  /// equations).
216  void update_constitutive_parameters(const double &nu,
217  const double &E=1.0)
218  {
219  //Set the three indepdent components
220  C[0] = 0.0;
221  double lambda=E*nu/((1.0+nu)*(1.0-2.0*nu));
222  double mu=E/(2.0*(1.0+nu));
223  this->set_lame_coefficients(lambda,mu);
224  }
225 
226 
227  ///Overload the independent coefficient function
228  inline double independent_component(const unsigned &i) const
229  {return C[StaticIndex[i]];}
230 
231 
232  private:
233 
234  //Set the values of the lame coefficients
235  void set_lame_coefficients(const double &lambda, const double &mu)
236  {
237  C[1] = lambda + 2.0*mu;
238  C[2] = lambda;
239  C[3] = mu;
240  }
241 
242  };
243 
244 
245 
246 //========================================================================
247 /// A general elasticity tensor that provides storage for all 21 independent
248 /// components
249 //===================================================================
251 {
252  //Storage for the independent components of the elasticity tensor
253  double C[21];
254 
255  public:
256 
257  /// \short Empy Constructor.
259  {
260  //Initialise all independent components to zero
261  for(unsigned i=0;i<21;i++) {C[i] = 0.0;}
262  }
263 
264  ///Overload the independent coefficient function
265  inline double independent_component(const unsigned &i) const
266  {return C[i];}
267 
268  /// \short Allow the values to be set
269  void set_value(const unsigned &i, const unsigned &j,
270  const unsigned &k, const unsigned &l, const double &value)
271  {
272  C[this->Index[i][j][k][l]] = value;
273  }
274 
275 
276  };
277 
278 
279 
280 
281 
282 }
283 #endif
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...
cstr elem_len * i
Definition: cfortran.h:607
static const unsigned StaticIndex[21]
Translation scheme for the 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...
void set_value(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l, const double &value)
Allow the values to be set.
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...
double independent_component(const unsigned &i) const
Overload the independent coefficient function.
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)...
double independent_component(const unsigned &i) const
Overload the independent coefficient function.
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)
IsotropicElasticityTensor(const double &nu)
Constructor. Passing in the value of the Poisson's ratio. Stresses and tractions in the governing equ...