poroelasticity/elasticity_tensor.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2025 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Header file for objects representing the fourth-rank elasticity tensor
// for linear elasticity problems
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_POROELASTICITY_TENSOR_HEADER
#define OOMPH_POROELASTICITY_TENSOR_HEADER
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "generic/oomph_utilities.h"
 
namespace oomph
{
  //=====================================================================
  /// A base class that represents the fourth-rank elasticity tensor
  /// \f$E_{ijkl}\f$ defined such that
  /// \f[\tau_{ij} = E_{ijkl} e_{kl},\f]
  /// where \f$e_{ij}\f$ is the infinitessimal (Cauchy) strain tensor
  /// and \f$\tau_{ij}\f$ is the stress tensor. The symmetries of the
  /// tensor are such that
  /// \f[E_{ijkl} = E_{jikl} = E_{ijlk} = E_{klij}\f]
  /// and thus there are relatively few independent components. These
  /// symmetries are included in the definition of the object so that
  /// non-physical symmetries cannot be accidentally imposed.
  //=====================================================================
  class ElasticityTensor
  {
    /// Translation table from the four indices to the corresponding
    /// independent component
    static const unsigned Index[3][3][3][3];
 
  protected:
    /// Member function that returns the i-th independent component of the
    /// elasticity tensor
    virtual inline double independent_component(const unsigned& i) const
    {
      return 0.0;
    }
 
 
    /// Helper range checking function
    /// (Note that this only captures over-runs in 3D but
    /// errors are likely to be caught in evaluation of the
    /// stress and strain tensors anyway...)
    void range_check(const unsigned& i,
                     const unsigned& j,
                     const unsigned& k,
                     const unsigned& l) const
    {
      if ((i > 2) || (j > 2) || (k > 2) || (l > 2))
      {
        std::ostringstream error_message;
        if (i > 2)
        {
          error_message << "Range Error : Index 1 " << i
                        << " is not in the range (0,2)";
        }
        if (j > 2)
        {
          error_message << "Range Error : Index 2 " << j
                        << " is not in the range (0,2)";
        }
 
        if (k > 2)
        {
          error_message << "Range Error : Index 2 " << k
                        << " is not in the range (0,2)";
        }
 
        if (l > 2)
        {
          error_message << "Range Error : Index 4 " << l
                        << " is not in the range (0,2)";
        }
 
        // Throw the error
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
 
    /// Empty Constructor
    ElasticityTensor() {}
 
  public:
    /// Empty virtual Destructor
    virtual ~ElasticityTensor() {}
 
  public:
    /// Return the appropriate independent component
    /// via the index translation scheme (const version).
    double operator()(const unsigned& i,
                      const unsigned& j,
                      const unsigned& k,
                      const unsigned& l) const
    {
      // Range check
#ifdef PARANOID
      range_check(i, j, k, l);
#endif
      return independent_component(Index[i][j][k][l]);
    }
  };
 
 
  //===================================================================
  /// An isotropic elasticity tensor defined in terms of Young's modulus
  /// and Poisson's ratio. The elasticity tensor is assumed to be
  /// non-dimensionalised on some reference value for Young's modulus
  /// so the value provided to the constructor (if any) is to be
  /// interpreted as the ratio of the actual Young's modulus to the
  /// Young's modulus used to non-dimensionalise the stresses/tractions
  /// in the governing equations.
  //===================================================================
  class IsotropicElasticityTensor : public ElasticityTensor
  {
    // Storage for the independent components of the elasticity tensor
    double C[4];
 
    // Translation scheme between the 21 independent components of the general
    // elasticity tensor and the isotropic case
    static const unsigned StaticIndex[21];
 
  public:
    /// Constructor. Passing in the values of the Poisson's ratio
    /// and Young's modulus (interpreted as the ratio of the actual
    /// Young's modulus to the Young's modulus (or other reference stiffness)
    /// used to non-dimensionalise stresses and tractions in the governing
    /// equations).
    IsotropicElasticityTensor(const double& nu, const double& E)
      : ElasticityTensor()
    {
      // Set the three independent components
      C[0] = 0.0;
      double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
      double mu = E / (2.0 * (1.0 + nu));
      this->set_lame_coefficients(lambda, mu);
    }
 
    /// Constructor. Passing in the value of the Poisson's ratio.
    /// Stresses and tractions in the governing equations are assumed
    /// to have been non-dimensionalised on Young's modulus.
    IsotropicElasticityTensor(const double& nu) : ElasticityTensor()
    {
      // Set the three independent components
      C[0] = 0.0;
 
      // reference value
      double E = 1.0;
      double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
      double mu = E / (2.0 * (1.0 + nu));
      this->set_lame_coefficients(lambda, mu);
    }
 
    /// Constructur. Passing in the values of the two lame
    /// coefficients directly (interpreted as the ratios of these
    /// quantities to a reference stiffness used to non-dimensionalised
    IsotropicElasticityTensor(const Vector<double>& lame)
    {
      // Set the three independent components
      C[0] = 0.0;
      this->set_lame_coefficients(lame[0], lame[1]);
    }
 
 
    /// Overload the independent coefficient function
    inline double independent_component(const unsigned& i) const
    {
      return C[StaticIndex[i]];
    }
 
 
  private:
    // Set the values of the lame coefficients
    void set_lame_coefficients(const double& lambda, const double& mu)
    {
      C[1] = lambda + 2.0 * mu;
      C[2] = lambda;
      C[3] = mu;
    }
  };
 
 
  //===================================================================
  /// An isotropic elasticity tensor defined in terms of Young's modulus
  /// and Poisson's ratio. The elasticity tensor is assumed to be
  /// non-dimensionalised on some reference value for Young's modulus
  /// so the value provided to the constructor (if any) is to be
  /// interpreted as the ratio of the actual Young's modulus to the
  /// Young's modulus used to non-dimensionalise the stresses/tractions
  /// in the governing equations.
  //===================================================================
  class DeviatoricIsotropicElasticityTensor : public ElasticityTensor
  {
    // Storage for the independent components of the elasticity tensor
    double C[3];
 
    // Storage for the first Lame parameter
    double Lambda;
 
    // Storage for the second Lame parameter
    double Mu;
 
    // Translation scheme between the 21 independent components of the general
    // elasticity tensor and the isotropic case
    static const unsigned StaticIndex[21];
 
  public:
    /// Constructor. For use with incompressibility. Requires no
    /// parameters since Poisson's ratio is fixed at 0.5 and lambda is set to a
    /// dummy value of 0 (since it would be infinite)
    DeviatoricIsotropicElasticityTensor() : ElasticityTensor()
    {
      C[0] = 0.0;
      double E = 1.0;
      double nu = 0.5;
      double lambda = 0.0;
      double mu = E / (2.0 * (1.0 + nu));
      this->set_lame_coefficients(lambda, mu);
    }
 
    /// Constructor. Passing in the values of the Poisson's ratio
    /// and Young's modulus (interpreted as the ratio of the actual
    /// Young's modulus to the Young's modulus (or other reference stiffness)
    /// used to non-dimensionalise stresses and tractions in the governing
    /// equations).
    DeviatoricIsotropicElasticityTensor(const double& nu, const double& E)
      : ElasticityTensor()
    {
      // Set the three independent components
      C[0] = 0.0;
      double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
      double mu = E / (2.0 * (1.0 + nu));
      this->set_lame_coefficients(lambda, mu);
    }
 
    /// Constructor. Passing in the value of the Poisson's ratio.
    /// Stresses and tractions in the governing equations are assumed
    /// to have been non-dimensionalised on Young's modulus.
    DeviatoricIsotropicElasticityTensor(const double& nu) : ElasticityTensor()
    {
      // Set the three independent components
      C[0] = 0.0;
 
      // reference value
      double E = 1.0;
      double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
      double mu = E / (2.0 * (1.0 + nu));
      this->set_lame_coefficients(lambda, mu);
    }
 
    /// Constructur. Passing in the values of the two lame
    /// coefficients directly (interpreted as the ratios of these
    /// quantities to a reference stiffness used to non-dimensionalised
    DeviatoricIsotropicElasticityTensor(const Vector<double>& lame)
    {
      // Set the three independent components
      C[0] = 0.0;
      this->set_lame_coefficients(lame[0], lame[1]);
    }
 
 
    /// Overload the independent coefficient function
    inline double independent_component(const unsigned& i) const
    {
      return C[StaticIndex[i]];
    }
 
    /// Accessor function for the first lame parameter
    const double& lambda() const
    {
      return Lambda;
    }
 
    /// Accessor function for the second lame parameter
    const double& mu() const
    {
      return Mu;
    }
 
  private:
    // Set the values of the lame coefficients
    void set_lame_coefficients(const double& lambda, const double& mu)
    {
      C[1] = 2.0 * mu;
      C[2] = mu;
 
      Lambda = lambda;
      Mu = mu;
    }
  };
 
} // namespace oomph
#endif