constitutive_laws.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 ConstitutiveLaw objects that will be used in all
// elasticity-type elements
 
#ifndef OOMPH_CONSTITUTIVE_LAWS_HEADER
#define OOMPH_CONSTITUTIVE_LAWS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB includes
#include "../generic/oomph_utilities.h"
#include "../generic/matrices.h"
 
namespace oomph
{
  //=====================================================================
  //====================================================================
  class StrainEnergyFunction
  {
  public:
    StrainEnergyFunction() {}
 
    virtual ~StrainEnergyFunction() {}
 
 
    virtual double W(const DenseMatrix<double>& gamma)
    {
      std::string error_message =
        "The strain-energy function as a function of the strain-tensor,\n";
      error_message +=
        "gamma, is not implemented for this strain energy function.\n";
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      return 0.0;
    }
 
 
    virtual double W(const Vector<double>& I)
    {
      std::string error_message =
        "The strain-energy function as a function of the strain\n ";
      error_message +=
        "invariants, I1, I2, I3, is not implemented for this strain\n ";
      error_message += "energy function\n";
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      return 0.0;
    }
 
 
    virtual void derivative(const DenseMatrix<double>& gamma,
                            DenseMatrix<double>& dWdgamma)
    {
      throw OomphLibError(
        "Sorry, the FD setup of dW/dgamma hasn't been implemented yet",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
 
    virtual void derivatives(Vector<double>& I, Vector<double>& dWdI)
    {
      // Calculate the derivatives of the strain-energy-function wrt the strain
      // invariants
      double FD_Jstep = 1.0e-8; // Usual comments about global stuff
      double energy = W(I);
 
      // Loop over the strain invariants
      for (unsigned i = 0; i < 3; i++)
      {
        // Store old value
        double I_prev = I[i];
        // Increase ith strain invariant
        I[i] += FD_Jstep;
        // Get the new value of the strain energy
        double energy_new = W(I);
        // Calculate the value of the derivative
        dWdI[i] = (energy_new - energy) / FD_Jstep;
        // Reset value of ith strain invariant
        I[i] = I_prev;
      }
    }
 
    virtual bool requires_incompressibility_constraint() = 0;
  };
 
 
 
 
  //=====================================================================
  //====================================================================
  class MooneyRivlin : public StrainEnergyFunction
  {
  public:
    MooneyRivlin(double* c1_pt, double* c2_pt)
      : StrainEnergyFunction(), C1_pt(c1_pt), C2_pt(c2_pt)
    {
    }
 
    virtual ~MooneyRivlin() {}
 
    double W(const DenseMatrix<double>& gamma)
    {
      return StrainEnergyFunction::W(gamma);
    }
 
    double W(const Vector<double>& I)
    {
      return (*C1_pt) * (I[0] - 3.0) + (*C2_pt) * (I[1] - 3.0);
    }
 
 
    void derivatives(Vector<double>& I, Vector<double>& dWdI)
    {
      dWdI[0] = (*C1_pt);
      dWdI[1] = (*C2_pt);
      dWdI[2] = 0.0;
    }
 
    bool requires_incompressibility_constraint()
    {
      return true;
    }
 
 
  private:
    double* C1_pt;
 
    double* C2_pt;
  };
 
 
 
 
  //=====================================================================
  //====================================================================
  class GeneralisedMooneyRivlin : public StrainEnergyFunction
  {
  public:
    GeneralisedMooneyRivlin(double* nu_pt, double* c1_pt)
      : StrainEnergyFunction(),
        Nu_pt(nu_pt),
        C1_pt(c1_pt),
        E_pt(new double(1.0)),
        Must_delete_e(true)
    {
    }
 
    GeneralisedMooneyRivlin(double* nu_pt, double* c1_pt, double* e_pt)
      : StrainEnergyFunction(),
        Nu_pt(nu_pt),
        C1_pt(c1_pt),
        E_pt(e_pt),
        Must_delete_e(false)
    {
    }
 
 
    virtual ~GeneralisedMooneyRivlin()
    {
      if (Must_delete_e) delete E_pt;
    }
 
    double W(const DenseMatrix<double>& gamma)
    {
      return StrainEnergyFunction::W(gamma);
    }
 
 
    double W(const Vector<double>& I)
    {
      double G = (*E_pt) / (2.0 * (1.0 + (*Nu_pt)));
      return 0.5 * ((*C1_pt) * (I[0] - 3.0) + (G - (*C1_pt)) * (I[1] - 3.0) +
                    ((*C1_pt) - 2.0 * G) * (I[2] - 1.0) +
                    (1.0 - (*Nu_pt)) * G * (I[2] - 1.0) * (I[2] - 1.0) /
                      (2.0 * (1.0 - 2.0 * (*Nu_pt))));
    }
 
 
    void derivatives(Vector<double>& I, Vector<double>& dWdI)
    {
      double G = (*E_pt) / (2.0 * (1.0 + (*Nu_pt)));
      dWdI[0] = 0.5 * (*C1_pt);
      dWdI[1] = 0.5 * (G - (*C1_pt));
      dWdI[2] = 0.5 * ((*C1_pt) - 2.0 * G +
                       2.0 * (1.0 - (*Nu_pt)) * G * (I[2] - 1.0) /
                         (2.0 * (1.0 - 2.0 * (*Nu_pt))));
    }
 
 
    bool requires_incompressibility_constraint()
    {
      return false;
    }
 
  private:
    double* Nu_pt;
 
    double* C1_pt;
 
    double* E_pt;
 
    bool Must_delete_e;
  };
 
 
 
 
  //===========================================================================
  //==========================================================================
  class ConstitutiveLaw
  {
  protected:
    bool is_matrix_square(const DenseMatrix<double>& M);
 
    bool are_matrices_of_equal_dimensions(const DenseMatrix<double>& M1,
                                          const DenseMatrix<double>& M2);
 
    void error_checking_in_input(const DenseMatrix<double>& g,
                                 const DenseMatrix<double>& G,
                                 DenseMatrix<double>& sigma);
 
    double calculate_contravariant(const DenseMatrix<double>& Gcov,
                                   DenseMatrix<double>& Gcontra);
 
    void calculate_d_contravariant_dG(const DenseMatrix<double>& Gcov,
                                      RankFourTensor<double>& dGcontra_dG,
                                      DenseMatrix<double>& d_detG_dG);
 
 
  public:
    ConstitutiveLaw() {}
 
 
    virtual ~ConstitutiveLaw() {}
 
 
    virtual void calculate_second_piola_kirchhoff_stress(
      const DenseMatrix<double>& g,
      const DenseMatrix<double>& G,
      DenseMatrix<double>& sigma) = 0;
 
    virtual void calculate_d_second_piola_kirchhoff_stress_dG(
      const DenseMatrix<double>& g,
      const DenseMatrix<double>& G,
      const DenseMatrix<double>& sigma,
      RankFourTensor<double>& d_sigma_dG,
      const bool& symmetrize_tensor = true);
 
 
    virtual void calculate_second_piola_kirchhoff_stress(
      const DenseMatrix<double>& g,
      const DenseMatrix<double>& G,
      DenseMatrix<double>& sigma_dev,
      DenseMatrix<double>& G_contra,
      double& Gdet)
    {
      throw OomphLibError(
        "Incompressible formulation not implemented for this constitutive law",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void calculate_d_second_piola_kirchhoff_stress_dG(
      const DenseMatrix<double>& g,
      const DenseMatrix<double>& G,
      const DenseMatrix<double>& sigma,
      const double& detG,
      const double& interpolated_solid_p,
      RankFourTensor<double>& d_sigma_dG,
      DenseMatrix<double>& d_detG_dG,
      const bool& symmetrize_tensor = true);
 
 
    virtual void calculate_second_piola_kirchhoff_stress(
      const DenseMatrix<double>& g,
      const DenseMatrix<double>& G,
      DenseMatrix<double>& sigma_dev,
      DenseMatrix<double>& Gcontra,
      double& gen_dil,
      double& inv_kappa)
    {
      throw OomphLibError(
        "Near-incompressible formulation not implemented for constitutive law",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void calculate_d_second_piola_kirchhoff_stress_dG(
      const DenseMatrix<double>& g,
      const DenseMatrix<double>& G,
      const DenseMatrix<double>& sigma,
      const double& gen_dil,
      const double& inv_kappa,
      const double& interpolated_solid_p,
      RankFourTensor<double>& d_sigma_dG,
      DenseMatrix<double>& d_gen_dil_dG,
      const bool& symmetrize_tensor = true);
 
 
    virtual bool requires_incompressibility_constraint() = 0;
  };
 
 
 
 
  //========================================================================
  //=========================================================================
  class GeneralisedHookean : public ConstitutiveLaw
  {
  public:
    GeneralisedHookean(double* nu_pt, double* e_pt)
      : ConstitutiveLaw(), Nu_pt(nu_pt), E_pt(e_pt), Must_delete_e(false)
    {
    }
 
    GeneralisedHookean(double* nu_pt)
      : ConstitutiveLaw(),
        Nu_pt(nu_pt),
        E_pt(new double(1.0)),
        Must_delete_e(true)
    {
    }
 
 
    virtual ~GeneralisedHookean()
    {
      if (Must_delete_e) delete E_pt;
    }
 
    void calculate_second_piola_kirchhoff_stress(const DenseMatrix<double>& g,
                                                 const DenseMatrix<double>& G,
                                                 DenseMatrix<double>& sigma);
 
 
    void calculate_second_piola_kirchhoff_stress(const DenseMatrix<double>& g,
                                                 const DenseMatrix<double>& G,
                                                 DenseMatrix<double>& sigma_dev,
                                                 DenseMatrix<double>& G_contra,
                                                 double& Gdet);
 
 
    void calculate_second_piola_kirchhoff_stress(const DenseMatrix<double>& g,
                                                 const DenseMatrix<double>& G,
                                                 DenseMatrix<double>& sigma_dev,
                                                 DenseMatrix<double>& Gcontra,
                                                 double& gen_dil,
                                                 double& inv_kappa);
 
 
    bool requires_incompressibility_constraint()
    {
      return false;
    }
 
  private:
    double* Nu_pt;
 
    double* E_pt;
 
    bool Must_delete_e;
  };
 
 
 
 
  //=====================================================================
  //=====================================================================
  class IsotropicStrainEnergyFunctionConstitutiveLaw : public ConstitutiveLaw
  {
  private:
    StrainEnergyFunction* Strain_energy_function_pt;
 
  public:
    IsotropicStrainEnergyFunctionConstitutiveLaw(
      StrainEnergyFunction* const& strain_energy_function_pt)
      : ConstitutiveLaw(), Strain_energy_function_pt(strain_energy_function_pt)
    {
    }
 
    void calculate_second_piola_kirchhoff_stress(const DenseMatrix<double>& g,
                                                 const DenseMatrix<double>& G,
                                                 DenseMatrix<double>& sigma);
 
 
    void calculate_second_piola_kirchhoff_stress(const DenseMatrix<double>& g,
                                                 const DenseMatrix<double>& G,
                                                 DenseMatrix<double>& sigma_dev,
                                                 DenseMatrix<double>& G_contra,
                                                 double& Gdet);
 
 
    void calculate_second_piola_kirchhoff_stress(const DenseMatrix<double>& g,
                                                 const DenseMatrix<double>& G,
                                                 DenseMatrix<double>& sigma_dev,
                                                 DenseMatrix<double>& Gcontra,
                                                 double& gen_dil,
                                                 double& inv_kappa);
 
 
    bool requires_incompressibility_constraint()
    {
      return Strain_energy_function_pt->requires_incompressibility_constraint();
    }
  };
 
} // namespace oomph
 
#endif