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-2022 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 generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "../generic/oomph_utilities.h"
 
namespace oomph
{
  //=====================================================================
  //=====================================================================
  class ElasticityTensor
  {
    static const unsigned Index[3][3][3][3];
 
  protected:
    virtual inline double independent_component(const unsigned& i) const
    {
      return 0.0;
    }
 
 
    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);
      }
    }
 
 
    ElasticityTensor() {}
 
  public:
    virtual ~ElasticityTensor() {}
 
  public:
    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]);
    }
  };
 
 
  //===================================================================
  //===================================================================
  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:
    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);
    }
 
    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);
    }
 
    IsotropicElasticityTensor(const Vector<double>& lame)
    {
      // Set the three independent components
      C[0] = 0.0;
      this->set_lame_coefficients(lame[0], lame[1]);
    }
 
 
    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;
    }
  };
 
 
  //===================================================================
  //===================================================================
  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:
    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);
    }
 
    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);
    }
 
    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);
    }
 
    DeviatoricIsotropicElasticityTensor(const Vector<double>& lame)
    {
      // Set the three independent components
      C[0] = 0.0;
      this->set_lame_coefficients(lame[0], lame[1]);
    }
 
 
    inline double independent_component(const unsigned& i) const
    {
      return C[StaticIndex[i]];
    }
 
    const double& lambda() const
    {
      return Lambda;
    }
 
    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