elastoVariationalFunctionals.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2015-2019 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef FEM_DIFFOPS_ELASTOVARIATIONALFUNCTIONALS_HH
#define FEM_DIFFOPS_ELASTOVARIATIONALFUNCTIONALS_HH
 
#include "fem/diffops/materialLaws.hh"
 
namespace Kaskade 
{
  namespace Elastomechanics
  {
    
    template <class HyperelasticEnergy, class StrainTensor>
    class HyperelasticVariationalFunctional
    {
    public:
      using Scalar = typename HyperelasticEnergy::Scalar;
      static int const dim = HyperelasticEnergy::dim;
      using Tensor = Dune::FieldMatrix<Scalar,dim,dim>;
      
      template <typename... Args>
      HyperelasticVariationalFunctional(Args&&... args)
      : energy(std::forward<Args>(args)...)
      {
      }
      
      HyperelasticVariationalFunctional(HyperelasticEnergy const& energy_, StrainTensor const& strain_)
      : du(0), energy(energy_), strain(strain_)
      { 
      }
      
      void setLinearizationPoint(Tensor const& du0)
      {
        du = du0;
        strain = StrainTensor(du);
        energy.setLinearizationPoint(strain.d0());
      }
      
      Scalar d0() const
      {
        return energy.d0();
      }
      
      Dune::FieldVector<Scalar,dim> d1(VariationalArg<Scalar,dim> const& arg) const
      {
        // TODO: I have doubts that this is the most efficient implementation...
        Dune::FieldVector<Scalar, dim> ret;
        
        // step trough all spatial dimensions and test with vectorial test function 
        // with values only in that direction.
        for (int i=0; i<dim; ++i) {
          Dune::FieldMatrix<Scalar,dim,dim> dv(0);
          dv[i] = arg.derivative[0];
          
          Dune::FieldMatrix<Scalar,dim,dim> dE = strain.d1(dv);
          ret[i] = energy.d1(dE);
        }
                
        return ret;
      }
      
      Tensor d2(VariationalArg<Scalar,dim> const &arg1, VariationalArg<Scalar,dim> const &arg2) const 
      {
        // TODO: I have doubts that this is the most efficient implementation...
        Dune::FieldMatrix<Scalar, dim, dim> ret;
        
        // In a double loop, step trough all spatial dimensions and test with vectorial test functions 
        // with values only in that direction.
        // TODO: due to symmetry, can we skip half the computation?
        for (int i=0; i<dim; ++i)
        {
          Dune::FieldMatrix<Scalar,dim,dim> dv1(0);
          dv1[i] = arg1.derivative[0];
          Dune::FieldMatrix<Scalar,dim,dim> dE1 = strain.d1(dv1);
          
          for (int j=0; j<dim; ++j)
          {
            Dune::FieldMatrix<Scalar,dim,dim> dv2(0);
            dv2[j] = arg2.derivative[0];
            Dune::FieldMatrix<Scalar,dim,dim> dE2 = strain.d1(dv2);
            Dune::FieldMatrix<Scalar,dim,dim> ddE = strain.d2(dv1,dv2);
            
            // according to product rule...
            ret[i][j] = energy.d2(dE1,dE2) + energy.d1(ddE);
          }
        }
        return ret;
      }
      
      HyperelasticEnergy& getHyperelasticEnergy() 
      {
        return energy;
      }
 
      HyperelasticEnergy const& getHyperelasticEnergy() const
      {
        return energy;
      }
 
      Tensor const& getDu() const
      {
        return du;
      }
 
      StrainTensor const& getStrain() const
      {
        return strain;
      }
 
      Tensor stress() const { return energy.stress(); }
      
      Tensor cauchyStress() const 
      { 
        Tensor F = du; 
        for (int i=0; i<dim; ++i)
          F[i][i] += 1; // F = I + ux
        
        return F*stress()*transpose(F) / F.determinant(); 
      }
      
    protected:
      Tensor             du;
      HyperelasticEnergy energy;
      StrainTensor       strain;
    };
    
    // ---------------------------------------------------------------------------------------------------------
    
    template <int dim_, class Scalar_=double>
    class LameNavier
    : public HyperelasticVariationalFunctional<MaterialLaws::StVenantKirchhoff<dim_,Scalar_>,LinearizedGreenLagrangeTensor<Scalar_,dim_>>
    {
      using Energy = MaterialLaws::StVenantKirchhoff<dim_,Scalar_>;
      using Base = HyperelasticVariationalFunctional<Energy,LinearizedGreenLagrangeTensor<Scalar_,dim_>>;
      
      
    public:
      using Scalar = Scalar_;
      
    private:
      Scalar lambda, mu;
      
    public:
      static int const dim = dim_;
      
      LameNavier()
      : LameNavier(ElasticModulus(1,1)) {}
      
      LameNavier(ElasticModulus const& moduli)
      : Base(moduli), lambda(moduli.lame()), mu(moduli.shear())
      {
      }
 
      Dune::FieldVector<Scalar,dim> d1(VariationalArg<Scalar,dim> const& arg) const
      {
        Dune::FieldVector<Scalar, dim> ret = arg.derivative[0];
        ret *= lambda * trace(this->strain.d0());
        ret += 2*mu * (this->strain.d0()*arg.derivative[0]);
 
        return ret;
      }
 
      typename Base::Tensor d2(VariationalArg<Scalar,dim> const &arg1, VariationalArg<Scalar,dim> const &arg2) const
      {
        Dune::FieldMatrix<Scalar, dim, dim> ret;
 
        for (int i=0; i<dim; ++i) {
          double trace_ev = arg1.derivative[0][i];
 
          for (int j=0; j<dim; ++j) {
            double trace_ew = arg2.derivative[0][j];
            double ct = arg1.derivative[0][j] * arg2.derivative[0][i];
            if (j == i) {
              for (int k = 0; k < dim; ++k) {
                ct += arg1.derivative[0][k] * arg2.derivative[0][k];
              }
            }
            ret[i][j] += lambda*trace_ev*trace_ew + mu*ct;
          }
        }
        return ret;
      }
 
      void setElasticModuli(ElasticModulus const& modulus)
      {
        lambda = modulus.lame();
        mu = modulus.shear();
        this->energy = Energy(modulus);
      }
      
      ElasticModulus getElasticModuli()
      {
        return ElasticModulus(lambda,mu);
      }
      
      Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> strainToStressMatrix() const
      {
        // The implementation follows Braess, Finite Elemente, Chapter VI.3
        Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> C(0.0);
        
        if (dim==3)
        {
          C[0][0] = lambda+2*mu; C[0][1] = lambda; C[0][2] = lambda;
          C[1][0] = lambda; C[1][1] = lambda+2*mu; C[1][2] = lambda;
          C[2][0] = lambda; C[2][1] = lambda; C[2][2] = lambda+2*mu;
          for (int i=3; i<6; ++i)
            C[i][i] = 2*mu;
        }
        else if (dim==2)       
        {
          C[0][0] = lambda+2*mu; C[0][1] = lambda; 
          C[1][0] = lambda; C[1][1] = lambda+2*mu; 
          C[2][2] = 2*mu;
        }
        return C;
      }
      
      static Dune::FieldVector<Scalar,dim> mv(Dune::FieldVector<Scalar,dim*(dim+1)/2> const& A, Dune::FieldVector<Scalar,dim> const& x)
      {
        Dune::FieldVector<Scalar,dim> result; 
        
        if (dim==3)
        {
          result[0] = A[0]*x[0] + A[3]*x[1] + A[4]*x[2];
          result[1] = A[3]*x[0] + A[1]*x[1] + A[5]*x[2];
          result[2] = A[4]*x[0] + A[5]*x[1] + A[2]*x[2];
        }
        else if (dim==2)
        {
          result[0] = A[0]*x[0] + A[2]*x[1];
          result[1] = A[2]*x[0] + A[1]*x[1];
        }
        else if (dim==1)
        {
          abort();
        }
        
        return result;
      }
      
      Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*dim> displacementDerivativeToStrainMatrix() const
      {
        Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*dim> E(0.0);
        
        if (dim==3)
        {
          E[0][0] = 2; // diagonal elements of epsilon
          E[1][4] = 2;
          E[2][8] = 2;
          E[3][1] = 1; E[3][3] = 1; // off-diagonal elements of epsilon
          E[4][2] = 1; E[4][6] = 1;
          E[5][5] = 1; E[5][7] = 1;
        }
        else if (dim==2)
        {
          E[0][0] = 2;
          E[1][3] = 2;
          E[2][1] = 1; E[2][2] = 1;
        } else if (dim==1)
          E[0][0] = 2;
        
        E *= 0.5;
        return E;
      }
    };
    
    // ---------------------------------------------------------------------------------------------------------
    
    template <class Scalar, int dim>
    class OrthotropicLameNavier {
    public:
      OrthotropicLameNavier(Dune::FieldMatrix<Scalar,dim,dim> const& orth_, Dune::FieldMatrix<Scalar,dim,dim> const& mat_);
      
      OrthotropicLameNavier(Dune::FieldMatrix<Scalar,dim,dim> const& orth_, 
                            Dune::FieldMatrix<Scalar,dim,dim> const& mat1, Dune::FieldVector<Scalar,dim*(dim-1)/2> const& mat2);
      
      OrthotropicLameNavier(ElasticModulus const& p);
      
      OrthotropicLameNavier()
      : OrthotropicLameNavier(ElasticModulus()) {}
      
      Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> const& stiffnessMatrix() const { return stiffness; }
      
      void setLinearizationPoint(Dune::FieldMatrix<Scalar,dim,dim> const& du0_)
      {
        du0 = du0_;
      }
      
      Scalar d0() const;
      
      Dune::FieldVector<Scalar,dim> d1(VariationalArg<Scalar,dim> const& arg) const;
      
      Dune::FieldMatrix<Scalar,dim,dim> d2(VariationalArg<Scalar,dim> const &arg1, VariationalArg<Scalar,dim> const &arg2) const;
      
      void setLinear(bool lin){
                  linear = lin;
          }
      
      
    private:
      Dune::FieldMatrix<Scalar,dim,dim> orth; // the orthonormal matrix mapping material to global coordinates
      Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> stiffness;
      Dune::FieldMatrix<Scalar,dim,dim> du0; // the linearization point of the displacement derivative
      bool linear;
    };
    
    // ---------------------------------------------------------------------------------------------------------
 
    template <class HyperelasticEnergy, class StrainTensor>
    class FirstPiolaKirchhoffStress
    {
    public:
      using EnergyFunctional = HyperelasticVariationalFunctional<HyperelasticEnergy, StrainTensor>;
      using Scalar = typename HyperelasticEnergy::Scalar;
      using Tensor = typename EnergyFunctional::Tensor;
      static constexpr int dim = HyperelasticEnergy::dim;
 
      FirstPiolaKirchhoffStress(EnergyFunctional const& ef) : energy(ef) {}
 
      Tensor deformationGradient() const 
      {
        return deformationDerivative();
      }
      
      Tensor deformationDerivative() const 
      {
        Tensor F = energy.getDu();
        for (int i=0; i<dim; ++i) 
          F[i][i] += 1; // F = I + du
        return F;
      }
      
      Tensor d0() const
      {
        // T(du) = F(du) * Sigma(du)     Sigma is second Pioala Kirchhoff stress
        return deformationDerivative()*energy.stress();
      }
 
      Tensor d1(VariationalArg<Scalar, dim, dim> const& arg) const {
        // T'[dv] = F'[dv] * Sigma(du) + F * Sigma'[dv]
        // F'[dv] = dv
        // Sigma'[dv] = W''(E(du))[E'(du)[dv]]    W is hyperelastic energy (depending on strain E)
        return arg.derivative*energy.stress() + deformationGradient()*spkStressD1(arg);
      }
 
      Tensor d2(VariationalArg<Scalar, dim, dim> const& argv, VariationalArg<Scalar, dim, dim> const& argw) const {
        // T''[dv,dw] = dv*W''(E)[E'[dw]] + dw*W''(E)[E'[dv]] + F*Sigma''[dv,dw]
        // Sigma''[dv,dw] = W'''(E)[E'[dv],E'[dw]] + W''(E)[E''[dv,dw]]
        Tensor result{};
        result += argv.derivative*spkStressD1(argw);
        result += argw.derivative*spkStressD1(argv);
        //result += deformationGradient()*energyD3DeDf(energy.getStrain().d1(argv.derivative), energy.getStrain().d1(argw.derivative));
        result += deformationGradient()*energyD2De(energy.getStrain().d2(argv.derivative, argw.derivative));
        return result;
      }
 
    protected:
      // directional derivative of second Piola Kirchhoff stress in direction of displacement gradient
      Tensor spkStressD1(VariationalArg<Scalar, dim, dim> const& arg) const {
        // Sigma'[dv] = W''(E(du))[E'(du)[dv]]
        return energyD2De(energy.getStrain().d1(arg.derivative));
      }
 
      // directional derivative of first derivative of stored energy in direction of strain e
      Tensor energyD2De(Tensor const& e) const {
        // W''(E(du))[E'(du)[dv]]
        Tensor sigma;
 
        for (int i=0; i<dim; ++i)
          for (int j=0; j<=i; ++j)
          {
            Tensor e2(0);
            e2[i][j] = 1;
            sigma[i][j] = energy.getHyperelasticEnergy().d2(e, e2);
          }
 
        for (int i=0; i<dim-1; ++i)
          for (int j=i+1; j<dim; ++j)
            sigma[i][j] = sigma[j][i];
 
        return sigma;
      }
 
//      Tensor energyD3DeDf(Tensor const& e, Tensor const& f) {
//        // Todo: implement, only needed if hyperelastic energy has non vanishing third derivatives (with respect to strain)
//        return Tensor{};
//      }
 
    public:
      EnergyFunctional const& energy;
    };
 
    // ---------------------------------------------------------------------------------------------------------
 
  }
}
 
#endif