transformations.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2016 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_transformations_h
#define dealii_transformations_h
 
#include <deal.II/base/config.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 Tensor<1, 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);
 
    template <int dim, typename Number>
    Tensor<1, dim, Number>
    basis_transformation(const Tensor<1, dim, Number> &V,
                         const Tensor<2, dim, Number> &B);
 
    template <int dim, typename Number>
    Tensor<2, dim, Number>
    basis_transformation(const Tensor<2, dim, Number> &T,
                         const Tensor<2, dim, Number> &B);
 
    template <int dim, typename Number>
    SymmetricTensor<2, dim, Number>
    basis_transformation(const SymmetricTensor<2, dim, Number> &T,
                         const Tensor<2, dim, Number>          &B);
 
    template <int dim, typename Number>
    Tensor<4, dim, Number>
    basis_transformation(const Tensor<4, dim, Number> &H,
                         const Tensor<2, dim, Number> &B);
 
    template <int dim, typename Number>
    SymmetricTensor<4, dim, Number>
    basis_transformation(const SymmetricTensor<4, dim, Number> &H,
                         const Tensor<2, dim, Number>          &B);
 
  } // namespace Transformations
} // namespace Physics
 
 
 
#ifndef DOXYGEN
 
 
 
template <typename Number>
Tensor<2, 2, Number>
Physics::Transformations::Rotations::rotation_matrix_2d(const Number &angle)
{
  // Make things work with AD types
  using std::cos;
  using std::sin;
 
  const Number rotation[2][2] = {{cos(angle), -sin(angle)},
                                 {sin(angle), cos(angle)}};
  return Tensor<2, 2>(rotation);
}
 
 
 
template <typename Number>
Tensor<2, 3, Number>
Physics::Transformations::Rotations::rotation_matrix_3d(
  const Tensor<1, 3, Number> &axis,
  const Number               &angle)
{
  // Make things work with AD types
  using std::abs;
  using std::cos;
  using std::sin;
 
  Assert(abs(axis.norm() - 1.0) < 1e-9,
         ExcMessage("The supplied axial vector is not a unit vector."));
  const Number c              = cos(angle);
  const Number s              = sin(angle);
  const Number t              = 1. - c;
  const Number 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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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 Physics::Transformations::basis_transformation(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;
}
 
 
template <int dim, typename Number>
inline Tensor<1, dim, Number>
Physics::Transformations::basis_transformation(const Tensor<1, dim, Number> &V,
                                               const Tensor<2, dim, Number> &B)
{
  return contract<1, 0>(B, V);
}
 
 
 
template <int dim, typename Number>
inline Tensor<2, dim, Number>
Physics::Transformations::basis_transformation(const Tensor<2, dim, Number> &T,
                                               const Tensor<2, dim, Number> &B)
{
  return contract<1, 0>(B, contract<1, 1>(T, B));
}
 
 
 
template <int dim, typename Number>
inline SymmetricTensor<2, dim, Number>
Physics::Transformations::basis_transformation(
  const SymmetricTensor<2, dim, Number> &T,
  const Tensor<2, dim, Number>          &B)
{
  Tensor<2, dim, Number> tmp_1;
  for (unsigned int i = 0; i < dim; ++i)
    for (unsigned int J = 0; J < dim; ++J)
      // Loop over I but complex.h defines a macro I, so use I_ instead
      for (unsigned int I_ = 0; I_ < dim; ++I_)
        tmp_1[i][J] += B[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] += B[j][J] * tmp_1[i][J];
 
  return out;
}
 
 
 
template <int dim, typename Number>
inline Tensor<4, dim, Number>
Physics::Transformations::basis_transformation(const Tensor<4, dim, Number> &H,
                                               const Tensor<2, dim, Number> &B)
{
  // 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>(
    B, contract<1, 1>(B, contract<2, 1>(contract<2, 1>(H, B), B)));
}
 
 
 
template <int dim, typename Number>
inline SymmetricTensor<4, dim, Number>
Physics::Transformations::basis_transformation(
  const SymmetricTensor<4, dim, Number> &H,
  const Tensor<2, dim, Number>          &B)
{
  // 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;
  // Loop over I but complex.h defines a macro I, so use I_ instead
  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] += B[j][J] * H[I_][J][K][L];
 
  // Push forward (outer) indices 0 and 3
  tmp = contract<1, 0>(B, contract<3, 1>(tmp, B));
 
  // 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] += B[k][K] * tmp[i][j][K][l];
 
  return out;
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif