linearalgebrahelper.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
#include "concepts.hh"
#include "traits.hh"
 
#include <iosfwd>
#include <random>
#include <ranges>
 
#include <dune/common/tuplevector.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/istl/bvector.hh>
 
#include <Eigen/Core>
 
#include <autodiff/forward/dual/dual.hpp>
 
namespace Ikarus {
 
template <typename Derived>
auto orthonormalizeMatrixColumns(const Eigen::MatrixBase<Derived>& A) {
  // Gram Schmidt Ortho
  auto Q = A.eval();
 
  Q.col(0).normalize();
 
  for (int colIndex = 1; colIndex < Q.cols(); colIndex++) {
    Q.col(colIndex) -= Q.leftCols(colIndex) * (Q.leftCols(colIndex).transpose() * A.col(colIndex));
    Q.col(colIndex).normalize();
  }
 
  return Q;
}
 
template <typename ValueType>
auto viewAsFlatEigenVector(Dune::BlockVector<ValueType>& blockedVector) {
  Eigen::Map<Eigen::VectorX<typename ValueType::field_type>> vec(&blockedVector.begin()->begin().operator*(),
                                                                 blockedVector.size() * blockedVector[0].size());
 
  return vec;
}
 
template <typename ValueType>
auto viewAsFlatEigenVector(const Dune::BlockVector<ValueType>& blockedVector) {
  Eigen::Map<const Eigen::VectorX<typename ValueType::field_type>> vec(&blockedVector.begin()->begin().operator*(),
                                                                       blockedVector.size() * blockedVector[0].size());
 
  return vec;
}
 
template <typename ValueType>
auto viewAsEigenMatrixAsDynFixed(Dune::BlockVector<ValueType>& blockedVector) {
  Eigen::Map<Eigen::Matrix<typename ValueType::field_type, Eigen::Dynamic, ValueType::valueSize, Eigen::RowMajor>> vec(
      &blockedVector.begin()->begin().operator*(), blockedVector.size(), blockedVector[0].size());
 
  return vec;
}
 
template <typename ValueType>
auto viewAsEigenMatrixAsDynFixed(const Dune::BlockVector<ValueType>& blockedVector) {
  Eigen::Map<const Eigen::Matrix<typename ValueType::field_type, Eigen::Dynamic, ValueType::valueSize, Eigen::RowMajor>>
      vec(&blockedVector.begin()->begin().operator*(), blockedVector.size(), blockedVector[0].size());
 
  return vec;
}
 
template <typename ValueType>
auto viewAsEigenMatrixFixedDyn(Dune::BlockVector<ValueType>& blockedVector) {
  Eigen::Map<Eigen::Matrix<typename ValueType::field_type, ValueType::valueSize, Eigen::Dynamic>> vec(
      &blockedVector.begin()->begin().operator*(), blockedVector[0].size(), blockedVector.size());
 
  return vec;
}
 
template <typename ValueType>
auto viewAsEigenMatrixFixedDyn(const Dune::BlockVector<ValueType>& blockedVector) {
  Eigen::Map<const Eigen::Matrix<typename ValueType::field_type, ValueType::valueSize, Eigen::Dynamic>> vec(
      &blockedVector.begin()->begin().operator*(), blockedVector[0].size(), blockedVector.size());
 
  return vec;
}
 
template <typename Type>
size_t correctionSize(const Dune::BlockVector<Type>& a)
requires requires { Type::correctionSize; }
{
  return a.size() * Type::correctionSize;
}
 
template <typename Type>
size_t valueSize(const Dune::BlockVector<Type>& a)
requires requires { Type::valueSize; }
{
  return a.size() * Type::valueSize;
}
 
template <typename Type, typename Derived>
Dune::BlockVector<Type>& operator+=(Dune::BlockVector<Type>& a, const Eigen::MatrixBase<Derived>& b)
requires(Ikarus::Concepts::AddAssignAble<Type, decltype(b.template segment<Type::correctionSize>(0))> and
         requires() { Type::correctionSize; })
{
  assert(correctionSize(a) == static_cast<size_t>(b.size()) && " The passed vector has wrong size");
  for (auto i = 0U; i < a.size(); ++i)
    a[i] += b.template segment<Type::correctionSize>(i * Type::correctionSize);
  return a;
}
 
template <typename Type, typename Derived>
Dune::BlockVector<Type>& operator-=(Dune::BlockVector<Type>& a, const Eigen::MatrixBase<Derived>& b)
requires(Ikarus::Concepts::AddAssignAble<Type, decltype(b.template segment<Type::correctionSize>(0))> and
         requires() { Type::correctionSize; })
{
  return a += (-b);
}
 
template <typename... Types, typename Derived>
Dune::TupleVector<Types...>& operator+=(Dune::TupleVector<Types...>& a, const Eigen::MatrixBase<Derived>& b) {
  using namespace Dune::Indices;
  size_t posStart = 0;
  Dune::Hybrid::forEach(Dune::Hybrid::integralRange(Dune::index_constant<sizeof...(Types)>()), [&](const auto i) {
    const size_t size = correctionSize(a[i]);
    a[i] += b(Eigen::seqN(posStart, size));
    posStart += size;
  });
 
  return a;
}
 
template <typename ManifoldPoint, typename Derived>
Dune::BlockVector<ManifoldPoint>& addInEmbedding(Dune::BlockVector<ManifoldPoint>& a,
                                                 const Eigen::MatrixBase<Derived>& b)
requires(Ikarus::Concepts::AddAssignAble<ManifoldPoint, decltype(b.template segment<ManifoldPoint::valueSize>(0))> and
         requires() { ManifoldPoint::valueSize; })
{
  assert(valueSize(a) == static_cast<size_t>(b.size()) && " The passed vector has wrong size");
  for (auto i = 0U; i < a.size(); ++i)
    a[i].addInEmbedding(b.template segment<ManifoldPoint::valueSize>(i * ManifoldPoint::valueSize));
  return a;
}
 
template <typename Derived>
requires(!std::floating_point<Derived>)
auto norm(const Eigen::MatrixBase<Derived>& v) {
  return v.norm();
}
 
auto norm(const std::floating_point auto& v) { return std::abs(v); }
 
template <typename Scalar, int size>
auto operator*(const Eigen::DiagonalMatrix<Scalar, size>& a, const Eigen::DiagonalMatrix<Scalar, size>& b) {
  return (a.diagonal().cwiseProduct(b.diagonal())).asDiagonal();
}
 
template <typename Scalar, int size>
auto operator+=(Eigen::DiagonalMatrix<Scalar, size>& a, const Eigen::DiagonalMatrix<Scalar, size>& b) {
  a.diagonal() += b.diagonal();
  return a;
}
 
template <typename Derived, typename Scalar, int size>
auto operator+(const Eigen::MatrixBase<Derived>& a, const Eigen::DiagonalMatrix<Scalar, size>& b) {
  auto c = a.derived().eval();
  c.diagonal() += b.diagonal();
  return c;
}
 
template <typename Derived, typename Scalar, int size>
auto operator+(const Eigen::DiagonalMatrix<Scalar, size>& a, const Eigen::MatrixBase<Derived>& b) {
  return b + a;
}
 
template <typename Scalar, int size>
auto operator-(const Eigen::DiagonalMatrix<Scalar, size>& a) {
  return (-a.diagonal()).asDiagonal();
}
 
template <typename Derived, typename Derived2>
auto operator+(const Eigen::MatrixBase<Derived>& a, const Eigen::DiagonalWrapper<Derived2>& b) {
  auto c = a.derived().eval();
  c.diagonal() += b.diagonal();
  return c;
}
 
template <typename Derived, typename Derived2>
auto operator+(const Eigen::DiagonalWrapper<Derived>& a, const Eigen::MatrixBase<Derived2>& b) {
  return b + a;
}
 
template <typename Scalar, int size>
std::ostream& operator<<(std::ostream& os, const Eigen::DiagonalMatrix<Scalar, size>& a) {
  os << Eigen::Matrix<Scalar, size, size>(a);
  return os;
}
 
template <typename Derived>
Derived sym(const Eigen::MatrixBase<Derived>& A) {
  return 0.5 * (A + A.transpose());
}
 
template <typename Derived>
Derived skew(const Eigen::MatrixBase<Derived>& A) {
  return 0.5 * (A - A.transpose());
}
 
template <typename Derived>
void printForMaple(const Eigen::EigenBase<Derived>& A) {
  Eigen::IOFormat mapleFmt(Eigen::FullPrecision, 0, ", ", "|\n", "<", ">", "<", ">");
  if constexpr (std::convertible_to<Derived, const Eigen::MatrixBase<Derived>&>) {
    std::cout << "\n" << A.derived().format(mapleFmt) << std::endl;
  } else { // branch for Eigen::DiagonalMatrix
    using Scalar = typename Derived::Scalar;
    using namespace Eigen;
    constexpr int diag_size = EIGEN_SIZE_MIN_PREFER_DYNAMIC(Derived::RowsAtCompileTime, Derived::ColsAtCompileTime);
    std::cout << "\n"
              << Eigen::Matrix<Scalar, diag_size, diag_size>(A.derived().diagonal().asDiagonal()).format(mapleFmt)
              << std::endl;
  }
}
 
template <typename FieldVectorT>
auto createRandomVector(typename FieldVectorT::value_type lower = -1, typename FieldVectorT::value_type upper = 1) {
  std::random_device rd;
  std::mt19937 mt(rd());
  std::uniform_real_distribution<typename FieldVectorT::value_type> dist(lower, upper);
  auto rand = [&dist, &mt]() { return dist(mt); };
  FieldVectorT vec;
  std::generate(vec.begin(), vec.end(), rand);
  return vec;
}
 
template <typename ScalarType>
Eigen::Matrix<ScalarType, 3, 3> skew(const Eigen::Vector<ScalarType, 3>& a) {
  Eigen::Matrix<ScalarType, 3, 3> A;
  A << 0, -a(2), a(1), a(2), 0, -a(0), -a(1), a(0), 0;
  return A;
}
 
namespace Impl {
  constexpr std::tuple<std::array<std::array<int, 2>, 1>, std::array<std::array<int, 2>, 3>,
                       std::array<std::array<int, 2>, 6>>
      voigtIndices = {{{{0, 0}}}, {{{0, 0}, {1, 1}, {0, 1}}}, {{{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}}}};
}
 
template <int dim>
constexpr auto voigtNotationContainer = std::get<dim - 1>(Impl::voigtIndices);
 
template <typename Derived, size_t sizeOfCondensedIndices>
auto staticCondensation(const Eigen::MatrixBase<Derived>& E,
                        const std::array<size_t, sizeOfCondensedIndices>& indices) {
  constexpr size_t colsFull = std::remove_cvref_t<Derived>::ColsAtCompileTime;
  constexpr size_t rowsFull = std::remove_cvref_t<Derived>::RowsAtCompileTime;
  static_assert(colsFull == rowsFull, "staticCondensation only allowed for square matrices");
  std::array<size_t, rowsFull - sizeOfCondensedIndices> remainingIndices{};
  std::ranges::set_difference(std::ranges::iota_view(size_t(0), size_t(colsFull)), indices, remainingIndices.begin());
 
  auto K11 = E(remainingIndices, remainingIndices);
  auto K12 = E(indices, remainingIndices);
  auto K22 = E(indices, indices);
 
  return (K11 - K12.transpose() * K22.inverse() * K12).eval();
}
 
template <typename Derived, size_t sizeOfCondensedIndices>
auto reduceMatrix(const Eigen::MatrixBase<Derived>& E, const std::array<size_t, sizeOfCondensedIndices>& indices) {
  constexpr size_t colsFull = std::remove_cvref_t<Derived>::ColsAtCompileTime;
  constexpr size_t rowsFull = std::remove_cvref_t<Derived>::RowsAtCompileTime;
  static_assert(colsFull == rowsFull, "reduceMatrix only allowed for square matrices");
  std::array<size_t, rowsFull - sizeOfCondensedIndices> remainingIndices{};
  std::ranges::set_difference(std::ranges::iota_view(size_t(0), size_t(colsFull)), indices, remainingIndices.begin());
 
  auto K11 = E(remainingIndices, remainingIndices);
 
  return K11.eval();
}
 
template <typename Derived, size_t sizeOfRemovedCols>
auto removeCol(const Eigen::MatrixBase<Derived>& E, const std::array<size_t, sizeOfRemovedCols>& indices) {
  constexpr size_t colsFull = std::remove_cvref_t<Derived>::ColsAtCompileTime;
  constexpr size_t rowsFull = std::remove_cvref_t<Derived>::RowsAtCompileTime;
  static_assert(colsFull == 1);
 
  std::array<size_t, rowsFull - sizeOfRemovedCols> remainingIndices{};
  std::ranges::set_difference(std::ranges::iota_view(size_t(0), size_t(rowsFull)), indices, remainingIndices.begin());
 
  return (E(remainingIndices)).eval();
}
 
template <typename ST, typename MaterialImpl>
auto toVoigtAndMaybeReduce(const Eigen::Matrix<ST, 3, 3>& E, [[maybe_unused]] const MaterialImpl& material,
                           bool isStrain = true) {
  if constexpr (!MaterialImpl::isReduced)
    return toVoigt(E, isStrain);
  else {
    auto ev = toVoigt(E, isStrain);
    static_assert(decltype(ev)::RowsAtCompileTime == 6 and decltype(ev)::ColsAtCompileTime == 1);
 
    auto evRed = removeCol(ev, MaterialImpl::fixedVoigtIndices);
    static_assert(decltype(evRed)::RowsAtCompileTime == 6 - MaterialImpl::fixedVoigtIndices.size() and
                  decltype(evRed)::ColsAtCompileTime == 1);
    return evRed;
  }
}
 
template <typename M, typename Derived>
decltype(auto) enlargeIfReduced(const Eigen::MatrixBase<Derived>& E) {
  if constexpr (!Concepts::EigenVector6<Derived> and Concepts::EigenVector<Derived>) {
    static_assert(M::isReduced, "You should only end up here, if your material is reduced");
 
    auto freeindices = M::MaterialImpl::freeVoigtIndices;
    auto p_E         = Eigen::Vector<typename M::MaterialImpl::ScalarType, 6>::Zero().eval();
    for (int ri = 0; auto i : freeindices)
      p_E(i) = E(ri++);
    return p_E;
 
  } else if constexpr (Concepts::EigenMatrix66<Derived> or Concepts::EigenMatrix33<Derived> or
                       Concepts::EigenVector6<Derived>) {
    return E.derived();
  } else {
    static_assert(M::isReduced, "You should only end up here, if your material is reduced");
 
    auto freeindices = M::MaterialImpl::freeVoigtIndices;
    auto p_E         = Eigen::Matrix<typename M::MaterialImpl::ScalarType, 6, 6>::Zero().eval();
    for (int ri = 0; auto i : freeindices) {
      for (int rj = 0; auto j : freeindices)
        p_E(i, j) = E(ri, rj++);
      ++ri;
    }
    return p_E;
  }
}
 
} // namespace Ikarus