From 684eab7b8a2fd6cf3d6dc55a6904aeff678a7d3a Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Thu, 28 Dec 2023 14:06:58 +0100
Subject: [PATCH] [ML] Vectorized tensor utilities
The deformation gradient and alike quantities are second order tensors
stored as a vector. Convenience operations on these tensors are
collected in the namespace.
The conversion to a Eigen matrix and back are used so far in tests only
for comparisons with Eigen operations.
---
MathLib/VectorizedTensor.cpp | 28 ++++
MathLib/VectorizedTensor.h | 229 +++++++++++++++++++++++++++++
Tests/MathLib/VectorizedTensor.cpp | 225 ++++++++++++++++++++++++++++
3 files changed, 482 insertions(+)
create mode 100644 MathLib/VectorizedTensor.cpp
create mode 100644 MathLib/VectorizedTensor.h
create mode 100644 Tests/MathLib/VectorizedTensor.cpp
diff --git a/MathLib/VectorizedTensor.cpp b/MathLib/VectorizedTensor.cpp
new file mode 100644
index 00000000000..c760468bee7
--- /dev/null
+++ b/MathLib/VectorizedTensor.cpp
@@ -0,0 +1,28 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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 "VectorizedTensor.h"
+
+namespace MathLib ::VectorizedTensor
+{
+
+bool isTensorConvertibleTo1d(Eigen::Matrix3d const& tensor)
+{
+ return (tensor(0, 1) == 0. && tensor(1, 0) == 0. && //
+ tensor(0, 2) == 0. && tensor(2, 0) == 0. && //
+ tensor(1, 2) == 0. && tensor(2, 1) == 0.);
+}
+
+bool isTensorConvertibleTo2d(Eigen::Matrix3d const& tensor)
+{
+ return (tensor(0, 2) == 0. && tensor(1, 2) == 0. && //
+ tensor(2, 0) == 0. && tensor(2, 1) == 0.);
+}
+
+} // namespace MathLib::VectorizedTensor
diff --git a/MathLib/VectorizedTensor.h b/MathLib/VectorizedTensor.h
new file mode 100644
index 00000000000..3a70bee3b0d
--- /dev/null
+++ b/MathLib/VectorizedTensor.h
@@ -0,0 +1,229 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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>
+#include <Eigen/LU>
+#include <type_traits>
+
+#include "BaseLib/Error.h"
+
+namespace MathLib
+{
+/// Namespace for vectorized second-order tensor types and operations on those.
+/// The tensors are for example, deformation gradients.
+namespace VectorizedTensor
+{
+/// Vectorized tensor size for given displacement dimension.
+constexpr int size(int const displacement_dim)
+{
+ if (displacement_dim == 1)
+ {
+ return 3;
+ }
+ if (displacement_dim == 2)
+ {
+ return 5;
+ }
+ if (displacement_dim == 3)
+ {
+ return 9;
+ }
+ OGS_FATAL(
+ "Cannot convert displacement dimension {} to vectorized tensor size.",
+ displacement_dim);
+}
+
+/// Displacement dimension of a vectorized tensor.
+template <typename VectorizedTensor>
+constexpr int dimension()
+{
+ static_assert(VectorizedTensor::ColsAtCompileTime == 1);
+ constexpr int rows = VectorizedTensor::RowsAtCompileTime;
+ if (rows == 3)
+ {
+ return 1;
+ }
+ if (rows == 5)
+ {
+ return 2;
+ }
+ if (rows == 9)
+ {
+ return 3;
+ }
+ OGS_FATAL(
+ "Cannot convert vectorized tensor of size {} to displacement "
+ "dimension.",
+ rows);
+}
+
+/// Vectorized tensor type for given displacement dimension.
+/// \note The Eigen vector is always a fixed size vector in contrast to a shape
+/// matrix policy types like GMatrixPolicyType::GradientVectorType.
+template <int DisplacementDim>
+using Type = Eigen::Matrix<double, size(DisplacementDim), 1, Eigen::ColMajor>;
+
+/// Vectorized identity tensor expression. Corresponds to 3x3 identity matrix in
+/// all dimensions.
+template <int DisplacementDim>
+constexpr auto identity()
+{
+ static_assert(
+ 0 < DisplacementDim && DisplacementDim <= 3,
+ "Identity is implemented only for displacement dimension 1, 2, or 3.");
+ return Type<DisplacementDim>::NullaryExpr(
+ size(DisplacementDim), 1,
+ [](Eigen::Index const row, [[maybe_unused]] Eigen::Index const col)
+ {
+ assert(col == 0);
+ if constexpr (DisplacementDim == 1)
+ { // All entries are diagonal entries.
+ return 1;
+ }
+ if constexpr (DisplacementDim == 2)
+ {
+ if (row == 0 || row == 3 || row == 4)
+ {
+ return 1;
+ }
+ }
+ if constexpr (DisplacementDim == 3)
+ {
+ if (row == 0 || row == 4 || row == 8)
+ {
+ return 1;
+ }
+ }
+ return 0;
+ });
+}
+
+/// Computes determinant of a vectorized tensor.
+template <typename Derived>
+double determinant(Eigen::MatrixBase<Derived> const& tensor)
+{
+ constexpr int displacement_dim = dimension<Derived>();
+ static_assert(0 < displacement_dim && displacement_dim <= 3,
+ "Vectorized tensor determinant is implemented only for "
+ "displacement dimension 1, 2, or 3.");
+
+ if constexpr (displacement_dim == 1)
+ {
+ return tensor[0] * tensor[1] * tensor[2];
+ }
+ if constexpr (displacement_dim == 2)
+ {
+ Eigen::Map<Eigen::Matrix2d const> const top_left{
+ tensor.derived().data()};
+
+ return top_left.determinant() * tensor[4];
+ }
+ if constexpr (displacement_dim == 3)
+ {
+ return Eigen::Map<Eigen::Matrix3d const>(tensor.derived().data())
+ .determinant();
+ }
+}
+
+/// Only a diagonal tensor can be converted to a 1d vectorized tensor.
+bool isTensorConvertibleTo1d(Eigen::Matrix3d const& tensor);
+
+/// Only a tensor of form
+/// \f$ \left( \begin{array}{ccc}
+/// a & b & 0 \\%
+/// c & d & 0 \\%
+/// 0 & 0 & e \\%
+/// \end{array} \right)\f$
+/// can be converted to a 2d vectorized tensor.
+bool isTensorConvertibleTo2d(Eigen::Matrix3d const& tensor);
+
+/// Converts a 3x3 matrix expression to a vectorized tensor if the conversion
+/// for that dimension is possible; see isTensorConvertibleTo1d() and
+/// isTensorConvertibleTo2d() functions.
+template <int DisplacementDim, typename Derived>
+Type<DisplacementDim> toVector(Eigen::MatrixBase<Derived> const& tensor)
+{
+ static_assert(0 < DisplacementDim && DisplacementDim <= 3,
+ "Conversion to displacement dimension other than 1, 2, or 3 "
+ "is not valid.");
+
+ constexpr int rows = Derived::RowsAtCompileTime;
+ constexpr int cols = Derived::ColsAtCompileTime;
+ static_assert(rows == 3 || rows == Eigen::Dynamic);
+ static_assert(cols == 3 || cols == Eigen::Dynamic);
+ if (tensor.rows() != 3 || tensor.cols() != 3)
+ {
+ OGS_FATAL(
+ "Incorrect tensor size, must be 3x3, but tensor is {:d}x{:d}.",
+ tensor.rows(), tensor.cols());
+ }
+
+ if constexpr (DisplacementDim == 1)
+ {
+ if (!isTensorConvertibleTo1d(tensor))
+ {
+ OGS_FATAL(
+ "Cannot convert a tensor with non-zero off-diagonal elements "
+ "to a 1d vectorized tensor representation.");
+ }
+ return tensor.diagonal();
+ }
+ if constexpr (DisplacementDim == 2)
+ {
+ if (!isTensorConvertibleTo2d(tensor))
+ {
+ OGS_FATAL(
+ "Cannot convert a tensor with non-zero elements at (0, 2), (1, "
+ "2), (2, 0), and (2, 1) positions to a 2d vectorized tensor "
+ "representation.");
+ }
+ Type<2> result;
+ result.template head<4>() =
+ tensor.template block<2, 2>(0, 0).reshaped();
+ result(4) = tensor(2, 2);
+ return result;
+ }
+ if constexpr (DisplacementDim == 3)
+ {
+ return tensor.reshaped();
+ }
+}
+
+/// Converts a vectorized tensor to a 3x3 matrix.
+template <int DisplacementDim>
+Eigen::Matrix3d toTensor(Type<DisplacementDim> const& tensor)
+{
+ static_assert(
+ DisplacementDim == 1 || DisplacementDim == 2 || DisplacementDim == 3,
+ "Conversion from displacement dimension other than 1, 2, or 3 "
+ "is not valid.");
+
+ using Matrix = Eigen::Matrix3d;
+ if constexpr (DisplacementDim == 1)
+ {
+ return tensor.asDiagonal();
+ }
+ if constexpr (DisplacementDim == 2)
+ {
+ Matrix m = Matrix::Zero();
+ m.template block<2, 2>(0, 0) =
+ Eigen::Map<Eigen::Matrix<double, 2, 2> const>(tensor.data());
+ m(2, 2) = tensor(4);
+ return m;
+ }
+ if constexpr (DisplacementDim == 3)
+ {
+ return Eigen::Map<Matrix const>(tensor.data());
+ }
+}
+
+} // namespace VectorizedTensor
+} // namespace MathLib
diff --git a/Tests/MathLib/VectorizedTensor.cpp b/Tests/MathLib/VectorizedTensor.cpp
new file mode 100644
index 00000000000..76822209a38
--- /dev/null
+++ b/Tests/MathLib/VectorizedTensor.cpp
@@ -0,0 +1,225 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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/VectorizedTensor.h"
+
+#include <gtest/gtest.h>
+
+#include <Eigen/LU>
+#include <range/v3/algorithm/for_each.hpp>
+#include <range/v3/view/cartesian_product.hpp>
+#include <range/v3/view/iota.hpp>
+
+#include "Tests/AutoCheckTools.h"
+
+using namespace MathLib;
+namespace ac = autocheck;
+
+TEST(VectorizedTensorTest, DimensionsAndSizes)
+{
+ namespace VT = MathLib::VectorizedTensor;
+
+ // Type implicitly uses size(), which is being tested.
+ static_assert(VT::dimension<VT::Type<1>>() == 1);
+ static_assert(VT::dimension<VT::Type<2>>() == 2);
+ static_assert(VT::dimension<VT::Type<3>>() == 3);
+}
+
+//
+// Self-tests of a tensor in matrix form to a vector form and back conversion.
+//
+template <int DisplacementDim>
+auto vectorizedTensorToMatrixMapping()
+{
+ static constexpr int size = VectorizedTensor::size(DisplacementDim);
+ ac::randomEigenMatrixGenerator<double, size, 1> generator;
+
+ auto to_tensor_and_back =
+ [](VectorizedTensor::Type<DisplacementDim> const& v)
+ {
+ return (v == VectorizedTensor::toVector<DisplacementDim>(
+ VectorizedTensor::toTensor<DisplacementDim>(v)));
+ };
+
+ return ac::check<VectorizedTensor::Type<DisplacementDim>>(
+ to_tensor_and_back, 1000, ac::make_arbitrary(generator),
+ ac::gtest_reporter());
+}
+
+template <int DisplacementDim, typename Generator>
+auto matrixToVectorizedTensorMapping(Generator const& generator)
+{
+ auto to_vector_and_back = [](Eigen::Matrix3d const& t)
+ {
+ return (t == VectorizedTensor::toTensor<DisplacementDim>(
+ VectorizedTensor::toVector<DisplacementDim>(t)));
+ };
+
+ return ac::check<Eigen::Matrix3d>(to_vector_and_back, 1000,
+ ac::make_arbitrary(generator),
+ ac::gtest_reporter());
+}
+
+struct random1dConvertibleEigenMatrixGenerator
+{
+ ac::generator<double> source;
+ using result_type = Eigen::Matrix<double, 3, 3>;
+
+ result_type operator()(std::size_t size = 0)
+ {
+ result_type rv = result_type::Zero();
+ rv.diagonal()(0) = fix(size, source)();
+ rv.diagonal()(1) = fix(size, source)();
+ rv.diagonal()(2) = fix(size, source)();
+ if (!VectorizedTensor::isTensorConvertibleTo1d(rv))
+ {
+ OGS_FATAL(
+ "VectorizedTensorTest; 1d test matrix generation failed.");
+ }
+ return rv;
+ }
+};
+
+struct random2dConvertibleEigenMatrixGenerator
+{
+ ac::generator<double> source;
+ using result_type = Eigen::Matrix<double, 3, 3>;
+
+ result_type operator()(std::size_t size = 0)
+ {
+ result_type rv = random1dConvertibleEigenMatrixGenerator{}(size);
+ rv(0, 1) = fix(size, source)();
+ rv(1, 0) = fix(size, source)();
+ if (!VectorizedTensor::isTensorConvertibleTo2d(rv))
+ {
+ OGS_FATAL(
+ "VectorizedTensorTest; 2d test matrix generation failed.");
+ }
+ return rv;
+ }
+};
+
+TEST(VectorizedTensorTest, SelfTestMapping1)
+{
+ vectorizedTensorToMatrixMapping<1>();
+ matrixToVectorizedTensorMapping<1>(
+ random1dConvertibleEigenMatrixGenerator{});
+}
+TEST(VectorizedTensorTest, SelfTestMapping2)
+{
+ vectorizedTensorToMatrixMapping<2>();
+ matrixToVectorizedTensorMapping<2>(
+ random2dConvertibleEigenMatrixGenerator{});
+}
+TEST(VectorizedTensorTest, SelfTestMapping3)
+{
+ vectorizedTensorToMatrixMapping<3>();
+ matrixToVectorizedTensorMapping<3>(
+ ac::randomEigenMatrixGenerator<double, 3, 3>{});
+}
+
+TEST(VectorizedTensorTest, DynamicSizeTensorToVector)
+{
+ // Testing only for dimension 1, because that code part is dimension
+ // independent.
+ EXPECT_NO_THROW(VectorizedTensor::toVector<1>(Eigen::MatrixXd::Zero(3, 3)));
+ EXPECT_THROW(VectorizedTensor::toVector<1>(Eigen::MatrixXd::Zero(3, 4)),
+ std::runtime_error);
+ EXPECT_THROW(VectorizedTensor::toVector<1>(Eigen::MatrixXd::Zero(4, 3)),
+ std::runtime_error);
+ EXPECT_THROW(VectorizedTensor::toVector<1>(Eigen::MatrixXd::Zero(3, 2)),
+ std::runtime_error);
+ EXPECT_THROW(VectorizedTensor::toVector<1>(Eigen::MatrixXd::Zero(2, 3)),
+ std::runtime_error);
+}
+
+TEST(VectorizedTensorTest, ConvertibleAndNonconvertibleTensors)
+{
+ auto one_ij = [](int i, int j)
+ {
+ Eigen::Matrix3d m = Eigen::Matrix3d::Zero();
+ m(i, j) = 1.;
+ return m;
+ };
+
+ // 1d
+ EXPECT_TRUE(
+ VectorizedTensor::isTensorConvertibleTo1d(Eigen::MatrixXd::Zero(3, 3)));
+ EXPECT_TRUE(VectorizedTensor::isTensorConvertibleTo1d(
+ Eigen::MatrixXd::Identity(3, 3)));
+ using namespace ranges::views;
+ ranges::for_each(
+ cartesian_product(iota(0, 3), iota(0, 3)),
+ [&](auto const ij)
+ {
+ auto const [i, j] = ij;
+ EXPECT_EQ(i == j,
+ VectorizedTensor::isTensorConvertibleTo1d(one_ij(i, j)));
+ });
+
+ // 2d
+ EXPECT_TRUE(
+ VectorizedTensor::isTensorConvertibleTo2d(Eigen::MatrixXd::Zero(3, 3)));
+ EXPECT_TRUE(VectorizedTensor::isTensorConvertibleTo2d(
+ Eigen::MatrixXd::Identity(3, 3)));
+ {
+ EXPECT_FALSE(VectorizedTensor::isTensorConvertibleTo2d(one_ij(0, 2)));
+ EXPECT_FALSE(VectorizedTensor::isTensorConvertibleTo2d(one_ij(1, 2)));
+ EXPECT_FALSE(VectorizedTensor::isTensorConvertibleTo2d(one_ij(2, 0)));
+ EXPECT_FALSE(VectorizedTensor::isTensorConvertibleTo2d(one_ij(2, 1)));
+ }
+}
+//
+// Identities of a vectorized tensors correspond to Eigen matrix identities.
+//
+TEST(VectorizedTensorTest, Identity)
+{
+ EXPECT_EQ(VectorizedTensor::toTensor<1>(VectorizedTensor::identity<1>()),
+ Eigen::Matrix3d::Identity());
+ EXPECT_EQ(VectorizedTensor::toTensor<2>(VectorizedTensor::identity<2>()),
+ Eigen::Matrix3d::Identity());
+ EXPECT_EQ(VectorizedTensor::toTensor<3>(VectorizedTensor::identity<3>()),
+ Eigen::Matrix3d::Identity());
+}
+
+//
+// Determinants of a vectorized tensor and corresponding Eigen::Matrix are
+// equal.
+//
+template <int DisplacementDim>
+auto vectorizedTensorDeterminant()
+{
+ static constexpr int size = VectorizedTensor::size(DisplacementDim);
+ ac::randomEigenMatrixGenerator<double, size, 1> generator;
+
+ auto f = [](VectorizedTensor::Type<DisplacementDim> const& v)
+ {
+ return std::abs(VectorizedTensor::determinant(v) -
+ VectorizedTensor::toTensor<DisplacementDim>(v)
+ .determinant()) <=
+ std::numeric_limits<double>::epsilon() *
+ std::pow(v.norm(), 3.07);
+ };
+
+ ac::check<VectorizedTensor::Type<DisplacementDim>>(
+ f, 1000, ac::make_arbitrary(generator), ac::gtest_reporter());
+}
+
+TEST(VectorizedTensorTest, Determinant1)
+{
+ vectorizedTensorDeterminant<1>();
+}
+TEST(VectorizedTensorTest, Determinant2)
+{
+ vectorizedTensorDeterminant<2>();
+}
+TEST(VectorizedTensorTest, Determinant3)
+{
+ vectorizedTensorDeterminant<3>();
+}
--
GitLab