diff --git a/MathLib/KelvinVector.cpp b/MathLib/KelvinVector.cpp
index 94ccd6482bf1dd0cc2cd7105190f89e804aca713..0c6718ddb8f9b2d73b066ed5cf3420ab59d19240 100644
--- a/MathLib/KelvinVector.cpp
+++ b/MathLib/KelvinVector.cpp
@@ -80,6 +80,37 @@ Eigen::Matrix<double, 3, 3> kelvinVectorToTensor(
return m;
}
+template <>
+KelvinVectorType<2> tensorToKelvin<2>(Eigen::Matrix<double, 3, 3> const& m)
+{
+ assert(std::abs(m(0, 1) - m(1, 0)) <
+ std::numeric_limits<double>::epsilon());
+ assert(m(0, 2) == 0);
+ assert(m(1, 2) == 0);
+ assert(m(2, 0) == 0);
+ assert(m(2, 1) == 0);
+
+ KelvinVectorType<2> v;
+ v << m(0, 0), m(1, 1), m(2, 2), m(0, 1) * std::sqrt(2.);
+ return v;
+}
+
+template <>
+KelvinVectorType<3> tensorToKelvin<3>(Eigen::Matrix<double, 3, 3> const& m)
+{
+ assert(std::abs(m(0, 1) - m(1, 0)) <
+ std::numeric_limits<double>::epsilon());
+ assert(std::abs(m(1, 2) - m(2, 1)) <
+ std::numeric_limits<double>::epsilon());
+ assert(std::abs(m(0, 2) - m(2, 0)) <
+ std::numeric_limits<double>::epsilon());
+
+ KelvinVectorType<3> v;
+ v << m(0, 0), m(1, 1), m(2, 2), m(0, 1) * std::sqrt(2.),
+ m(1, 2) * std::sqrt(2.), m(0, 2) * std::sqrt(2.);
+ return v;
+}
+
template <>
Eigen::Matrix<double, 4, 1> kelvinVectorToSymmetricTensor(
Eigen::Matrix<double, 4, 1, Eigen::ColMajor, 4, 1> const& v)
@@ -124,6 +155,5 @@ kelvinVectorToSymmetricTensor(Eigen::Matrix<double,
v.size());
}
}
-
} // namespace KelvinVector
} // namespace MathLib
diff --git a/MathLib/KelvinVector.h b/MathLib/KelvinVector.h
index d57a1eac3b84d9ec53af4bf6106616fa20a9d8fb..bf77b34534b7311b5461e92b36c458941bfe7b01 100644
--- a/MathLib/KelvinVector.h
+++ b/MathLib/KelvinVector.h
@@ -9,6 +9,7 @@
#pragma once
#include <Eigen/Dense>
+#include "BaseLib/Error.h"
namespace MathLib
{
@@ -125,6 +126,12 @@ Eigen::Matrix<double, 3, 3> kelvinVectorToTensor(Eigen::Matrix<double,
KelvinVectorSize,
1> const& v);
+/// Conversion of a 3x3 matrix to a Kelvin vector.
+/// Only implementations for KelvinVectorSize 4 and 6 are provided.
+template <int DisplacementDim>
+KelvinVectorType<DisplacementDim> tensorToKelvin(
+ Eigen::Matrix<double, 3, 3> const& m);
+
/// Conversion of a Kelvin vector to a short vector representation of a
/// symmetric 3x3 matrix.
///
@@ -132,6 +139,8 @@ Eigen::Matrix<double, 3, 3> kelvinVectorToTensor(Eigen::Matrix<double,
/// In the 3D case the entries for the xx, yy, zz, xy, yz, and xz components in
/// that particular order are stored.
///
+/// This is opposite of the symmetricTensorToKelvinVector()
+///
/// Only implementations for KelvinVectorSize 4 and 6, and dynamic size vectors
/// are provided.
template <int KelvinVectorSize>
@@ -143,6 +152,50 @@ kelvinVectorToSymmetricTensor(Eigen::Matrix<double,
KelvinVectorSize,
1> const& v);
+/// Conversion of a short vector representation of a
+/// symmetric 3x3 matrix to a Kelvin vector.
+///
+/// This is opposite of the kelvinVectorToSymmetricTensor()
+///
+/// Only implementations for KelvinVectorSize 4 and 6, and dynamic size vectors
+/// are provided.
+template <typename Derived>
+Eigen::Matrix<double, Eigen::MatrixBase<Derived>::RowsAtCompileTime, 1>
+symmetricTensorToKelvinVector(Eigen::MatrixBase<Derived> const& v)
+{
+ static_assert(
+ (Eigen::MatrixBase<Derived>::ColsAtCompileTime == 1) ||
+ (Eigen::MatrixBase<Derived>::ColsAtCompileTime == Eigen::Dynamic),
+ "KelvinVector must be a column vector");
+ if (v.cols() != 1)
+ {
+ OGS_FATAL(
+ "KelvinVector must be a column vector, but input has %d columns.",
+ v.cols());
+ }
+
+ Eigen::Matrix<double, Eigen::MatrixBase<Derived>::RowsAtCompileTime, 1>
+ result;
+ if (v.rows() == 4)
+ {
+ result.resize(4, 1);
+ result << v[0], v[1], v[2], v[3] * std::sqrt(2.);
+ }
+ else if (v.rows() == 6)
+ {
+ result.resize(6, 1);
+ result << v[0], v[1], v[2], v[3] * std::sqrt(2.), v[4] * std::sqrt(2.),
+ v[5] * std::sqrt(2.);
+ }
+ else
+ {
+ OGS_FATAL(
+ "Symmetric tensor to Kelvin vector conversion expected an input "
+ "vector of size 4 or 6, but a vector of size %d was given.",
+ v.size());
+ }
+ return result;
+}
} // namespace KelvinVector
} // namespace MathLib
diff --git a/Tests/MaterialLib/KelvinVector.cpp b/Tests/MathLib/KelvinVector.cpp
similarity index 73%
rename from Tests/MaterialLib/KelvinVector.cpp
rename to Tests/MathLib/KelvinVector.cpp
index b7be7b8115346045c651e153dc9dd9feea8475bc..2c84be49035ef68630e20f941a5b4ea460bea9d7 100644
--- a/Tests/MaterialLib/KelvinVector.cpp
+++ b/Tests/MathLib/KelvinVector.cpp
@@ -15,39 +15,6 @@
using namespace MathLib::KelvinVector;
namespace ac = autocheck;
-template <int Size>
-using KelvinVector = Eigen::Matrix<double, Size, 1, Eigen::ColMajor, Size, 1>;
-
-template <int Size>
-KelvinVector<Size> tensorToKelvin(Eigen::Matrix<double, 3, 3> const& m);
-
-template <>
-KelvinVector<4> tensorToKelvin(Eigen::Matrix<double, 3, 3> const& m)
-{
- EXPECT_NEAR(m(0, 1), m(1, 0), std::numeric_limits<double>::epsilon());
- EXPECT_EQ(m(0, 2), 0);
- EXPECT_EQ(m(1, 2), 0);
- EXPECT_EQ(m(2, 0), 0);
- EXPECT_EQ(m(2, 1), 0);
-
- KelvinVector<4> v;
- v << m(0, 0), m(1, 1), m(2, 2), m(0, 1) * std::sqrt(2.);
- return v;
-}
-
-template <>
-KelvinVector<6> tensorToKelvin(Eigen::Matrix<double, 3, 3> const& m)
-{
- EXPECT_NEAR(m(0, 1), m(1, 0), std::numeric_limits<double>::epsilon());
- EXPECT_NEAR(m(1, 2), m(2, 1), std::numeric_limits<double>::epsilon());
- EXPECT_NEAR(m(0, 2), m(2, 0), std::numeric_limits<double>::epsilon());
-
- KelvinVector<6> v;
- v << m(0, 0), m(1, 1), m(2, 2), m(0, 1) * std::sqrt(2.),
- m(1, 2) * std::sqrt(2.), m(0, 2) * std::sqrt(2.);
- return v;
-}
-
struct MaterialLibSolidsKelvinVector4 : public ::testing::Test
{
static const int size = 4;
@@ -71,23 +38,23 @@ struct MaterialLibSolidsKelvinVector6 : public ::testing::Test
TEST_F(MaterialLibSolidsKelvinVector4, SelfTestMappingKelvinToTensor)
{
- auto f = [](KelvinVector<4> const& v) {
- return (v - tensorToKelvin<4>(kelvinVectorToTensor(v))).norm() <=
+ auto f = [](KelvinVectorType<2> const& v) {
+ return (v - tensorToKelvin<2>(kelvinVectorToTensor(v))).norm() <=
2 * std::numeric_limits<double>::epsilon() * v.norm();
};
- ac::check<KelvinVector<4>>(
+ ac::check<KelvinVectorType<2>>(
f, 1000, ac::make_arbitrary(kelvinVectorGenerator), gtest_reporter);
}
TEST_F(MaterialLibSolidsKelvinVector6, SelfTestMappingKelvinToTensor)
{
- auto f = [](KelvinVector<6> const& v) {
- return (v - tensorToKelvin<6>(kelvinVectorToTensor(v))).norm() <=
+ auto f = [](KelvinVectorType<3> const& v) {
+ return (v - tensorToKelvin<3>(kelvinVectorToTensor(v))).norm() <=
1.5 * std::numeric_limits<double>::epsilon() * v.norm();
};
- ac::check<KelvinVector<6>>(
+ ac::check<KelvinVectorType<3>>(
f, 1000, ac::make_arbitrary(kelvinVectorGenerator), gtest_reporter);
}
@@ -97,27 +64,27 @@ TEST_F(MaterialLibSolidsKelvinVector6, SelfTestMappingKelvinToTensor)
TEST_F(MaterialLibSolidsKelvinVector4, Determinant)
{
- auto f = [](KelvinVector<4> const& v) {
+ auto f = [](KelvinVectorType<2> const& v) {
return std::abs(Invariants<4>::determinant(v) -
kelvinVectorToTensor(v).determinant()) <=
std::numeric_limits<double>::epsilon() *
std::pow(v.norm(), 3.07);
};
- ac::check<KelvinVector<4>>(
+ ac::check<KelvinVectorType<2>>(
f, 1000, ac::make_arbitrary(kelvinVectorGenerator), gtest_reporter);
}
TEST_F(MaterialLibSolidsKelvinVector6, Determinant)
{
- auto f = [](KelvinVector<6> const& v) {
+ auto f = [](KelvinVectorType<3> const& v) {
return std::abs(Invariants<6>::determinant(v) -
kelvinVectorToTensor(v).determinant()) <=
std::numeric_limits<double>::epsilon() *
std::pow(v.norm(), 3.07);
};
- ac::check<KelvinVector<6>>(
+ ac::check<KelvinVectorType<3>>(
f, 1000, ac::make_arbitrary(kelvinVectorGenerator), gtest_reporter);
}
@@ -127,18 +94,18 @@ TEST_F(MaterialLibSolidsKelvinVector6, Determinant)
TEST_F(MaterialLibSolidsKelvinVector4, Inverse)
{
- auto f = [](KelvinVector<4> const& v) {
+ auto f = [](KelvinVectorType<2> const& v) {
auto const error =
- (inverse(v) - tensorToKelvin<4>(kelvinVectorToTensor(v).inverse()))
+ (inverse(v) - tensorToKelvin<2>(kelvinVectorToTensor(v).inverse()))
.norm();
// The error is only weekly depending on the input vector norm.
return error < 1e-6 && error < 1e-8 * std::pow(v.norm(), 1.4);
};
- ac::check<KelvinVector<4>>(
+ ac::check<KelvinVectorType<2>>(
f, 1000,
ac::make_arbitrary(kelvinVectorGenerator)
- .discard_if([](KelvinVector<4> const& v) {
+ .discard_if([](KelvinVectorType<2> const& v) {
// only invertable matrices
return (std::abs(kelvinVectorToTensor(v).determinant()) == 0);
}),
@@ -147,18 +114,18 @@ TEST_F(MaterialLibSolidsKelvinVector4, Inverse)
TEST_F(MaterialLibSolidsKelvinVector6, Inverse)
{
- auto f = [](KelvinVector<6> const& v) {
+ auto f = [](KelvinVectorType<3> const& v) {
auto const error =
- (inverse(v) - tensorToKelvin<6>(kelvinVectorToTensor(v).inverse()))
+ (inverse(v) - tensorToKelvin<3>(kelvinVectorToTensor(v).inverse()))
.norm();
// The error is only weekly depending on the input vector norm.
return error < 1e-6 && error < 1e-8 * std::pow(v.norm(), 1.4);
};
- ac::check<KelvinVector<6>>(
+ ac::check<KelvinVectorType<3>>(
f, 1000,
ac::make_arbitrary(kelvinVectorGenerator)
- .discard_if([](KelvinVector<6> const& v) {
+ .discard_if([](KelvinVectorType<3> const& v) {
// only invertable matrices
return (std::abs(kelvinVectorToTensor(v).determinant()) == 0);
}),