elasto.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) 2012-2023 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef ELASTO_HH
#define ELASTO_HH
 
#include <cassert>
 
#include "dune/common/config.h"
#include "dune/common/fmatrix.hh"
#include "dune/grid/config.h"
 
#include "fem/diffops/materialProperties.hh"
#include "fem/fixdune.hh"
#include "fem/integration.hh"
 
#include "linalg/determinant.hh"
 
#include "utilities/power.hh"
 
namespace Kaskade 
{
  namespace Elastomechanics
  {
    // ---------------------------------------------------------------------------------------------------------
    
    class ElasticModulus
    {
    public:
      ElasticModulus(): lambda(1), mu(1) {}
      
      ElasticModulus(double lambda_, double mu_): lambda(lambda_), mu(mu_) {}
 
      ElasticModulus(std::string const& name)
      : ElasticModulus(ElasticModulus::material(name))
      {}
      
      ElasticModulus& setYoungPoisson(double E, double nu);
      
      double shear() const { return mu; }
      
      double bulk() const { return lambda + 2*mu/3; }
      
      double young() const { return mu*(3*lambda+2*mu)/(lambda+mu); }
      
      double lame() const { return lambda; }
      
      double poisson() const { return lambda/2/(lambda+mu); }
      
      static ElasticModulus const& material(std::string const& name);
      
      static std::map<std::string,ElasticModulus> const& materials();
      
    private:
      double lambda, mu;
    };
    
    double yieldStrength(std::string const& name);
    
    
    // ---------------------------------------------------------------------------------------------------------
    
    namespace ElastoDetails
    {
      // Computes the dimension of a matrix from the length of its Voigt representation vector.
      constexpr int dim(int n) { return (sqrti(1+8*n)-1)/2; }
    }
    // ---------------------------------------------------------------------------------------------------------
    
    template <class Scalar, int n>
    Dune::FieldVector<Scalar,n*(n+1)/2> pack(Dune::FieldMatrix<Scalar,n,n> const& e, Scalar s=2.0);
    
    template <class Scalar, int n>
    Dune::FieldMatrix<Scalar,ElastoDetails::dim(n),ElastoDetails::dim(n)> unpack(Dune::FieldVector<Scalar,n> const& c, Scalar s=2.0);
    
    template <int d, class Scalar>
    Dune::FieldMatrix<Scalar,d,d> deviatoricPart(Dune::FieldMatrix<Scalar,d,d> const& s)
    {
      return s - (trace(s)/d) * unitMatrix<Scalar,d>();
    }
    
    // --------------------------------------------------------------------------------------------
    
    template <int dim, class Scalar>
    class DetIpm1
    {
    public:
      using Matrix = Dune::FieldMatrix<Scalar,dim,dim>;
    
      DetIpm1(Matrix const& A_): A(A_) {}
      
      Scalar d0() const;
      Scalar d1(Matrix const& dA) const;
      Scalar d2(Matrix const& dA1, Matrix const& dA2) const;
    private:
      Matrix A;
    };
    
    class Pstable
    {
    public:
      Pstable(double p_, double x_);
      double d0() const { return f; }
      double d1(double dx) const { return df*dx; }
      double d2(double dx1, double dx2) const { return ddf*dx1*dx2; }
    private:
      double f, df, ddf;
    };
    
    // --------------------------------------------------------------------------------------------
    
    template <int dim, class Scalar=double>
    class ShiftedInvariants
    {
    public:
      using Tensor = Dune::FieldMatrix<Scalar,dim,dim>;
      using Invariants = Dune::FieldVector<Scalar,dim>;
      
      ShiftedInvariants(Tensor const& A_): A(A_), det(A) {}
      Invariants d0() const;
      Invariants d1(Tensor const& dA) const;
      Invariants d2(Tensor const& dA1, Tensor const& dA2) const;
      
    public:
      Tensor A;
      DetIpm1<dim,Scalar> det;
    };
  
    // --------------------------------------------------------------------------------------------
 
    template <class Scalar, int dim, bool byValue=true>
    class LinearizedGreenLagrangeTensor
    {
    public:
      typedef Dune::FieldMatrix<Scalar,dim,dim> Tensor;
      
      explicit LinearizedGreenLagrangeTensor(Tensor const& du_): du(du_) {}
      
      LinearizedGreenLagrangeTensor(): du(0) {}
      
      
      void setDisplacementDerivative(Tensor const& du_)
      {
        du = du_;
      }
      
      Tensor d0() const { return static_cast<Scalar>(0.5)*(du+transpose(du)); }
      
      Tensor d1(Tensor const& dv) const { return static_cast<Scalar>(0.5)*(dv+transpose(dv)); }
      
      Tensor d2(Tensor const&, Tensor const&) const { return Tensor(0.0); }
      
    private:
      std::conditional_t<byValue,Tensor,Tensor const&> du;
    };
    
    template <class Scalar, int dim, bool byValue=true>
    class GreenLagrangeTensor
    {
    public:
      typedef Dune::FieldMatrix<Scalar,dim,dim> Tensor;
      
      GreenLagrangeTensor(Tensor const& du_): du(du_) {}
      
      GreenLagrangeTensor(): du(0) {}
      
      [[gnu::deprecated("create a new LinearizedGreenLagrangeTensor if needed")]]
      void setDisplacementDerivative(Tensor const& du_)
      {
        du = du_;
      }
      
      Tensor d0() const { return static_cast<Scalar>(0.5)*(du+transpose(du)+transpose(du)*du); }
      
      Tensor d1(Tensor const& dv) const { return static_cast<Scalar>(0.5)*(dv+transpose(dv)+transpose(dv)*du + transpose(du)*dv); }
      
      Tensor d2(Tensor const& dv, Tensor const& dw) const { return static_cast<Scalar>(0.5)*(transpose(dv)*dw + transpose(dw)*dv); }
      
    private:
      std::conditional_t<byValue,Tensor,Tensor const&> du;
    };
    
    
    template <class Scalar, int dim, bool byValue=true>
    class SurfaceGreenLagrangeTensor
    {
    public:
      typedef Dune::FieldMatrix<Scalar,dim+1,dim> InputTensor;
      typedef Dune::FieldMatrix<Scalar,dim,dim> OutputTensor;
      
      SurfaceGreenLagrangeTensor(InputTensor const& du_): du(du_), dI(0) 
      {
        for(int i=0; i < dim; ++i)
          dI[i][i] = 1;
//         std::cout << "dI=" << dI << "\n";
//         std::cout << "transpose(dI)=" << transpose(dI) << "\n";
      }
      
      SurfaceGreenLagrangeTensor(): du(0), dI(0) 
      {
        for(int i=0; i < dim; ++i)
          dI[i][i] = 1;
      }
      
      
      OutputTensor d0() const { return static_cast<Scalar>(0.5)*(transpose(dI)*du+transpose(du)*dI+transpose(du)*du); }
      
      OutputTensor d1(InputTensor const& dv) const { return static_cast<Scalar>(0.5)*(transpose(dI)*dv+transpose(dv)*dI+transpose(dv)*du + transpose(du)*dv); }
      
      OutputTensor d2(InputTensor const& dv, InputTensor const& dw) const { return static_cast<Scalar>(0.5)*(transpose(dv)*dw + transpose(dw)*dv); }
      
    private:
      std::conditional_t<byValue,InputTensor,InputTensor const&> du;
      InputTensor dI;
    };
    
    template <class Scalar, int dim, bool byValue=true>
    class ExtendedGreenLagrangeTensor
    {
    public:
      typedef Dune::FieldMatrix<Scalar,dim+1,dim> InputTensor;
      typedef Dune::FieldMatrix<Scalar,dim+1,dim+1> OutputTensor;
      
      ExtendedGreenLagrangeTensor(InputTensor const& du_): du(du_), duext(extend(du_))
      {
      }
      
      ExtendedGreenLagrangeTensor(): du(0), duext(0) 
      {        
      }
      
      
      OutputTensor d0() const 
      {
        return static_cast<Scalar>(0.5)*(duext+transpose(duext)+transpose(duext)*duext);         
      }
      
      OutputTensor d1(InputTensor const& dv) const 
      {
        OutputTensor dvext = extend(dv);
        return static_cast<Scalar>(0.5)*(dvext+transpose(dvext)+transpose(dvext)*duext + transpose(duext)*dvext); 
      }
      
      OutputTensor d2(InputTensor const& dv, InputTensor const& dw) const 
      {
        OutputTensor dvext = extend(dv);
        OutputTensor dwext = extend(dw); 
        return static_cast<Scalar>(0.5)*(transpose(dvext)*dwext + transpose(dwext)*dvext);         
      }
      
    private:
      
      OutputTensor extend(InputTensor const& d) const
      {
        OutputTensor dext(0);
        for(int i=0; i < dim+1; i++)
          for(int j=0; j < dim; j++)
            dext[i][j] = d[i][j];
        
        dext[0][dim] = -dext[dim][0];
        dext[1][dim] = -dext[dim][1];
        return dext;
      }
      
      std::conditional_t<byValue,InputTensor,InputTensor const&> du;
      OutputTensor duext;
    };
    
    
    template <int dim, class Scalar=double>
    class IsochoricGreenLagrangeTensor
    {
    public:
      typedef Dune::FieldMatrix<Scalar,dim,dim> Tensor;
      
      IsochoricGreenLagrangeTensor(Tensor const& E_): E(E_), det(2*E), z(-1.0/dim,det.d0())
      {
      }
      
      IsochoricGreenLagrangeTensor(): E(0), z(0.0,0.0) {}
      
      Tensor d0() const { 
        return E + z.d0()*(E + static_cast<Scalar>(0.5)*unitMatrix<Scalar,dim>());
      }
      
      Tensor d1(Tensor const& dE) const 
      { 
        return dE + z.d1(det.d1(2*dE))*(E + static_cast<Scalar>(0.5)*unitMatrix<Scalar,dim>()) + z.d0()*dE;
      }
      
      Tensor d2(Tensor const& dE1, Tensor const& dE2) const 
      { 
        return z.d2(det.d1(2*dE1,2*dE2))*(E + static_cast<Scalar>(0.5)*unitMatrix<Scalar,dim>()) + z.d1(det.d1(2*dE1))*dE2 + z.d1(det.d1(2*dE2))*dE1;
      }
      
      DetIpm1<dim,Scalar> const& determinant() const
      {
        return det;
      }
      
    private:
      Tensor E;
      DetIpm1<dim,Scalar> det;
      Pstable z;
    };
    
    // --------------------------------------------------------------------------------------------
    
    
    template <class Displacement, class StrainTensor>
    class StrainView
    {
    public:
      using Space = typename Displacement::Space;
      using ValueType = typename StrainTensor::Tensor;
      
      
      StrainView(Displacement const& u_)
      : u(u_) {}
      
      Space const& space() const
      {
        return u.space();
      }
      
      ValueType value(typename Space::Evaluator const& evaluator) const
      {
        return StrainTensor(u.derivative(evaluator)).d0();
      }
      
    private:
      Displacement const& u;
    };
    
    template <class Displacement>
    auto makeLinearizedGreenLagrangeStrainView(Displacement const& u)
    {
      return StrainView<Displacement,LinearizedGreenLagrangeTensor<typename Displacement::Scalar,Displacement::Space::GridView::dimension>>(u);
    }
    
    template <class Displacement>
    auto makeGreenLagrangeStrainView(Displacement const& u)
    {
      return StrainView<Displacement,GreenLagrangeTensor<typename Displacement::Scalar,Displacement::Space::GridView::dimension>>(u);
    }
    
    // ---------------------------------------------------------------------------------------------------------
    
    template <class Real, int dim>
    std::array<Dune::FieldVector<Real,dim>,dim> principalDirections(Dune::FieldMatrix<Real,dim,dim> const& A);
    
    // ---------------------------------------------------------------------------------------------------------
    
    
    namespace Elasto_Details {
      
      template <class Real, int n> 
      Real maxOrientationPreservingStepsize(Dune::FieldMatrix<Real,n,n> const& A, Dune::FieldMatrix<Real,n,n> const& dA, Real eps);
      
      template <class Function>
      struct Limiter {
        static int const dim = Function::components;
        
        Limiter(Function const& y_, Function const& dy_, double epsilon_)
        : y(y_), dy(dy_), epsilon(epsilon_), step(1), I(unitMatrix<typename Function::Scalar,dim>()) { mindet = std::numeric_limits<double>::max(); }
        
        template <class Cell, class QuadraturePoint, class Evaluator>
        void operator()(Cell const&, QuadraturePoint const&, Evaluator const& eval) {
          auto yx = y.derivative(eval);
          auto dyx = dy.derivative(eval);
          auto tmp = determinant(I+yx);
          
          if (tmp <= 0)
            std::cerr << "#### AIEEE: given determinant nonpositive: " << tmp << " !\n";
          
          mindet = std::min(mindet,tmp);
          while (step>0 && determinant(I+yx+step*dyx)<epsilon*tmp)
            step = 0.7*step - 1e-4;
          
          for (int i=1; i<10; ++i)
            if (determinant(I+yx+i*step/10*dyx)<epsilon*tmp)
              std::cerr << "##### Aieee: at step " << step << " * " << i << "/10 we have " << determinant(I+yx+i*step/10*dyx) << " < " << epsilon*tmp << "\n";
        }
        
        template <class Evaluator>
        int order(Evaluator const& eval) const { return y.order(eval); }
        
        Function const& y, dy;
        double epsilon, step, mindet;
        Dune::FieldMatrix<typename Function::Scalar,dim,dim> I;
      };
      
    }
    
    template <class Function>
    double orientationPreservingStepsize(Function const& y, Function const& dy, double epsilon)
    {
      
      Elasto_Details::Limiter<Function> g(y,dy,epsilon);
      forEachQP(g,y.space());
                     std::cout << "minimal determinant: " << g.mindet << "\n";
      return std::max(g.step,0.0);
    }
    
  // ---------------------------------------------------------------------------------------------------------
  }
  
  // TODO: remove this backwards compatibility hack as soon as all users switched over
  using namespace Elastomechanics;
  
}
 
 
  
 
#endif