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