standard_tensors.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__elasticity__standard_tensors_h
#define dealii__elasticity__standard_tensors_h


#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>

DEAL_II_NAMESPACE_OPEN

namespace Physics
{

  namespace Elasticity
  {

    template <int dim>
    class StandardTensors
    {
    public:


      static const SymmetricTensor<2, dim> I;

      static const SymmetricTensor<4, dim> S;

      static const SymmetricTensor<4, dim> IxI;



      static const SymmetricTensor<4, dim> dev_P;

      template <typename Number>
      static SymmetricTensor<4, dim, Number>
      Dev_P (const Tensor<2, dim, Number> &F);

      template <typename Number>
      static SymmetricTensor<4, dim, Number>
      Dev_P_T (const Tensor<2, dim, Number> &F);


      template <typename Number>
      static SymmetricTensor<2, dim, Number>
      ddet_F_dC (const Tensor<2, dim, Number> &F);



      template <typename Number>
      static SymmetricTensor<4, dim, Number>
      dC_inv_dC (const Tensor<2, dim, Number> &F);

    };

  }
}



#ifndef DOXYGEN

// ------------------------- inline functions ------------------------


template <int dim>
template <typename Number>
inline
SymmetricTensor<4, dim, Number>
Physics::Elasticity::StandardTensors<dim>::Dev_P (const Tensor<2, dim, Number> &F)
{
  const Number det_F = determinant(F);
  Assert(det_F > Number(0.0),
         ExcMessage("Deformation gradient has a negative determinant."));
  const Tensor<2,dim,Number> C_ns = transpose(F)*F;
  const SymmetricTensor<2,dim,Number> C = symmetrize(C_ns);
  const SymmetricTensor<2,dim,Number> C_inv = symmetrize(invert(C_ns));

  // See Wriggers p46 equ 3.125 (but transpose indices)
  SymmetricTensor<4,dim,Number> Dev_P = outer_product(C,C_inv);  // Dev_P = C_x_C_inv
  Dev_P /= -dim;                                                 // Dev_P = -[1/dim]C_x_C_inv
  Dev_P += SymmetricTensor<4,dim,Number>(S);                     // Dev_P = S - [1/dim]C_x_C_inv
  Dev_P *= std::pow(det_F, -2.0/dim);                            // Dev_P = J^{-2/dim} [S - [1/dim]C_x_C_inv]

  return Dev_P;
}



template <int dim>
template <typename Number>
inline
SymmetricTensor<4, dim, Number>
Physics::Elasticity::StandardTensors<dim>::Dev_P_T (const Tensor<2, dim, Number> &F)
{
  const Number det_F = determinant(F);
  Assert(det_F > Number(0.0),
         ExcMessage("Deformation gradient has a negative determinant."));
  const Tensor<2,dim,Number> C_ns = transpose(F)*F;
  const SymmetricTensor<2,dim,Number> C = symmetrize(C_ns);
  const SymmetricTensor<2,dim,Number> C_inv = symmetrize(invert(C_ns));

  // See Wriggers p46 equ 3.125 (not transposed)
  SymmetricTensor<4,dim,Number> Dev_P_T = outer_product(C_inv,C);  // Dev_P = C_inv_x_C
  Dev_P_T /= -dim;                                                 // Dev_P = -[1/dim]C_inv_x_C
  Dev_P_T += SymmetricTensor<4,dim,Number>(S);                     // Dev_P = S - [1/dim]C_inv_x_C
  Dev_P_T *= std::pow(det_F, -2.0/dim);                            // Dev_P = J^{-2/dim} [S - [1/dim]C_inv_x_C]

  return Dev_P_T;
}



template <int dim>
template <typename Number>
inline
SymmetricTensor<2, dim, Number>
Physics::Elasticity::StandardTensors<dim>::ddet_F_dC (const Tensor<2, dim, Number> &F)
{
  return Number(0.5)*determinant(F)*symmetrize(invert(transpose(F)*F));
}



template <int dim>
template <typename Number>
inline
SymmetricTensor<4, dim, Number>
Physics::Elasticity::StandardTensors<dim>::dC_inv_dC (const Tensor<2, dim, Number> &F)
{
  const SymmetricTensor<2,dim,Number> C_inv = symmetrize(invert(transpose(F)*F));

  SymmetricTensor<4,dim,Number> dC_inv_dC;
  for (unsigned int A=0; A<dim; ++A)
    for (unsigned int B=A; B<dim; ++B)
      for (unsigned int C=0; C<dim; ++C)
        for (unsigned int D=C; D<dim; ++D)
          dC_inv_dC[A][B][C][D] -= 0.5*(C_inv[A][C] * C_inv[B][D] + C_inv[A][D] * C_inv[B][C] );

  return dC_inv_dC;
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif