version 0.4.4
linearalgebrahelper.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers ikarus@ibb.uni-stuttgart.de
2// SPDX-License-Identifier: LGPL-3.0-or-later
3
9#pragma once
10#include <iosfwd>
11#include <random>
12#include <ranges>
13
14#include <dune/common/tuplevector.hh>
15#include <dune/functions/functionspacebases/lagrangebasis.hh>
16#include <dune/istl/bvector.hh>
17
18#include <Eigen/Core>
19
20#include <autodiff/forward/dual/dual.hpp>
21
23
24namespace Ikarus {
25
33template <typename Derived>
34auto orthonormalizeMatrixColumns(const Eigen::MatrixBase<Derived>& A) {
35 // Gram Schmidt Ortho
36 auto Q = A.eval();
37
38 Q.col(0).normalize();
39
40 for (int colIndex = 1; colIndex < Q.cols(); colIndex++) {
41 Q.col(colIndex) -= Q.leftCols(colIndex) * (Q.leftCols(colIndex).transpose() * A.col(colIndex));
42 Q.col(colIndex).normalize();
43 }
44
45 return Q;
46}
47
55template <typename ValueType>
56auto viewAsFlatEigenVector(Dune::BlockVector<ValueType>& blockedVector) {
57 Eigen::Map<Eigen::VectorX<typename ValueType::field_type>> vec(&blockedVector.begin()->begin().operator*(),
58 blockedVector.size() * blockedVector[0].size());
59
60 return vec;
61}
62
70template <typename ValueType>
71auto viewAsFlatEigenVector(const Dune::BlockVector<ValueType>& blockedVector) {
72 Eigen::Map<const Eigen::VectorX<typename ValueType::field_type>> vec(&blockedVector.begin()->begin().operator*(),
73 blockedVector.size() * blockedVector[0].size());
74
75 return vec;
76}
77
86template <typename ValueType>
87auto viewAsEigenMatrixAsDynFixed(Dune::BlockVector<ValueType>& blockedVector) {
88 Eigen::Map<Eigen::Matrix<typename ValueType::field_type, Eigen::Dynamic, ValueType::valueSize, Eigen::RowMajor>> vec(
89 &blockedVector.begin()->begin().operator*(), blockedVector.size(), blockedVector[0].size());
90
91 return vec;
92}
93
102template <typename ValueType>
103auto viewAsEigenMatrixAsDynFixed(const Dune::BlockVector<ValueType>& blockedVector) {
104 Eigen::Map<const Eigen::Matrix<typename ValueType::field_type, Eigen::Dynamic, ValueType::valueSize, Eigen::RowMajor>>
105 vec(&blockedVector.begin()->begin().operator*(), blockedVector.size(), blockedVector[0].size());
106
107 return vec;
108}
109
118template <typename ValueType>
119auto viewAsEigenMatrixFixedDyn(Dune::BlockVector<ValueType>& blockedVector) {
120 Eigen::Map<Eigen::Matrix<typename ValueType::field_type, ValueType::valueSize, Eigen::Dynamic>> vec(
121 &blockedVector.begin()->begin().operator*(), blockedVector[0].size(), blockedVector.size());
122
123 return vec;
124}
125
133template <typename ValueType>
134auto viewAsEigenMatrixFixedDyn(const Dune::BlockVector<ValueType>& blockedVector) {
135 Eigen::Map<const Eigen::Matrix<typename ValueType::field_type, ValueType::valueSize, Eigen::Dynamic>> vec(
136 &blockedVector.begin()->begin().operator*(), blockedVector[0].size(), blockedVector.size());
137
138 return vec;
139}
140
148template <typename Type>
149size_t correctionSize(const Dune::BlockVector<Type>& a)
150requires requires { Type::correctionSize; }
151{
152 return a.size() * Type::correctionSize;
153}
154
162template <typename Type>
163size_t valueSize(const Dune::BlockVector<Type>& a)
164requires requires { Type::valueSize; }
165{
166 return a.size() * Type::valueSize;
167}
168
178template <typename Type, typename Derived>
179Dune::BlockVector<Type>& operator+=(Dune::BlockVector<Type>& a, const Eigen::MatrixBase<Derived>& b)
180requires(Ikarus::Concepts::AddAssignAble<Type, decltype(b.template segment<Type::correctionSize>(0))> and
181 requires() { Type::correctionSize; })
182{
183 assert(correctionSize(a) == static_cast<size_t>(b.size()) && " The passed vector has wrong size");
184 for (auto i = 0U; i < a.size(); ++i)
185 a[i] += b.template segment<Type::correctionSize>(i * Type::correctionSize);
186 return a;
187}
188
198template <typename Type, typename Derived>
199Dune::BlockVector<Type>& operator-=(Dune::BlockVector<Type>& a, const Eigen::MatrixBase<Derived>& b)
200requires(Ikarus::Concepts::AddAssignAble<Type, decltype(b.template segment<Type::correctionSize>(0))> and
201 requires() { Type::correctionSize; })
202{
203 return a += (-b);
204}
205
215template <typename... Types, typename Derived>
216Dune::TupleVector<Types...>& operator+=(Dune::TupleVector<Types...>& a, const Eigen::MatrixBase<Derived>& b) {
217 using namespace Dune::Indices;
218 size_t posStart = 0;
219 Dune::Hybrid::forEach(Dune::Hybrid::integralRange(Dune::index_constant<sizeof...(Types)>()), [&](const auto i) {
220 const size_t size = correctionSize(a[i]);
221 a[i] += b(Eigen::seqN(posStart, size));
222 posStart += size;
223 });
224
225 return a;
226}
227
237template <typename ManifoldPoint, typename Derived>
238Dune::BlockVector<ManifoldPoint>& addInEmbedding(Dune::BlockVector<ManifoldPoint>& a,
239 const Eigen::MatrixBase<Derived>& b)
240requires(Ikarus::Concepts::AddAssignAble<ManifoldPoint, decltype(b.template segment<ManifoldPoint::valueSize>(0))> and
241 requires() { ManifoldPoint::valueSize; })
242{
243 assert(valueSize(a) == static_cast<size_t>(b.size()) && " The passed vector has wrong size");
244 for (auto i = 0U; i < a.size(); ++i)
245 a[i].addInEmbedding(b.template segment<ManifoldPoint::valueSize>(i * ManifoldPoint::valueSize));
246 return a;
247}
248
256template <typename Derived>
257requires(!std::floating_point<Derived>)
258auto norm(const Eigen::MatrixBase<Derived>& v) {
259 return v.norm();
260}
261
270template <typename Derived>
271requires(!std::floating_point<Derived>)
272auto floatingPointNorm(const Eigen::MatrixBase<Derived>& v) {
274 return v.template cast<autodiff::detail::NumericType<typename Derived::Scalar>>().norm();
275 else
276 return v.norm();
277}
278
285auto norm(const std::floating_point auto& v) { return std::abs(v); }
286
293auto floatingPointNorm(const std::floating_point auto& v) { return std::abs(v); }
294
304template <typename Scalar, int size>
305auto operator*(const Eigen::DiagonalMatrix<Scalar, size>& a, const Eigen::DiagonalMatrix<Scalar, size>& b) {
306 return (a.diagonal().cwiseProduct(b.diagonal())).asDiagonal();
307}
308
318template <typename Scalar, int size>
319auto operator+=(Eigen::DiagonalMatrix<Scalar, size>& a, const Eigen::DiagonalMatrix<Scalar, size>& b) {
320 a.diagonal() += b.diagonal();
321 return a;
322}
323
334template <typename Derived, typename Scalar, int size>
335auto operator+(const Eigen::MatrixBase<Derived>& a, const Eigen::DiagonalMatrix<Scalar, size>& b) {
336 auto c = a.derived().eval();
337 c.diagonal() += b.diagonal();
338 return c;
339}
340
351template <typename Derived, typename Scalar, int size>
352auto operator+(const Eigen::DiagonalMatrix<Scalar, size>& a, const Eigen::MatrixBase<Derived>& b) {
353 return b + a;
354}
355
364template <typename Scalar, int size>
365auto operator-(const Eigen::DiagonalMatrix<Scalar, size>& a) {
366 return (-a.diagonal()).asDiagonal();
367}
368
378template <typename Derived, typename Derived2>
379auto operator+(const Eigen::MatrixBase<Derived>& a, const Eigen::DiagonalWrapper<Derived2>& b) {
380 auto c = a.derived().eval();
381 c.diagonal() += b.diagonal();
382 return c;
383}
384
394template <typename Derived, typename Derived2>
395auto operator+(const Eigen::DiagonalWrapper<Derived>& a, const Eigen::MatrixBase<Derived2>& b) {
396 return b + a;
397}
398
408template <typename Scalar, int size>
409std::ostream& operator<<(std::ostream& os, const Eigen::DiagonalMatrix<Scalar, size>& a) {
410 os << Eigen::Matrix<Scalar, size, size>(a);
411 return os;
412}
413
421template <typename Derived>
422Derived sym(const Eigen::MatrixBase<Derived>& A) {
423 return 0.5 * (A + A.transpose());
424}
425
433template <typename Derived>
434Derived skew(const Eigen::MatrixBase<Derived>& A) {
435 return 0.5 * (A - A.transpose());
436}
437
444template <typename Derived>
446 Eigen::IOFormat mapleFmt(Eigen::FullPrecision, 0, ", ", "|\n", "<", ">", "<", ">");
447 if constexpr (std::convertible_to<Derived, const Eigen::MatrixBase<Derived>&>) {
448 std::cout << "\n" << A.derived().format(mapleFmt) << std::endl;
449 } else { // branch for Eigen::DiagonalMatrix
450 using Scalar = typename Derived::Scalar;
451 using namespace Eigen;
452 constexpr int diag_size = EIGEN_SIZE_MIN_PREFER_DYNAMIC(Derived::RowsAtCompileTime, Derived::ColsAtCompileTime);
453 std::cout << "\n"
454 << Eigen::Matrix<Scalar, diag_size, diag_size>(A.derived().diagonal().asDiagonal()).format(mapleFmt)
455 << std::endl;
456 }
457}
458
467template <typename FieldVectorT>
468auto createRandomVector(typename FieldVectorT::value_type lower = -1, typename FieldVectorT::value_type upper = 1) {
469 std::random_device rd;
470 std::mt19937 mt(rd());
471 std::uniform_real_distribution<typename FieldVectorT::value_type> dist(lower, upper);
472 auto rand = [&dist, &mt]() { return dist(mt); };
473 FieldVectorT vec;
474 std::generate(vec.begin(), vec.end(), rand);
475 return vec;
476}
477
485template <typename ScalarType>
486Eigen::Matrix<ScalarType, 3, 3> skew(const Eigen::Vector<ScalarType, 3>& a) {
487 Eigen::Matrix<ScalarType, 3, 3> A;
488 A << 0, -a(2), a(1), a(2), 0, -a(0), -a(1), a(0), 0;
489 return A;
490}
491
492namespace Impl {
493 constexpr std::tuple<std::array<std::array<int, 2>, 1>, std::array<std::array<int, 2>, 3>,
494 std::array<std::array<int, 2>, 6>>
495 voigtIndices = {{{{0, 0}}}, {{{0, 0}, {1, 1}, {0, 1}}}, {{{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}}}};
496}
497
506template <int dim>
507constexpr auto voigtNotationContainer = std::get<dim - 1>(Impl::voigtIndices);
508
521template <typename Derived, size_t sizeOfCondensedIndices>
522auto staticCondensation(const Eigen::MatrixBase<Derived>& E,
523 const std::array<size_t, sizeOfCondensedIndices>& indices) {
524 constexpr size_t colsFull = std::remove_cvref_t<Derived>::ColsAtCompileTime;
525 constexpr size_t rowsFull = std::remove_cvref_t<Derived>::RowsAtCompileTime;
526 static_assert(colsFull == rowsFull, "staticCondensation only allowed for square matrices");
527 std::array<size_t, rowsFull - sizeOfCondensedIndices> remainingIndices{};
528 std::ranges::set_difference(std::ranges::iota_view(size_t(0), size_t(colsFull)), indices, remainingIndices.begin());
529
530 auto K11 = E(remainingIndices, remainingIndices);
531 auto K12 = E(indices, remainingIndices);
532 auto K22 = E(indices, indices);
533
534 return (K11 - K12.transpose() * K22.inverse() * K12).eval();
535}
536
537template <typename Derived, size_t sizeOfCondensedIndices>
538auto reduceMatrix(const Eigen::MatrixBase<Derived>& E, const std::array<size_t, sizeOfCondensedIndices>& indices) {
539 constexpr size_t colsFull = std::remove_cvref_t<Derived>::ColsAtCompileTime;
540 constexpr size_t rowsFull = std::remove_cvref_t<Derived>::RowsAtCompileTime;
541 static_assert(colsFull == rowsFull, "reduceMatrix only allowed for square matrices");
542 std::array<size_t, rowsFull - sizeOfCondensedIndices> remainingIndices{};
543 std::ranges::set_difference(std::ranges::iota_view(size_t(0), size_t(colsFull)), indices, remainingIndices.begin());
544
545 auto K11 = E(remainingIndices, remainingIndices);
546
547 return K11.eval();
548}
549
562template <typename Derived, size_t sizeOfRemovedCols>
563auto removeCol(const Eigen::MatrixBase<Derived>& E, const std::array<size_t, sizeOfRemovedCols>& indices) {
564 constexpr size_t colsFull = std::remove_cvref_t<Derived>::ColsAtCompileTime;
565 constexpr size_t rowsFull = std::remove_cvref_t<Derived>::RowsAtCompileTime;
566 static_assert(colsFull == 1);
567
568 std::array<size_t, rowsFull - sizeOfRemovedCols> remainingIndices{};
569 std::ranges::set_difference(std::ranges::iota_view(size_t(0), size_t(rowsFull)), indices, remainingIndices.begin());
570
571 return (E(remainingIndices)).eval();
572}
573
588template <typename ST, typename MaterialImpl>
589auto toVoigtAndMaybeReduce(const Eigen::Matrix<ST, 3, 3>& E, [[maybe_unused]] const MaterialImpl& material,
590 bool isStrain = true) {
591 if constexpr (!MaterialImpl::isReduced)
592 return toVoigt(E, isStrain);
593 else {
594 auto ev = toVoigt(E, isStrain);
595 static_assert(decltype(ev)::RowsAtCompileTime == 6 and decltype(ev)::ColsAtCompileTime == 1);
596
597 auto evRed = removeCol(ev, MaterialImpl::fixedVoigtIndices);
598 static_assert(decltype(evRed)::RowsAtCompileTime == 6 - MaterialImpl::fixedVoigtIndices.size() and
599 decltype(evRed)::ColsAtCompileTime == 1);
600 return evRed;
601 }
602}
603
616template <typename M, typename Derived>
617decltype(auto) enlargeIfReduced(const Eigen::MatrixBase<Derived>& E) {
618 if constexpr (!Concepts::EigenVector6<Derived> and Concepts::EigenVector<Derived>) {
619 static_assert(M::isReduced, "You should only end up here, if your material is reduced");
620
621 auto freeindices = M::MaterialImpl::freeVoigtIndices;
622 auto p_E = Eigen::Vector<typename M::MaterialImpl::ScalarType, 6>::Zero().eval();
623 for (int ri = 0; auto i : freeindices)
624 p_E(i) = E(ri++);
625 return p_E;
626
627 } else if constexpr (Concepts::EigenMatrix66<Derived> or Concepts::EigenMatrix33<Derived> or
628 Concepts::EigenVector6<Derived>) {
629 return E.derived();
630 } else {
631 static_assert(M::isReduced, "You should only end up here, if your material is reduced");
632
633 auto freeindices = M::MaterialImpl::freeVoigtIndices;
634 auto p_E = Eigen::Matrix<typename M::MaterialImpl::ScalarType, 6, 6>::Zero().eval();
635 for (int ri = 0; auto i : freeindices) {
636 for (int rj = 0; auto j : freeindices)
637 p_E(i, j) = E(ri, rj++);
638 ++ri;
639 }
640 return p_E;
641 }
642}
643
644} // namespace Ikarus
auto operator+(const Eigen::MatrixBase< Derived > &a, const Eigen::DiagonalMatrix< Scalar, size > &b)
Eigen::Matrix + Eigen::DiagonalMatrix addition missing in Eigen.
Definition: linearalgebrahelper.hh:335
Dune::BlockVector< Type > & operator+=(Dune::BlockVector< Type > &a, const Eigen::MatrixBase< Derived > &b)
Enables the += operator for Dune::BlockVector += Eigen::Vector.
Definition: linearalgebrahelper.hh:179
auto norm(const Eigen::MatrixBase< Derived > &v)
Adding free norm function to Eigen types.
Definition: linearalgebrahelper.hh:258
Derived skew(const Eigen::MatrixBase< Derived > &A)
Returns the skew part of a matrix.
Definition: linearalgebrahelper.hh:434
auto staticCondensation(const Eigen::MatrixBase< Derived > &E, const std::array< size_t, sizeOfCondensedIndices > &indices)
Performs static condensation on a square matrix.
Definition: linearalgebrahelper.hh:522
auto removeCol(const Eigen::MatrixBase< Derived > &E, const std::array< size_t, sizeOfRemovedCols > &indices)
Removes specified columns from a matrix.
Definition: linearalgebrahelper.hh:563
decltype(auto) enlargeIfReduced(const Eigen::MatrixBase< Derived > &E)
Enlarges a matrix if it reduced in the context of material laws, i.e., VanishingStress If the materia...
Definition: linearalgebrahelper.hh:617
auto toVoigtAndMaybeReduce(const Eigen::Matrix< ST, 3, 3 > &E, const MaterialImpl &material, bool isStrain=true)
Converts a 3x3 matrix to Voigt notation, possibly reducing it based on material properties.
Definition: linearalgebrahelper.hh:589
void printForMaple(const Eigen::EigenBase< Derived > &A)
Method to print the matrix in a format that can directly be copied to Maple.
Definition: linearalgebrahelper.hh:445
auto viewAsFlatEigenVector(Dune::BlockVector< ValueType > &blockedVector)
View Dune::BlockVector as an Eigen::Vector.
Definition: linearalgebrahelper.hh:56
size_t correctionSize(const Dune::BlockVector< Type > &a)
Returns the total correction size of a block vector with a Manifold as the underlying type.
Definition: linearalgebrahelper.hh:149
auto floatingPointNorm(const Eigen::MatrixBase< Derived > &v)
Adding free floatingPointNorm function to Eigen types this is an indirection since otherwise norm fai...
Definition: linearalgebrahelper.hh:272
auto createRandomVector(typename FieldVectorT::value_type lower=-1, typename FieldVectorT::value_type upper=1)
Creates a random vector of the specified type within a given range.
Definition: linearalgebrahelper.hh:468
constexpr auto voigtNotationContainer
Container for Voigt notation indices based on dimension.1D: 0,0 2D: 0,0; 1,1; 0,1 3D: 0,...
Definition: linearalgebrahelper.hh:507
Dune::BlockVector< Type > & operator-=(Dune::BlockVector< Type > &a, const Eigen::MatrixBase< Derived > &b)
Enables the -= operator for Dune::BlockVector += Eigen::Vector.
Definition: linearalgebrahelper.hh:199
auto operator-(const Eigen::DiagonalMatrix< Scalar, size > &a)
Unary minus for Eigen::DiagonalMatrix.
Definition: linearalgebrahelper.hh:365
auto orthonormalizeMatrixColumns(const Eigen::MatrixBase< Derived > &A)
Orthonormalizes all Matrix columns using Gram-Schmidt Orthogonalization.
Definition: linearalgebrahelper.hh:34
Derived sym(const Eigen::MatrixBase< Derived > &A)
Returns the symmetric part of a matrix.
Definition: linearalgebrahelper.hh:422
size_t valueSize(const Dune::BlockVector< Type > &a)
Returns the total value size of a block vector with a Manifold as the underlying type.
Definition: linearalgebrahelper.hh:163
auto viewAsEigenMatrixFixedDyn(Dune::BlockVector< ValueType > &blockedVector)
View Dune::BlockVector as an Eigen::Matrix with fixed rows depending on the size of the ValueType and...
Definition: linearalgebrahelper.hh:119
std::ostream & operator<<(std::ostream &os, const Eigen::DiagonalMatrix< Scalar, size > &a)
Output stream operator for Eigen::DiagonalMatrix.
Definition: linearalgebrahelper.hh:409
auto operator*(const Eigen::DiagonalMatrix< Scalar, size > &a, const Eigen::DiagonalMatrix< Scalar, size > &b)
Eigen::DiagonalMatrix Product Missing in Eigen.
Definition: linearalgebrahelper.hh:305
Dune::BlockVector< ManifoldPoint > & addInEmbedding(Dune::BlockVector< ManifoldPoint > &a, const Eigen::MatrixBase< Derived > &b)
Enables the addition in the embedding space of a vector in the space M^n, where M is a manifold with ...
Definition: linearalgebrahelper.hh:238
auto viewAsEigenMatrixAsDynFixed(Dune::BlockVector< ValueType > &blockedVector)
View Dune::BlockVector as an Eigen::Matrix with dynamic rows and fixed columns depending on the size ...
Definition: linearalgebrahelper.hh:87
constexpr Eigen::Index toVoigt(Eigen::Index i, Eigen::Index j) noexcept
Converts 2D indices to Voigt notation index.
Definition: tensorutils.hh:182
Definition: assemblermanipulatorbuildingblocks.hh:22
auto reduceMatrix(const Eigen::MatrixBase< Derived > &E, const std::array< size_t, sizeOfCondensedIndices > &indices)
Definition: linearalgebrahelper.hh:538
Definition: truncatedconjugategradient.hh:24
Definition: utils/concepts.hh:30
Concept defining the requirements for types that support in-place addition.
Definition: utils/concepts.hh:309
Concept defining the requirements for Eigen vectors.
Definition: utils/concepts.hh:365
Concept to check if the underlying scalar type is a dual type.
Definition: utils/concepts.hh:625
Several concepts.