materialLaws.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_MATERIALLAWS_HH
#define FEM_DIFFOPS_MATERIALLAWS_HH
 
#include <cassert>
 
#include "dune/common/config.h"
#include "dune/common/fmatrix.hh"
 
#include "fem/diffops/elasto.hh"
#include "linalg/determinant.hh"
 
namespace Kaskade
{
  namespace Elastomechanics
  {
    namespace MaterialLaws 
    {
      template <int dim_, class Scalar_, class MaterialLaw>
      class MaterialLawBase
      {
      public:
        using Scalar = Scalar_;
        static int const dim = dim_;
        using Tensor = Dune::FieldMatrix<Scalar,dim,dim>; 
        using VoigtTensor = Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2>;
        
        Tensor stress() const
        {
          Tensor sigma;
          
          for (int i=0; i<dim; ++i)
            for (int j=0; j<=i; ++j)
            {
              Tensor E(0); 
              E[i][j] = 1;
              sigma[i][j] = static_cast<MaterialLaw const*>(this)->d1(E);
            }
            
          for (int i=0; i<dim-1; ++i)
            for (int j=i+1; j<dim; ++j)
              sigma[i][j] = sigma[j][i];
            
          return sigma;
        }
        
        VoigtTensor strainToStressMatrix() const
        {
          int const n = dim*(dim+1)/2;
          Dune::FieldMatrix<Scalar,n,n> C;
          
          // For the Voigt notation it holds that x^T C y = W''[unpack(x),unpack(y)]. 
          // Due to symmetry of C, we only compute the lower half directly...
          for (int i=0; i<n; ++i)
          {
            Dune::FieldVector<Scalar,n> x;
            x[i] = 1;
            auto X = unpack(x);
            
            for (int j=0; j<=i; ++j)
            {
              Dune::FieldVector<Scalar,n> y;
              y[j] = 1;
              auto Y = unpack(y);
              
              C[i][j] = static_cast<MaterialLaw const*>(this)->d2(X,Y);
            }
          }
          
          // ... and copy it to the upper half.
          for (int i=0; i<n-1; ++i)
            for (int j=i+1; j<n; ++j)
              C[i][j] = C[j][i];
            
          return C;
        }
      };
      
      // ---------------------------------------------------------------------------------------------------------
      
      template <class Material>
      class InvariantsMaterialLaw: public MaterialLawBase<Material::dim,typename Material::Scalar,InvariantsMaterialLaw<Material>>
      {
      public:
        using Scalar = typename Material::Scalar;
        static int const dim = Material::dim;
        using Tensor = Dune::FieldMatrix<Scalar,dim,dim>;
        
        template <class ... Args>
        InvariantsMaterialLaw(Args... args)
        : invariants(Tensor(0.0))
        , material(args...)
        {}
        
        void setLinearizationPoint(Tensor const& e)
        {
          invariants = ShiftedInvariants<dim,Scalar>(2.0*e);
          material.setLinearizationPoint(invariants.d0());
        }
        
        Scalar d0() const
        {
          return material.d0();
        }
        
        Scalar d1(Tensor const& e1) const
        {
          return 2*material.d1(invariants.d1(e1));
        }
        
        Scalar d2(Tensor const& e1, Tensor const& e2) const
        {
          return 4*(material.d2(invariants.d1(e1),invariants.d1(e2)) + material.d1(invariants.d2(e1,e2)));
        }
        
      private:
        ShiftedInvariants<dim,Scalar> invariants;
        Material                      material;
      };
      
      template <class Material, class DeviatoricPenalty>
      class CompressibleInvariantsMaterialLaw
      : public MaterialLawBase<Material::dim,typename Material::Scalar,CompressibleInvariantsMaterialLaw<Material,DeviatoricPenalty>>
      {
      public:
        using Scalar = typename Material::Scalar;
        static int const dim = Material::dim;
        using Tensor = Dune::FieldMatrix<Scalar,dim,dim>;
        
        template <class ... Args>
        CompressibleInvariantsMaterialLaw(Tensor const& E, Args... args): bE(E), material(bE.d0(),args...), deviator(bE.determinant().d0()) {}
        
        void setLinearizationPoint(Tensor const& e)
        {
          bE = IsochoricGreenLagrangeTensor<dim,Scalar>(e);
          deviator = DeviatoricPenalty(bE.determinant().d0());
          material.setLinearizationPoint(bE.d0());
        }
        
        Scalar d0() const
        {
          return material.d0() + deviator.d0();
        }
        
        Scalar d1(Tensor const& e1) const
        {
          return material.d1(bE.d1(e1)) + deviator.d1(bE.determinant().d1(e1));
        }
        
        Scalar d2(Tensor const& e1, Tensor const& e2) const
        {
          return material.d2(bE.d1(e1),bE.d1(e2)) + material.d1(bE.d2(e1,e2))
                 + deviator.d2(bE.determinant().d1(e1),bE.determinant().d1(e2)) + deviator.d1(bE.determinant().d2(e1,e2));
        }
        
      private:
        IsochoricGreenLagrangeTensor<dim,Scalar> bE;
        InvariantsMaterialLaw<Material>          material;
        DeviatoricPenalty                        deviator;
      };
      
      // ---------------------------------------------------------------------------------------------------------
      // ---------------------------------------------------------------------------------------------------------
      
      
      template <int dim, class Scalar=double>
      class StVenantKirchhoff: public MaterialLawBase<dim,Scalar,StVenantKirchhoff<dim,Scalar>>
      {
        using Det = Determinant<dim,Scalar>;
 
      public:
        
        using Tensor = Dune::FieldMatrix<Scalar,dim,dim>;
        
        StVenantKirchhoff() = default;
        
        StVenantKirchhoff(ElasticModulus const& moduli, Scalar alpha_=0)
        : lambda(moduli.lame()), mu(moduli.shear()), alpha(alpha_)
        , detC(unitMatrix<Scalar,dim>())
        {
        }
        
        StVenantKirchhoff(ElasticModulus const& moduli, Tensor const& e_)
        : lambda(moduli.lame()), mu(moduli.shear()), e(e_), tre(trace(e))
        , detC(unitMatrix<Scalar,dim>()+2*e)
        {
        }
        
        void setLinearizationPoint(Tensor const& e_)
        {
          e = e_;
          tre = trace(e);
          detC = Det(unitMatrix<Scalar,dim>()+2*e);
        }
        
        Scalar d0() const
        {
          Scalar w = lambda*tre*tre/2 + mu*contraction(e,e);
          if (alpha > 0)
          {
            auto detc = detC.d0();
            if (detc <= 0)
              w = std::numeric_limits<Scalar>::infinity();
            else
              w += alpha*(1/detc-1 + 2*tre);
          }
          return w;
        }
        
        Scalar d1(Tensor const& e1) const
        {
          Scalar dw = lambda*tre*trace(e1) + 2*mu*contraction(e,e1);
          if (alpha > 0)
            dw += alpha*(-std::pow(detC.d0(),-2)*detC.d1(e1)*2 + 2*trace(e1));
          return dw;
        }
        
        Scalar d2(Tensor const& e1, Tensor const& e2) const
        {
          Scalar ddw = lambda*trace(e1)*trace(e2) + 2*mu*contraction(e1,e2);
          if (alpha > 0)
            ddw += alpha*(  2*std::pow(detC.d0(),-3)*detC.d1(e1)*2*detC.d1(e2)*2
                          - std::pow(detC.d0(),-2)*detC.d2(e1,e2)*4 );
          return ddw;
        }
 
        Scalar d3(Tensor const& e1, Tensor const& e2, Tensor const& e3) const
        {
          if (alpha > 0)
            return alpha * ( -6*std::pow(detC.d0(),-4)*detC.d1(e1)*2*detC.d1(e2)*2*detC.d1(e3)*2
                             + 2*std::pow(detC.d0(),-3)*detC.d2(e1,e3)*4*detC.d1(e2)*2
                             + 2*std::pow(detC.d0(),-3)*detC.d1(e1)*2*detC.d2(e2,e3)*4
                             + 2*std::pow(detC.d0(),-3)*detC.d1(e3)*detC.d2(e1,e2)*4
                             - std::pow(detC.d0(),-2)*detC.d3(e1,e2,e3)*8 );
          else
            return 0;
        }
 
        Scalar d4(Tensor const& , Tensor const& , Tensor const&, Tensor const& ) const
        {
          // That's even easier for a quadratic energy...
          if (alpha > 0)
            abort(); // not yet implemented
          return 0;
        }
 
        private:
          Scalar lambda = 1;            // Lame parameter
          Scalar mu = 1;                // Lame parameter
          Scalar alpha = 0;             // compression penalty factor
          Tensor e = 0;                 // the strain tensor
          Scalar tre = 0;               // trace of strain tensor
          Determinant<dim,Scalar> detC; // determinant of Green-Lagrange tensor
      };
      
      // ---------------------------------------------------------------------------------------------------------
      
    template <int dim, class Scalar=double>
    class OrthotropicLinearMaterial : public MaterialLawBase<dim,Scalar,OrthotropicLinearMaterial<dim,Scalar>>
    {
      
      using Tensor = Dune::FieldMatrix<Scalar,dim,dim>;
    
    public:
      
      OrthotropicLinearMaterial(Dune::FieldMatrix<Scalar,dim,dim> const& orth_, Dune::FieldMatrix<Scalar,dim,dim> const& mat_)
      : orth(orth_), orthT(transpose(orth_)),  e(0)
      {
        Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> compliance(0);
        if (dim==3)
        {
          // see http://en.wikipedia.org/wiki/Orthotropic_material
          compliance[0][0] = 1 / mat_[0][0];             // 1/E_11
          compliance[0][1] = - mat_[1][0] / mat_[0][0];   // - nu_12 / E_11 = - nu_21 / E_22
          compliance[0][2] = - mat_[2][0] / mat_[0][0];   // - nu_13 / E_11 = - nu_31 / E_33
          compliance[1][0] = compliance[0][1];          // symmetry
          compliance[1][1] = 1 / mat_[1][1];
          compliance[1][2] = - mat_[2][1] / mat_[1][1];   // - nu_23 / E_22 = - nu_32 / E_33
          compliance[2][0] = compliance[0][2];          // symmetry
          compliance[2][1] = compliance[1][2];          // symmetry
          compliance[2][2] = 1 / mat_[2][2];
          
          compliance[3][3] = 1 / mat_[1][2];             // 1 / G_23
          compliance[4][4] = 1 / mat_[0][2];             // 1 / G_13
          compliance[5][5] = 1 / mat_[0][1];             // 1 / G_12
        } 
        else if (dim==2)
        {
          compliance[0][0] = 1 / mat_[0][0];             // 1/E_11;
          compliance[0][1] = - mat_[1][0] / mat_[0][0];   // - nu_12 / E_11 = - nu_21 / E_22
          compliance[1][0] = compliance[0][1];          // symmetry
          compliance[1][1] = 1 / mat_[1][1];             // 1/E_22
          
          compliance[2][2] = 1/mat_[0][1];               // 1 / G_12
        }
        
        // stiffness tensor is inverse of compliance tensor
        stiffness = compliance;
        stiffness.invert();
      }
      
      OrthotropicLinearMaterial(Dune::FieldMatrix<Scalar,dim,dim> const& orth_, 
                            Dune::FieldMatrix<Scalar,dim,dim> const& mat1, Dune::FieldVector<Scalar,dim*(dim-1)/2> const& mat2)
      : orth(orth_), orthT(transpose(orth_)), stiffness(0), e(0)
      {
        for (int i=0; i<dim; ++i)
          for (int j=0; j<dim; ++j)
            stiffness[i][j] = mat1[i][j];
        for (int i=0; i<dim*(dim-1)/2; ++i)
          stiffness[dim+i][dim+i] = mat2[i];
      }
      
      OrthotropicLinearMaterial(ElasticModulus const& p)
      : orth(unitMatrix<Scalar,dim>()), orthT(unitMatrix<Scalar,dim>()), stiffness(0), e(0) 
      {
        Scalar lambda = p.lame(), mu = p.shear();
        if (dim==2)
        {
          stiffness[0][0] = stiffness[1][1] = 2*mu + lambda;
          stiffness[1][0] = stiffness[0][1] = lambda;
          stiffness[2][2] = mu;
        }
        if (dim==3)
          abort();
      }
      
      OrthotropicLinearMaterial()
      : OrthotropicLinearMaterial(ElasticModulus()) {}
      
      Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> const& stiffnessMatrix() const { return stiffness; }
      
      void setLinearizationPoint(Tensor const& e_)
      {
        e = e_;
      }
      
      Scalar d0() const
      {
        auto eps = pack(orthT*e*orth);
        auto sigma = stiffness * eps;
        return 0.5 * (sigma*eps); 
      }
      
      Scalar d1(Tensor const& e1) const
      {
        auto eps  = pack(orthT*e*orth);
        auto eps1 = pack(orthT*e1*orth);
        auto sigma = stiffness * eps1;
        return sigma*eps; 
      }
      
      Scalar d2(Tensor const& e1, Tensor const &e2) const
      {
        auto eps1 = pack(orthT*e1*orth);
        auto eps2 = pack(orthT*e2*orth);
        auto sigma = stiffness * eps2;
        return sigma*eps1; 
      }
      
      Scalar d3(Tensor const& e1, Tensor const& e2, Tensor const &e3) const
      {
        return 0.0;
      }
      
      Scalar d4(Tensor const& e1, Tensor const& e2, Tensor const &e3, Tensor const &e4) const
      {
        return 0.0;
      }
      
    private:
      Dune::FieldMatrix<Scalar,dim,dim> orth; // the orthonormal matrix mapping material to global coordinates
      Dune::FieldMatrix<Scalar,dim,dim> orthT; // the transposed orthonormal matrix mapping material to global coordinates (for convenience)
      Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> stiffness;
      Dune::FieldMatrix<Scalar,dim,dim> e; // the linearization point 
    };
      // ---------------------------------------------------------------------------------------------------------
      
      template <int dim, class Scalar=double>
    class OrthotropicNonLinearMaterial : public MaterialLawBase<dim,Scalar,OrthotropicLinearMaterial<dim,Scalar>>
    {
      
      using Tensor = Dune::FieldMatrix<Scalar,dim,dim>;
    
    public:
      
      OrthotropicNonLinearMaterial(Dune::FieldMatrix<Scalar,dim,dim> const& orth_, Dune::FieldMatrix<Scalar,dim,dim> const& mat_)
      : orth(orth_), orthT(transpose(orth_)),  e(0), a(0), G(0)
      {
                 
 
 
                Scalar D = 0.75 * (2/(mat_[0][0]*mat_[1][1]) + 2/(mat_[1][1]*mat_[2][2]) + 2/(mat_[2][2]*mat_[0][0]) 
                               - 1/(mat_[0][0]*mat_[0][0]) - 1/(mat_[1][1]*mat_[1][1]) - 1/(mat_[2][2]*mat_[2][2]));
 
                                //Materialparameters aij
a[0][0] = 1/D * (2/(3*mat_[1][1]) + 2/(3*mat_[2][2]) - 1/(3*mat_[0][0]));    
a[0][1] = 1/D * (1/(6*mat_[0][0]) + 1/(6*mat_[1][1]) - 5/(6*mat_[2][2]));       
a[0][2] = 1/D * (1/(6*mat_[0][0]) + 1/(6*mat_[2][2]) - 5/(6*mat_[1][1]));
                
a[1][0] = 1/D * (1/(6*mat_[1][1]) + 1/(6*mat_[0][0]) - 5/(6*mat_[2][2]));       
a[1][1] = 1/D * (2/(3*mat_[2][2]) + 2/(3*mat_[0][0]) - 1/(3*mat_[1][1]));       
a[1][2] = 1/D * (1/(6*mat_[1][1]) + 1/(6*mat_[2][2]) - 5/(6*mat_[0][0]));
                
a[2][0] = 1/D * (1/(6*mat_[2][2]) + 1/(6*mat_[0][0]) - 5/(6*mat_[1][1]));   
a[2][1] = 1/D * (1/(6*mat_[2][2]) + 1/(6*mat_[1][1]) - 5/(6*mat_[0][0]));       
a[2][2] = 1/D * (2/(3*mat_[1][1]) + 2/(3*mat_[0][0]) - 1/(3*mat_[2][2]));
        
        //std::cout << a << "\n";
        
                //Materialparameters Gij
                                                                G[0][1] = mat_[0][1];   G[0][2] = mat_[0][2]; 
                G[1][0] = mat_[0][1];                                                   G[1][2] = mat_[1][2]; 
                G[2][0] = mat_[0][2];   G[2][1] = mat_[1][2]; 
                     
        // stiffness tensor see http://en.wikipedia.org/wiki/Orthotropic_material
        stiffness[0][0] = a[0][0] ; stiffness[0][1] = a[0][1] ; stiffness[0][2] = a[0][2] ;
        stiffness[1][0] = a[1][0] ; stiffness[1][1] = a[1][1] ; stiffness[1][2] = a[1][2] ;
        stiffness[2][0] = a[2][0] ; stiffness[2][1] = a[2][1] ; stiffness[2][2] = a[2][2] ;
        
        stiffness[3][3] = mat_[1][2];
        stiffness[4][4] = mat_[0][2];
        stiffness[5][5] = mat_[0][1];         
          
        //std::cout << stiffness << "\n";  
          
        // Structural tensor 
        x[0] =1;x[1] =0;x[2] =0;  
        y[0] =0;y[1] =1;y[2] =0;  
        z[0] =0;z[1] =0;z[2] =1;      
                        
        L[0] = outerProduct(orth*x,orth*x);
        L[1] = outerProduct(orth*y,orth*y);
        L[2] = outerProduct(orth*z,orth*z);
        
        /*
        std::cout << "L0"<< "\n" <<L[0] << "\n" << std::endl;
        std::cout << "L1"<< "\n" <<L[1] << "\n" << std::endl;
        std::cout << "L2"<< "\n" <<L[2] << "\n" << std::endl;
                */
        
       }
      
      OrthotropicNonLinearMaterial(Dune::FieldMatrix<Scalar,dim,dim> const& orth_, 
                            Dune::FieldMatrix<Scalar,dim,dim> const& mat1, Dune::FieldVector<Scalar,dim*(dim-1)/2> const& mat2)
      : orth(orth_), orthT(transpose(orth_)), stiffness(0), e(0)
      {
        for (int i=0; i<dim; ++i)
          for (int j=0; j<dim; ++j)
            stiffness[i][j] = mat1[i][j];
        for (int i=0; i<dim*(dim-1)/2; ++i)
          stiffness[dim+i][dim+i] = mat2[i];
      }
      
      OrthotropicNonLinearMaterial(ElasticModulus const& p)
      : orth(unitMatrix<Scalar,dim>()), orthT(unitMatrix<Scalar,dim>()), stiffness(0), e(0) 
      {
        Scalar lambda = p.lame(), mu = p.shear();
        if (dim==2)
        {
          stiffness[0][0] = stiffness[1][1] = 2*mu + lambda;
          stiffness[1][0] = stiffness[0][1] = lambda;
          stiffness[2][2] = mu;
        }
        if (dim==3)
          abort();
      }
      
      OrthotropicNonLinearMaterial()
      : OrthotropicNonLinearMaterial(ElasticModulus()) {}
      
      Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> const& stiffnessMatrix() const { return stiffness; }
      
      void setLinearizationPoint(Tensor const& e_)
      {
        e = e_;
      }
      
      Scalar d0() const
      {
                        
            Scalar shear = 0 , eps = 0;
            Dune::FieldMatrix<Scalar, 3 ,3> A(0),B(0);
        
                for(unsigned int i = 0 ; i < dim ; ++i)
                {
                        for(unsigned int j = 0 ; j < dim ; ++j)
                        {
                                //eps += a[i][j]*contraction(transpose(e),L[i])*contraction(transpose(e),L[j]);
                                eps += a[i][j]*trace(e*(L[i]))*trace(e*(L[j]));
                        }
                }               
                                
                for(unsigned int i = 0 ; i < dim ; ++i)
                {
                        auto A = e*L[i];
                        for(unsigned int j = 0 ; j < dim ; ++j)
                        {
                                if(i==j)
                                {continue;}             
                                
                                auto B = e.rightmultiplyany(L[j]);
                                
                                shear += G[i][j]*trace(A.rightmultiplyany(B));
                        }
                }                       
                
 
        return 0.5 * eps + shear;
      }
      
      Scalar d1(Tensor const& e1) const
      {
            Scalar shear = 0 , eps = 0;
            Dune::FieldMatrix<Scalar, 3 ,3> A(0), B(0), C(0), D(0);
                  
                for(unsigned int i = 0 ; i < dim ; ++i)
                {
                        for(unsigned int j = 0 ; j < dim ; ++j)
                        {
                                eps += a[i][j]*(trace(e*(L[i]))*trace(e1*(L[j])) 
                                                           +trace(e1*(L[i]))*trace(e*(L[j])));
                        }
                }       
                
                for(unsigned int i = 0 ; i < dim ; ++i)
                {
                        auto A = e.rightmultiplyany(L[i]);
                        auto C = e1.rightmultiplyany(L[i]);
                        
                        for(unsigned int j = 0 ; j < dim ; ++j)
                        {
                                if(i==j)
                                {continue;}
                                                                
                                auto B = e1.rightmultiplyany(L[j]);
                                auto D = e.rightmultiplyany(L[j]);
                                
                                shear += G[i][j]*(trace(A.rightmultiplyany(B))+trace(C.rightmultiplyany(D)));
                        }
                }
 
                
                
 
        return 0.5 * eps + shear;
      }
      
      Scalar d2(Tensor const& e1, Tensor const &e2) const
      {
            
            Scalar shear = 0 , eps = 0;
            Dune::FieldMatrix<Scalar, 3 ,3> A(0),B(0),C(0),D(0);
                  
                for(unsigned int i = 0 ; i < dim ; ++i)
                {
                        for(unsigned int j = 0 ; j < dim ; ++j)
                        {
                                eps += a[i][j]*(trace(e2*L[i])*trace(e1*L[j]) 
                                                           +trace(e1*L[i])*trace(e2 *L[j]));
                        }
                }       
                
                for(unsigned int i = 0 ; i < dim ; ++i)
                {
                        auto A = e2*L[i];
                        auto C = e1*L[i];
                        
                        for(unsigned int j = 0 ; j < dim ; ++j)
                        {
                                if(i==j)
                                {continue;}                             
                                
                                auto B = e1*L[j];
                                auto D = e2*L[j];
                                
                                shear += G[i][j]*(trace(A.rightmultiplyany(B))+trace(C.rightmultiplyany(D)));
                        }
                }       
 
        return 0.5 * eps + shear;
      }
      
      Scalar d3(Tensor const& e1, Tensor const& e2, Tensor const &e3) const
      {
        return 0.0;
      }
      
      Scalar d4(Tensor const& e1, Tensor const& e2, Tensor const &e3, Tensor const &e4) const
      {
        return 0.0;
      }
      
    private:
      Dune::FieldMatrix<Scalar,dim,dim> orth; // the orthonormal matrix mapping material to global coordinates
      Dune::FieldMatrix<Scalar,dim,dim> orthT; // the transposed orthonormal matrix mapping material to global coordinates (for convenience)
      Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> stiffness;
      Dune::FieldMatrix<Scalar,dim,dim> e; // the linearization point 
      Dune::FieldMatrix<Scalar,dim,dim> a; // Materialparameters
          Dune::FieldMatrix<Scalar,dim,dim> G; // Materialparameters
          Dune::FieldVector< Dune::FieldMatrix<Scalar,dim,dim> , 3 > L;//Structural tensor
          Dune::FieldVector< Scalar , dim > x,y,z;
 
    };
      // ---------------------------------------------------------------------------------------------------------
      
      template <int dimension>
      class MooneyRivlin
      {
      public:
        using Scalar = double;
        
        static int const dim = dimension;
        
        using Invariants = Dune::FieldVector<Scalar,dim>;
        
        MooneyRivlin(double c1_, double c2_, double c3_=0.0)
        : c1(c1_), c2(dim==3?c2_:0), c3(c3_)
        {}
        
        void setLinearizationPoint(Invariants const& i)
        {
          i1 = i[0];
          i2 = i[1];
          id = i[dim-1];
        }
        
        Scalar d0() const
        {
          return c1*i1 + c2*i2 + c3*id*id/(1+id);
        }
        
        Scalar d1(Invariants const& di1) const
        {
          return c1*di1[0] + c2*di1[1] + c3*id*(2+id)/square(1+id)*di1[dim-1];
        }
        
        Scalar d2(Invariants const& di1, Invariants const& di2) const
        {
          return 2 / power(1+id,3) *di1[dim-1]*di2[dim-1];
        }
        
      private:
        double i1, i2, id, c1, c2, c3;
      };
      
      // ------------------------------------------------------------------------------------------
 
      template <int dimension>
      class NeoHookean: public MooneyRivlin<dimension>
      {
      public:
        using typename MooneyRivlin<dimension>::Invariants;
 
        NeoHookean(ElasticModulus const& moduli): MooneyRivlin<dimension>(moduli.shear()/2,0)
        {}
      };
 
      // ------------------------------------------------------------------------------------------
      
      template <int dimension>
      class BlatzKo
      {
      public:
        using Scalar = double;
        
        static int const dim = dimension;
        
        using Invariants = Dune::FieldVector<Scalar,dim>;
        
        BlatzKo(ElasticModulus const& moduli)
        : mu(moduli.shear()), a(2*moduli.poisson()/(1-2*moduli.poisson())), pstable(1,0)
        {
          assert(a>0);
        }
        
        void setLinearizationPoint(Invariants const& i)
        {
          i1 = i[0];
          i3 = i[dimension-1];
          pstable = Pstable(-a/2,i3);
        }
        
        Scalar d0() const
        {
          return mu*(i1 + 2*pstable.d0()/a)/2;
        }
        
        Scalar d1(Invariants const& di1) const
        {
          return mu*(di1[0] + 2*pstable.d1(di1[dimension-1])/a)/2;
        }
        
        Scalar d2(Invariants const& di1, Invariants const& di2) const
        {
          return mu*pstable.d2(di1[dimension-1],di2[dimension-1])/a;
        }
        
      private:
        double i1, i3, mu, a;
        Pstable pstable;
      };
      
      
      // ---------------------------------------------------------------------------------------------------------
      
      template <class HyperelasticEnergy>
      class ViscoPlasticEnergy: public HyperelasticEnergy
      {
      public:
        using Scalar = typename HyperelasticEnergy::Scalar;
        static int const dim = HyperelasticEnergy::dim;
        
        using Tensor = Dune::FieldMatrix<Scalar,dim,dim>;
        
        ViscoPlasticEnergy() = default;
        
        template<typename... Args>
        ViscoPlasticEnergy(Tensor const& evp_, Args&&... args): HyperelasticEnergy(std::forward<Args>(args)...), evp(evp_)
        {}
        
        void setLinearizationPoint(Tensor const& e)
        {
          HyperelasticEnergy::setLinearizationPoint(e-evp);
        }
                
      private:
        Tensor             evp;
      };
      
      // ---------------------------------------------------------------------------------------------------------
      
      template <int dim, class Scalar>
      Scalar vonMisesStress(Dune::FieldMatrix<Scalar,dim,dim> const& stress);
 
      // ---------------------------------------------------------------------------------------------------------
 
      template <int dim, class Scalar=double>
      Dune::FieldVector<Scalar,dim*(dim+1)/2> 
      duvautLionsJ2Flow(Dune::FieldMatrix<Scalar,dim,dim> const& cauchyStress, Dune::FieldMatrix<Scalar,dim*(dim+1)/2,dim*(dim+1)/2> const& C, double tau, double sigmaY);
      
      // ---------------------------------------------------------------------------------------------------------
      
    }
  }
}
 
#endif