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 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 Material, typename Derived>
  decltype(auto) enlargeIfReduced(const Eigen::MatrixBase<Derived>& E) {
    if constexpr (!Concepts::EigenVector6<Derived> and Concepts::EigenVector<Derived>) {
      static_assert(Material::isReduced, "You should only end up here, if your material is reduced");
 
      auto freeindices = Material::MaterialImplType::freeVoigtIndices;
      auto p_E         = Eigen::Vector<typename Material::MaterialImplType::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(Material::isReduced, "You should only end up here, if your material is reduced");
 
      auto freeindices = Material::MaterialImpl::freeVoigtIndices;
      auto p_E         = Eigen::Matrix<typename Material::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