diff --git a/MathLib/EigenBlockMatrixView.h b/MathLib/EigenBlockMatrixView.h new file mode 100644 index 0000000000000000000000000000000000000000..55ee6a4867db7cef936afe8b1316295bb4735893 --- /dev/null +++ b/MathLib/EigenBlockMatrixView.h @@ -0,0 +1,56 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#pragma once + +#include <Eigen/Core> + +namespace MathLib +{ +template <int D, typename M> +class EigenBlockMatrixViewFunctor +{ +public: + static constexpr int rows = M::RowsAtCompileTime == Eigen::Dynamic + ? Eigen::Dynamic + : D * M::RowsAtCompileTime; + static constexpr int cols = M::ColsAtCompileTime == Eigen::Dynamic + ? Eigen::Dynamic + : D * M::ColsAtCompileTime; + using Scalar = typename M::Scalar; + using Matrix = Eigen::Matrix<Scalar, rows, cols, Eigen::ColMajor>; + + constexpr explicit EigenBlockMatrixViewFunctor(const M& matrix) + : matrix_(matrix){}; + + constexpr const Scalar operator()(Eigen::Index row, Eigen::Index col) const + { + if (row / matrix_.rows() != col / matrix_.cols()) + { + return 0; + } + return matrix_(row % matrix_.rows(), col % matrix_.cols()); + } + +private: + const typename M::Nested& matrix_; +}; + +template <int D, typename M> +constexpr Eigen::CwiseNullaryOp< + EigenBlockMatrixViewFunctor<D, M>, + typename EigenBlockMatrixViewFunctor<D, M>::Matrix> +eigenBlockMatrixView(const Eigen::MatrixBase<M>& matrix) +{ + using Matrix = typename EigenBlockMatrixViewFunctor<D, M>::Matrix; + return Matrix::NullaryExpr( + D * matrix.rows(), D * matrix.cols(), + EigenBlockMatrixViewFunctor<D, M>(matrix.derived())); +} +} // namespace MathLib diff --git a/Tests/MathLib/EigenBlockMatrixView.cpp b/Tests/MathLib/EigenBlockMatrixView.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5bd37b43d3d9e88400251a4a300bb80cba287f9a --- /dev/null +++ b/Tests/MathLib/EigenBlockMatrixView.cpp @@ -0,0 +1,90 @@ +/** + * \file + * \copyright + * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + */ + +#include "MathLib/EigenBlockMatrixView.h" + +#include <gtest/gtest.h> + +#include <Eigen/Core> + +TEST(MathLib_EigenBlockMatriView, 1by1) +{ + Eigen::Matrix<double, 4, 4> expected = + Eigen::Vector4d::Constant(42).asDiagonal(); + + // fixed size + Eigen::Matrix<double, 1, 1> A{42}; + ASSERT_EQ(expected, MathLib::eigenBlockMatrixView<4>(A)); + + // dynamic size + Eigen::MatrixXd B(1, 1); + B << 42; + ASSERT_TRUE(expected == MathLib::eigenBlockMatrixView<4>(B)); +} + +TEST(MathLib_EigenBlockMatriView, 1byN) +{ + Eigen::Matrix<double, 6, 3> expected; + // clang-format off + expected << 1, 0, 0, + 2, 0, 0, + 0, 1, 0, + 0, 2, 0, + 0, 0, 1, + 0, 0, 2; + // clang-format on + + Eigen::Vector2d v{1, 2}; + ASSERT_TRUE(expected == MathLib::eigenBlockMatrixView<3>(v)); + + // dynamic + Eigen::MatrixXd w(2, 1); + w << 1, 2; + ASSERT_TRUE(expected == MathLib::eigenBlockMatrixView<3>(w)); +} + +TEST(MathLib_EigenBlockMatriView, Nby1) +{ + Eigen::Matrix<double, 3, 6> expected; + // clang-format off + expected << 1, 2, 0, 0, 0, 0, + 0, 0, 1, 2, 0, 0, + 0, 0, 0, 0, 1, 2; + // clang-format on + + Eigen::RowVector2d v{1, 2}; + ASSERT_TRUE(expected == MathLib::eigenBlockMatrixView<3>(v)); + + // dynamic + Eigen::MatrixXd w(1, 2); + w << 1, 2; + ASSERT_TRUE(expected == MathLib::eigenBlockMatrixView<3>(w)); +} + +TEST(MathLib_EigenBlockMatriView, NbyM) +{ + Eigen::Matrix<double, 6, 6> expected; + // clang-format off + expected << 1, 2, 0, 0, 0, 0, + 3, 4, 0, 0, 0, 0, + 0, 0, 1, 2, 0, 0, + 0, 0, 3, 4, 0, 0, + 0, 0, 0, 0, 1, 2, + 0, 0, 0, 0, 3, 4; + // clang-format on + + Eigen::Matrix2d A; + A << 1, 2, 3, 4; + ASSERT_TRUE(expected == MathLib::eigenBlockMatrixView<3>(A)); + + // dynamic + Eigen::MatrixXd B(2, 2); + B << 1, 2, 3, 4; + ASSERT_TRUE(expected == MathLib::eigenBlockMatrixView<3>(B)); +}