transformations.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2017 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_transformations_h
#define dealii_transformations_h

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

DEAL_II_NAMESPACE_OPEN


namespace Physics
{
  namespace Transformations
  {
    namespace Rotations
    {

      template <typename Number>
      Tensor<2, 2, Number>
      rotation_matrix_2d(const Number &angle);


      template <typename Number>
      Tensor<2, 3, Number>
      rotation_matrix_3d(const Point<3, Number> &axis, const Number &angle);


    } // namespace Rotations

    namespace Contravariant
    {

      template <int dim, typename Number>
      Tensor<1, dim, Number>
      push_forward(const Tensor<1, dim, Number> &V,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      Tensor<2, dim, Number>
      push_forward(const Tensor<2, dim, Number> &T,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<2, dim, Number>
      push_forward(const SymmetricTensor<2, dim, Number> &T,
                   const Tensor<2, dim, Number> &         F);

      template <int dim, typename Number>
      Tensor<4, dim, Number>
      push_forward(const Tensor<4, dim, Number> &H,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<4, dim, Number>
      push_forward(const SymmetricTensor<4, dim, Number> &H,
                   const Tensor<2, dim, Number> &         F);



      template <int dim, typename Number>
      Tensor<1, dim, Number>
      pull_back(const Tensor<1, dim, Number> &v,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      Tensor<2, dim, Number>
      pull_back(const Tensor<2, dim, Number> &t,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<2, dim, Number>
      pull_back(const SymmetricTensor<2, dim, Number> &t,
                const Tensor<2, dim, Number> &         F);

      template <int dim, typename Number>
      Tensor<4, dim, Number>
      pull_back(const Tensor<4, dim, Number> &h,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<4, dim, Number>
      pull_back(const SymmetricTensor<4, dim, Number> &h,
                const Tensor<2, dim, Number> &         F);

    } // namespace Contravariant

    namespace Covariant
    {

      template <int dim, typename Number>
      Tensor<1, dim, Number>
      push_forward(const Tensor<1, dim, Number> &V,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      Tensor<2, dim, Number>
      push_forward(const Tensor<2, dim, Number> &T,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<2, dim, Number>
      push_forward(const SymmetricTensor<2, dim, Number> &T,
                   const Tensor<2, dim, Number> &         F);

      template <int dim, typename Number>
      Tensor<4, dim, Number>
      push_forward(const Tensor<4, dim, Number> &H,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<4, dim, Number>
      push_forward(const SymmetricTensor<4, dim, Number> &H,
                   const Tensor<2, dim, Number> &         F);



      template <int dim, typename Number>
      Tensor<1, dim, Number>
      pull_back(const Tensor<1, dim, Number> &v,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      Tensor<2, dim, Number>
      pull_back(const Tensor<2, dim, Number> &t,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<2, dim, Number>
      pull_back(const SymmetricTensor<2, dim, Number> &t,
                const Tensor<2, dim, Number> &         F);

      template <int dim, typename Number>
      Tensor<4, dim, Number>
      pull_back(const Tensor<4, dim, Number> &h,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<4, dim, Number>
      pull_back(const SymmetricTensor<4, dim, Number> &h,
                const Tensor<2, dim, Number> &         F);

    } // namespace Covariant

    namespace Piola
    {

      template <int dim, typename Number>
      Tensor<1, dim, Number>
      push_forward(const Tensor<1, dim, Number> &V,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      Tensor<2, dim, Number>
      push_forward(const Tensor<2, dim, Number> &T,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<2, dim, Number>
      push_forward(const SymmetricTensor<2, dim, Number> &T,
                   const Tensor<2, dim, Number> &         F);

      template <int dim, typename Number>
      Tensor<4, dim, Number>
      push_forward(const Tensor<4, dim, Number> &H,
                   const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<4, dim, Number>
      push_forward(const SymmetricTensor<4, dim, Number> &H,
                   const Tensor<2, dim, Number> &         F);



      template <int dim, typename Number>
      Tensor<1, dim, Number>
      pull_back(const Tensor<1, dim, Number> &v,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      Tensor<2, dim, Number>
      pull_back(const Tensor<2, dim, Number> &t,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<2, dim, Number>
      pull_back(const SymmetricTensor<2, dim, Number> &t,
                const Tensor<2, dim, Number> &         F);

      template <int dim, typename Number>
      Tensor<4, dim, Number>
      pull_back(const Tensor<4, dim, Number> &h,
                const Tensor<2, dim, Number> &F);

      template <int dim, typename Number>
      SymmetricTensor<4, dim, Number>
      pull_back(const SymmetricTensor<4, dim, Number> &h,
                const Tensor<2, dim, Number> &         F);

    } // namespace Piola


    template <int dim, typename Number>
    Tensor<1, dim, Number>
    nansons_formula(const Tensor<1, dim, Number> &N,
                    const Tensor<2, dim, Number> &F);

  } // namespace Transformations
} // namespace Physics



#ifndef DOXYGEN

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

namespace internal
{
  namespace Physics
  {
    template <int dim, typename Number>
    inline Tensor<1, dim, Number>
    transformation_contraction(const Tensor<1, dim, Number> &V,
                               const Tensor<2, dim, Number> &F)
    {
      return contract<1, 0>(F, V);
    }



    template <int dim, typename Number>
    inline Tensor<2, dim, Number>
    transformation_contraction(const Tensor<2, dim, Number> &T,
                               const Tensor<2, dim, Number> &F)
    {
      return contract<1, 0>(F, contract<1, 1>(T, F));
    }



    template <int dim, typename Number>
    inline ::SymmetricTensor<2, dim, Number>
    transformation_contraction(const ::SymmetricTensor<2, dim, Number> &T,
                               const Tensor<2, dim, Number> &                 F)
    {
      Tensor<2, dim, Number> tmp_1;
      for (unsigned int i = 0; i < dim; ++i)
        for (unsigned int J = 0; J < dim; ++J)
          for (unsigned int I = 0; I < dim; ++I)
            tmp_1[i][J] += F[i][I] * T[I][J];

      ::SymmetricTensor<2, dim, Number> out;
      for (unsigned int i = 0; i < dim; ++i)
        for (unsigned int j = i; j < dim; ++j)
          for (unsigned int J = 0; J < dim; ++J)
            out[i][j] += F[j][J] * tmp_1[i][J];

      return out;
    }



    template <int dim, typename Number>
    inline Tensor<4, dim, Number>
    transformation_contraction(const Tensor<4, dim, Number> &H,
                               const Tensor<2, dim, Number> &F)
    {
      // This contraction order and indexing might look a bit dubious, so a
      // quick explanation as to what's going on is probably in order:
      //
      // When the contract() function operates on the inner indices, the
      // result has the inner index and outer index transposed, i.e.
      // contract<2,1>(H,F) implies
      // T_{IJLk} = (H_{IJMN} F_{mM}) \delta_{mL} \delta_{Nk}
      // rather than T_{IJkL} (the desired result).
      // So, in effect, contraction of the 3rd (inner) index with F as the
      // second argument results in its transposition with respect to its
      // adjacent neighbor. This is due to the position of the argument F,
      // leading to the free index being on the right hand side of the result.
      // However, given that we can do two transformations from the LHS of H
      // and two from the right we can undo the otherwise erroneous
      // swapping of the outer indices upon application of the second
      // sets of contractions.
      //
      // Note: Its significantly quicker (in 3d) to push forward
      // each index individually
      return contract<1, 1>(
        F, contract<1, 1>(F, contract<2, 1>(contract<2, 1>(H, F), F)));
    }



    template <int dim, typename Number>
    inline ::SymmetricTensor<4, dim, Number>
    transformation_contraction(const ::SymmetricTensor<4, dim, Number> &H,
                               const Tensor<2, dim, Number> &                 F)
    {
      // The first and last transformation operations respectively
      // break and recover the symmetry properties of the tensors.
      // We also want to perform a minimal number of operations here
      // and avoid some complications related to the transposition of
      // tensor indices when contracting inner indices using the contract()
      // function. (For an explanation of the contraction operations,
      // please see the note in the equivalent function for standard
      // Tensors.) So what we'll do here is manually perform the first
      // and last contractions that break/recover the tensor symmetries
      // on the inner indices, and use the contract() function only on
      // the outer indices.
      //
      // Note: Its significantly quicker (in 3d) to push forward
      // each index individually

      // Push forward (inner) index 1
      Tensor<4, dim, Number> tmp;
      for (unsigned int I = 0; I < dim; ++I)
        for (unsigned int j = 0; j < dim; ++j)
          for (unsigned int K = 0; K < dim; ++K)
            for (unsigned int L = 0; L < dim; ++L)
              for (unsigned int J = 0; J < dim; ++J)
                tmp[I][j][K][L] += F[j][J] * H[I][J][K][L];

      // Push forward (outer) indices 0 and 3
      tmp = contract<1, 0>(F, contract<3, 1>(tmp, F));

      // Push forward (inner) index 2
      ::SymmetricTensor<4, dim, Number> out;
      for (unsigned int i = 0; i < dim; ++i)
        for (unsigned int j = i; j < dim; ++j)
          for (unsigned int k = 0; k < dim; ++k)
            for (unsigned int l = k; l < dim; ++l)
              for (unsigned int K = 0; K < dim; ++K)
                out[i][j][k][l] += F[k][K] * tmp[i][j][K][l];

      return out;
    }
  } // namespace Physics
} // namespace internal



template <typename Number>
Tensor<2, 2, Number>
Physics::Transformations::Rotations::rotation_matrix_2d(const Number &angle)
{
  const double rotation[2][2] = {{std::cos(angle), -std::sin(angle)},
                                 {std::sin(angle), std::cos(angle)}};
  return Tensor<2, 2>(rotation);
}



template <typename Number>
Tensor<2, 3, Number>
Physics::Transformations::Rotations::rotation_matrix_3d(
  const Point<3, Number> &axis,
  const Number &          angle)
{
  Assert(std::abs(axis.norm() - 1.0) < 1e-9,
         ExcMessage("The supplied axial vector is not a unit vector."));
  const Number c              = std::cos(angle);
  const Number s              = std::sin(angle);
  const Number t              = 1. - c;
  const double rotation[3][3] = {{t * axis[0] * axis[0] + c,
                                  t * axis[0] * axis[1] - s * axis[2],
                                  t * axis[0] * axis[2] + s * axis[1]},
                                 {t * axis[0] * axis[1] + s * axis[2],
                                  t * axis[1] * axis[1] + c,
                                  t * axis[1] * axis[2] - s * axis[0]},
                                 {t * axis[0] * axis[2] - s * axis[1],
                                  t * axis[1] * axis[2] + s * axis[0],
                                  t * axis[2] * axis[2] + c}};
  return Tensor<2, 3, Number>(rotation);
}



template <int dim, typename Number>
inline Tensor<1, dim, Number>
Physics::Transformations::Contravariant::push_forward(
  const Tensor<1, dim, Number> &V,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(V, F);
}



template <int dim, typename Number>
inline Tensor<2, dim, Number>
Physics::Transformations::Contravariant::push_forward(
  const Tensor<2, dim, Number> &T,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(T, F);
}



template <int dim, typename Number>
inline SymmetricTensor<2, dim, Number>
Physics::Transformations::Contravariant::push_forward(
  const SymmetricTensor<2, dim, Number> &T,
  const Tensor<2, dim, Number> &         F)
{
  return internal::Physics::transformation_contraction(T, F);
}



template <int dim, typename Number>
inline Tensor<4, dim, Number>
Physics::Transformations::Contravariant::push_forward(
  const Tensor<4, dim, Number> &H,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(H, F);
}



template <int dim, typename Number>
inline SymmetricTensor<4, dim, Number>
Physics::Transformations::Contravariant::push_forward(
  const SymmetricTensor<4, dim, Number> &H,
  const Tensor<2, dim, Number> &         F)
{
  return internal::Physics::transformation_contraction(H, F);
}



template <int dim, typename Number>
inline Tensor<1, dim, Number>
Physics::Transformations::Contravariant::pull_back(
  const Tensor<1, dim, Number> &v,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(v, invert(F));
}



template <int dim, typename Number>
inline Tensor<2, dim, Number>
Physics::Transformations::Contravariant::pull_back(
  const Tensor<2, dim, Number> &t,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(t, invert(F));
}



template <int dim, typename Number>
inline SymmetricTensor<2, dim, Number>
Physics::Transformations::Contravariant::pull_back(
  const SymmetricTensor<2, dim, Number> &t,
  const Tensor<2, dim, Number> &         F)
{
  return internal::Physics::transformation_contraction(t, invert(F));
}



template <int dim, typename Number>
inline Tensor<4, dim, Number>
Physics::Transformations::Contravariant::pull_back(
  const Tensor<4, dim, Number> &h,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(h, invert(F));
}



template <int dim, typename Number>
inline SymmetricTensor<4, dim, Number>
Physics::Transformations::Contravariant::pull_back(
  const SymmetricTensor<4, dim, Number> &h,
  const Tensor<2, dim, Number> &         F)
{
  return internal::Physics::transformation_contraction(h, invert(F));
}



template <int dim, typename Number>
inline Tensor<1, dim, Number>
Physics::Transformations::Covariant::push_forward(
  const Tensor<1, dim, Number> &V,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(V, transpose(invert(F)));
}



template <int dim, typename Number>
inline Tensor<2, dim, Number>
Physics::Transformations::Covariant::push_forward(
  const Tensor<2, dim, Number> &T,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(T, transpose(invert(F)));
}



template <int dim, typename Number>
inline SymmetricTensor<2, dim, Number>
Physics::Transformations::Covariant::push_forward(
  const SymmetricTensor<2, dim, Number> &T,
  const Tensor<2, dim, Number> &         F)
{
  return internal::Physics::transformation_contraction(T, transpose(invert(F)));
}



template <int dim, typename Number>
inline Tensor<4, dim, Number>
Physics::Transformations::Covariant::push_forward(
  const Tensor<4, dim, Number> &H,
  const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(H, transpose(invert(F)));
}



template <int dim, typename Number>
inline SymmetricTensor<4, dim, Number>
Physics::Transformations::Covariant::push_forward(
  const SymmetricTensor<4, dim, Number> &H,
  const Tensor<2, dim, Number> &         F)
{
  return internal::Physics::transformation_contraction(H, transpose(invert(F)));
}



template <int dim, typename Number>
inline Tensor<1, dim, Number>
Physics::Transformations::Covariant::pull_back(const Tensor<1, dim, Number> &v,
                                               const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(v, transpose(F));
}



template <int dim, typename Number>
inline Tensor<2, dim, Number>
Physics::Transformations::Covariant::pull_back(const Tensor<2, dim, Number> &t,
                                               const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(t, transpose(F));
}



template <int dim, typename Number>
inline SymmetricTensor<2, dim, Number>
Physics::Transformations::Covariant::pull_back(
  const SymmetricTensor<2, dim, Number> &t,
  const Tensor<2, dim, Number> &         F)
{
  return internal::Physics::transformation_contraction(t, transpose(F));
}



template <int dim, typename Number>
inline Tensor<4, dim, Number>
Physics::Transformations::Covariant::pull_back(const Tensor<4, dim, Number> &h,
                                               const Tensor<2, dim, Number> &F)
{
  return internal::Physics::transformation_contraction(h, transpose(F));
}



template <int dim, typename Number>
inline SymmetricTensor<4, dim, Number>
Physics::Transformations::Covariant::pull_back(
  const SymmetricTensor<4, dim, Number> &h,
  const Tensor<2, dim, Number> &         F)
{
  return internal::Physics::transformation_contraction(h, transpose(F));
}



template <int dim, typename Number>
inline Tensor<1, dim, Number>
Physics::Transformations::Piola::push_forward(const Tensor<1, dim, Number> &V,
                                              const Tensor<2, dim, Number> &F)
{
  return Number(1.0 / determinant(F)) * Contravariant::push_forward(V, F);
}



template <int dim, typename Number>
inline Tensor<2, dim, Number>
Physics::Transformations::Piola::push_forward(const Tensor<2, dim, Number> &T,
                                              const Tensor<2, dim, Number> &F)
{
  return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
}



template <int dim, typename Number>
inline SymmetricTensor<2, dim, Number>
Physics::Transformations::Piola::push_forward(
  const SymmetricTensor<2, dim, Number> &T,
  const Tensor<2, dim, Number> &         F)
{
  return Number(1.0 / determinant(F)) * Contravariant::push_forward(T, F);
}



template <int dim, typename Number>
inline Tensor<4, dim, Number>
Physics::Transformations::Piola::push_forward(const Tensor<4, dim, Number> &H,
                                              const Tensor<2, dim, Number> &F)
{
  return Number(1.0 / determinant(F)) * Contravariant::push_forward(H, F);
}



template <int dim, typename Number>
inline SymmetricTensor<4, dim, Number>
Physics::Transformations::Piola::push_forward(
  const SymmetricTensor<4, dim, Number> &H,
  const Tensor<2, dim, Number> &         F)
{
  return Number(1.0 / determinant(F)) * Contravariant::push_forward(H, F);
}



template <int dim, typename Number>
inline Tensor<1, dim, Number>
Physics::Transformations::Piola::pull_back(const Tensor<1, dim, Number> &v,
                                           const Tensor<2, dim, Number> &F)
{
  return Number(determinant(F)) * Contravariant::pull_back(v, F);
}



template <int dim, typename Number>
inline Tensor<2, dim, Number>
Physics::Transformations::Piola::pull_back(const Tensor<2, dim, Number> &t,
                                           const Tensor<2, dim, Number> &F)
{
  return Number(determinant(F)) * Contravariant::pull_back(t, F);
}



template <int dim, typename Number>
inline SymmetricTensor<2, dim, Number>
Physics::Transformations::Piola::pull_back(
  const SymmetricTensor<2, dim, Number> &t,
  const Tensor<2, dim, Number> &         F)
{
  return Number(determinant(F)) * Contravariant::pull_back(t, F);
}



template <int dim, typename Number>
inline Tensor<4, dim, Number>
Physics::Transformations::Piola::pull_back(const Tensor<4, dim, Number> &h,
                                           const Tensor<2, dim, Number> &F)
{
  return Number(determinant(F)) * Contravariant::pull_back(h, F);
}



template <int dim, typename Number>
inline SymmetricTensor<4, dim, Number>
Physics::Transformations::Piola::pull_back(
  const SymmetricTensor<4, dim, Number> &h,
  const Tensor<2, dim, Number> &         F)
{
  return Number(determinant(F)) * Contravariant::pull_back(h, F);
}



template <int dim, typename Number>
inline Tensor<1, dim, Number>
Physics::Transformations::nansons_formula(const Tensor<1, dim, Number> &N,
                                          const Tensor<2, dim, Number> &F)
{
  return cofactor(F) * N;
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif