tensorutils.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <numeric>
#include <ranges>
#include <unsupported/Eigen/CXX11/Tensor>
 
#include <dune/common/promotiontraits.hh>
 
#include <ikarus/utils/concepts.hh>
#include <ikarus/utils/math.hh>
namespace Ikarus {
 
template <typename Derived, typename T, auto rank>
Eigen::Tensor<typename Derived::Scalar, rank> tensorView(const Eigen::EigenBase<Derived>& matrix,
                                                         const std::array<T, rank>& dims) {
  return Eigen::TensorMap<const Eigen::TensorFixedSize<
      const typename Derived::Scalar, Eigen::Sizes<Derived::RowsAtCompileTime, Derived::ColsAtCompileTime>>>(
      matrix.derived().eval().data(), dims);
}
 
auto dyadic(const auto& A_ij, const auto& B_kl) {
  Eigen::array<Eigen::IndexPair<long>, 0> empty_index_list = {};
  return A_ij.contract(B_kl, empty_index_list).eval();
}
 
template <typename ST, int size>
auto dyadic(const Eigen::Vector<ST, size>& a, const Eigen::Vector<ST, size>& b) {
  return (a * b.transpose()).eval();
}
 
template <typename ScalarType = double, int dim = 3>
auto symmetricIdentityFourthOrder() {
  Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim, dim, dim>> idTensor;
  for (int i = 0; i < dim; ++i)
    for (int j = 0; j < dim; ++j)
      for (int k = 0; k < dim; ++k)
        for (int l = 0; l < dim; ++l)
          idTensor(i, j, k, l) = 0.5 * ((i == k) * (j == l) + (i == l) * (j == k));
  return idTensor;
}
 
template <typename ScalarType = double, int dim = 3>
auto symmetricFourthOrder(const auto& A, const auto& B) {
  Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim, dim, dim>> res;
  for (int i = 0; i < dim; ++i)
    for (int j = 0; j < dim; ++j)
      for (int k = 0; k < dim; ++k)
        for (int l = 0; l < dim; ++l)
          res(i, j, k, l) = 0.5 * (A(i, k) * B(j, l) + A(i, l) * B(j, k));
  return res;
}
 
template <typename ScalarType = double, int dim = 3>
auto identityFourthOrder() {
  Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim, dim, dim>> idTensor;
  idTensor.setZero();
  for (int i = 0; i < dim; ++i)
    for (int k = 0; k < dim; ++k)
      idTensor(i, i, k, k) = 1.0;
  return idTensor;
}
 
template <typename AType, typename BType>
auto fourthOrderIKJL(const Eigen::MatrixBase<AType>& A, const Eigen::MatrixBase<BType>& B) {
  static_assert(AType::RowsAtCompileTime == BType::RowsAtCompileTime);
  static_assert(AType::ColsAtCompileTime == BType::ColsAtCompileTime);
  using ScalarResultType = typename Dune::PromotionTraits<typename AType::Scalar, typename BType::Scalar>::PromotedType;
  constexpr int dim      = AType::RowsAtCompileTime;
  Eigen::TensorFixedSize<ScalarResultType, Eigen::Sizes<dim, dim, dim, dim>> res;
  for (int i = 0; i < dim; ++i)
    for (int j = 0; j < dim; ++j)
      for (int k = 0; k < dim; ++k)
        for (int l = 0; l < dim; ++l)
          res(i, j, k, l) = A(i, k) * B(j, l);
  return res;
}
 
template <typename ScalarType, long int dim>
auto symTwoSlots(const Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim, dim, dim>>& t,
                 const std::array<size_t, 2>& slots) {
  std::array<size_t, 4> shuffleSlot;
  std::iota(shuffleSlot.begin(), shuffleSlot.end(), 0);    // create 0,1,2,3 array
  std::swap(shuffleSlot[slots[0]], shuffleSlot[slots[1]]); // swap  the given slots
  return (0.5 * (t + t.shuffle(shuffleSlot))).eval();
}
 
constexpr Eigen::Index toVoigt(Eigen::Index i, Eigen::Index j) noexcept {
  if (i == j) // _00 -> 0, _11 -> 1,  _22 -> 2
    return i;
  if ((i == 1 and j == 2) or (i == 2 and j == 1)) // _12 and _21 --> 3
    return 3;
  if ((i == 0 and j == 2) or (i == 2 and j == 0)) // _02 and _20 --> 4
    return 4;
  if ((i == 0 and j == 1) or (i == 1 and j == 0)) // _01 and _10 --> 5
    return 5;
  assert(i < 3 and j < 3 && "For Voigt notation the indices need to be 0,1 or 2.");
  __builtin_unreachable();
}
 
template <typename ScalarType = double>
Eigen::Matrix<ScalarType, 6, 6> toVoigt(const Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3, 3, 3>>& ft) {
  Eigen::Matrix<ScalarType, 6, 6> mat;
  for (Eigen::Index i = 0; i < 3; ++i)
    for (Eigen::Index j = 0; j < 3; ++j)
      for (Eigen::Index k = 0; k < 3; ++k)
        for (Eigen::Index l = 0; l < 3; ++l)
          mat(toVoigt(i, j), toVoigt(k, l)) = ft(i, j, k, l);
  return mat;
}
 
template <typename ST, int size, int Options, int maxSize>
requires((size > 0 and size <= 3) or (maxSize > 0 and maxSize <= 3 and size == Eigen::Dynamic))
auto toVoigt(const Eigen::Matrix<ST, size, size, Options, maxSize, maxSize>& E, bool isStrain = true) {
  constexpr bool isFixedSized   = (size != Eigen::Dynamic);
  const ST possibleStrainFactor = isStrain ? 2.0 : 1.0;
 
  const size_t inputSize = isFixedSized ? size : E.rows();
  auto EVoigt            = [&]() {
    if constexpr (isFixedSized) {
      Eigen::Vector<ST, (size * (size + 1)) / 2> EVoigt;
      EVoigt.template head<size>() = E.diagonal();
      return EVoigt;
    } else {
      Eigen::Matrix<ST, Eigen::Dynamic, 1, Options, 6, 1> EVoigt;
      EVoigt.resize((inputSize * (inputSize + 1)) / 2);
      EVoigt.template head(inputSize) = E.diagonal();
      return EVoigt;
    }
  }();
 
  if (inputSize == 2)
    EVoigt(2) = E(0, 1) * possibleStrainFactor;
  else if (inputSize == 3) {
    EVoigt(inputSize)     = E(1, 2) * possibleStrainFactor;
    EVoigt(inputSize + 1) = E(0, 2) * possibleStrainFactor;
    EVoigt(inputSize + 2) = E(0, 1) * possibleStrainFactor;
  }
  return EVoigt;
}
 
template <typename ST, int size, int Options, int maxSize>
requires((size == 1 or size == 3 or size == 6) or
         ((maxSize == 1 or maxSize == 3 or maxSize == 6) and size == Eigen::Dynamic))
auto fromVoigt(const Eigen::Matrix<ST, size, 1, Options, maxSize, 1>& EVoigt, bool isStrain = true) {
  constexpr bool isFixedSized   = (size != Eigen::Dynamic);
  const ST possibleStrainFactor = isStrain ? 0.5 : 1.0;
 
  const size_t inputSize = isFixedSized ? size : EVoigt.size();
  const size_t matrixSize =
      isFixedSized ? (-1 + ct_sqrt(1 + 8 * size)) / 2 : (-1 + static_cast<int>(std::sqrt(1 + 8 * inputSize))) / 2;
 
  auto E = [&]() {
    if constexpr (isFixedSized) {
      Eigen::Matrix<ST, matrixSize, matrixSize> E;
      E.diagonal() = EVoigt.template head<matrixSize>();
      return E;
    } else {
      Eigen::Matrix<ST, Eigen::Dynamic, Eigen::Dynamic, Options, 3, 3> E;
      E.resize(matrixSize, matrixSize);
      E.diagonal() = EVoigt.template head(matrixSize);
      return E;
    }
  }();
 
  if (matrixSize == 2) {
    E(0, 1) = E(1, 0) = EVoigt(2) * possibleStrainFactor;
  } else if (matrixSize == 3) {
    E(2, 1) = E(1, 2) = EVoigt(matrixSize) * possibleStrainFactor;
    E(2, 0) = E(0, 2) = EVoigt(matrixSize + 1) * possibleStrainFactor;
    E(1, 0) = E(0, 1) = EVoigt(matrixSize + 2) * possibleStrainFactor;
  }
 
  return E;
}
 
constexpr std::array<size_t, 2> fromVoigt(size_t i) {
  if (i < 3) // _00 -> 0, _11 -> 1,  _22 -> 2
    return {i, i};
  else if (i == 3)
    return {1, 2};
  else if (i == 4)
    return {0, 2};
  else if (i == 5)
    return {0, 1};
  else {
    assert(i < 6 && "For Voigt notation the indices need to be 0 and 5.");
    __builtin_unreachable();
  }
}
 
template <typename ScalarType>
auto fromVoigt(const Eigen::Matrix<ScalarType, 6, 6>& CVoigt) {
  Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3, 3, 3>> C;
  // size_t iR=0,jR=0;
  for (size_t i = 0; i < 6; ++i) {
    for (size_t j = 0; j < 6; ++j) {
      auto firstIndices                                                       = fromVoigt(i);
      auto secondIndices                                                      = fromVoigt(j);
      C(firstIndices[0], firstIndices[1], secondIndices[0], secondIndices[1]) = CVoigt(i, j);
    }
  }
  return C;
}
 
template <typename Geometry>
Eigen::Matrix3d calcTransformationMatrix2D(const Geometry& geometry) {
  const auto& referenceElement = Dune::ReferenceElements<double, 2>::general(geometry.type());
  const auto quadPos0          = referenceElement.position(0, 0);
 
  const auto jacobianinvT0 = toEigen(geometry.jacobianInverseTransposed(quadPos0));
  const auto detJ0         = geometry.integrationElement(quadPos0);
 
  auto jaco = (jacobianinvT0).inverse().eval();
  auto J11  = jaco(0, 0);
  auto J12  = jaco(0, 1);
  auto J21  = jaco(1, 0);
  auto J22  = jaco(1, 1);
 
  Eigen::Matrix3d T0{
      {      J11 * J11,       J12 * J12,             J11 * J12},
      {      J21 * J21,       J22 * J22,             J21 * J22},
      {2.0 * J11 * J21, 2.0 * J12 * J22, J21 * J12 + J11 * J22}
  };
 
  return T0.inverse() * detJ0;
}
 
template <typename Geometry>
Eigen::Matrix<double, 6, 6> calcTransformationMatrix3D(const Geometry& geometry) {
  const auto& referenceElement = Dune::ReferenceElements<double, 3>::general(geometry.type());
  const auto quadPos0          = referenceElement.position(0, 0);
 
  const auto jacobianinvT0 = toEigen(geometry.jacobianInverseTransposed(quadPos0));
  const auto detJ0         = geometry.integrationElement(quadPos0);
 
  auto jaco = (jacobianinvT0).inverse().eval();
  auto J11  = jaco(0, 0);
  auto J12  = jaco(0, 1);
  auto J13  = jaco(0, 2);
  auto J21  = jaco(1, 0);
  auto J22  = jaco(1, 1);
  auto J23  = jaco(1, 2);
  auto J31  = jaco(2, 0);
  auto J32  = jaco(2, 1);
  auto J33  = jaco(2, 2);
 
  // clang-format off
  Eigen::Matrix<double, 6, 6> T0  {
    {J11 * J11,       J12 * J12,       J13 * J13,             J11 * J12,             J11 * J13,             J12 * J13},
    {J21 * J21,       J22 * J22,       J23 * J23,             J21 * J22,             J21 * J23,             J22 * J23},
    {J31 * J31,       J32 * J32,       J33 * J33,             J31 * J32,             J31 * J33,             J32 * J33},
    {2.0 * J11 * J21, 2.0 * J12 * J22, 2.0 * J13 * J23, J11 * J22 + J21 * J12, J11 * J23 + J21 * J13, J12 * J23 + J22 * J13},
    {2.0 * J11 * J31, 2.0 * J12 * J32, 2.0 * J13 * J33, J11 * J32 + J31 * J12, J11 * J33 + J31 * J13, J12 * J33 + J32 * J13},
    {2.0 * J31 * J21, 2.0 * J32 * J22, 2.0 * J33 * J23, J31 * J22 + J21 * J32, J31 * J23 + J21 * J33, J32 * J23 + J22 * J33}
  };
  // clang-format on
 
  return T0.inverse() * detJ0;
}
 
} // namespace Ikarus