standard_tensors.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2020 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_elasticity_standard_tensors_h
#define dealii_elasticity_standard_tensors_h


#include <deal.II/base/config.h>

#include <deal.II/base/numbers.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 DEAL_II_CONSTEXPR const SymmetricTensor<2, dim> I
#ifdef DEAL_II_HAVE_CXX14_CONSTEXPR
        = unit_symmetric_tensor<dim>()
#endif
        ;

      static DEAL_II_CONSTEXPR const SymmetricTensor<4, dim> S
#ifdef DEAL_II_HAVE_CXX14_CONSTEXPR
        = identity_tensor<dim>()
#endif
        ;

      static DEAL_II_CONSTEXPR const SymmetricTensor<4, dim> IxI
#ifdef DEAL_II_HAVE_CXX14_CONSTEXPR
        = outer_product(unit_symmetric_tensor<dim>(),
                        unit_symmetric_tensor<dim>())
#endif
        ;



      static DEAL_II_CONSTEXPR const SymmetricTensor<4, dim> dev_P
#ifdef DEAL_II_HAVE_CXX14_CONSTEXPR
        = deviator_tensor<dim>()
#endif
        ;

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

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


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



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

    };

  } // namespace Elasticity
} // namespace Physics



#ifndef DOXYGEN

// --------------------- inline functions and constants -------------------


template <int dim>
template <typename Number>
DEAL_II_CONSTEXPR inline SymmetricTensor<4, dim, Number>
Physics::Elasticity::StandardTensors<dim>::Dev_P(
  const Tensor<2, dim, Number> &F)
{
  const Number det_F = determinant(F);
  Assert(numbers::value_is_greater_than(det_F, 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>
DEAL_II_CONSTEXPR 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(numbers::value_is_greater_than(det_F, 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>
DEAL_II_CONSTEXPR 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>
DEAL_II_CONSTEXPR 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