/**
* \file
* \copyright
* Copyright (c) 2012-2021, 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 <gtest/gtest.h>
#include "MathLib/KelvinVector.h"
#include "Tests/AutoCheckTools.h"
using namespace MathLib::KelvinVector;
namespace ac = autocheck;
struct MaterialLibSolidsKelvinVector4 : public ::testing::Test
{
static const int size = 4;
ac::randomEigenMatrixGenerator<double, size, 1> kelvinVectorGenerator;
ac::gtest_reporter gtest_reporter;
};
struct MaterialLibSolidsKelvinVector6 : public ::testing::Test
{
static const int size = 6;
ac::randomEigenMatrixGenerator<double, size, 1> kelvinVectorGenerator;
ac::gtest_reporter gtest_reporter;
};
//
// Self-test of Kelvin to tensor and back conversion, which should be identity
// up to numerical errors.
//
TEST_F(MaterialLibSolidsKelvinVector4, SelfTestMappingKelvinToTensor)
{
auto f = [](KelvinVectorType<2> const& v) {
return (v - tensorToKelvin<2>(kelvinVectorToTensor(v))).norm() <=
2 * std::numeric_limits<double>::epsilon() * v.norm();
};
ac::check<KelvinVectorType<2>>(
f, 1000, ac::make_arbitrary(kelvinVectorGenerator), gtest_reporter);
}
TEST_F(MaterialLibSolidsKelvinVector6, SelfTestMappingKelvinToTensor)
{
auto f = [](KelvinVectorType<3> const& v) {
return (v - tensorToKelvin<3>(kelvinVectorToTensor(v))).norm() <=
1.5 * std::numeric_limits<double>::epsilon() * v.norm();
};
ac::check<KelvinVectorType<3>>(
f, 1000, ac::make_arbitrary(kelvinVectorGenerator), gtest_reporter);
}
//
// Determinants of a Kelvin vector and corresponding Eigen::Matrix are equal.
//
TEST_F(MaterialLibSolidsKelvinVector4, Determinant)
{
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<KelvinVectorType<2>>(
f, 1000, ac::make_arbitrary(kelvinVectorGenerator), gtest_reporter);
}
TEST_F(MaterialLibSolidsKelvinVector6, Determinant)
{
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<KelvinVectorType<3>>(
f, 1000, ac::make_arbitrary(kelvinVectorGenerator), gtest_reporter);
}
//
// Inverse of a Kelvin vector and corresponding Eigen::Matrix are equal.
//
TEST_F(MaterialLibSolidsKelvinVector4, Inverse)
{
auto f = [](KelvinVectorType<2> const& v) {
auto const error =
(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<KelvinVectorType<2>>(
f, 1000,
ac::make_arbitrary(kelvinVectorGenerator)
.discard_if([](KelvinVectorType<2> const& v) {
// only invertable matrices
return (std::abs(kelvinVectorToTensor(v).determinant()) == 0);
}),
gtest_reporter);
}
TEST_F(MaterialLibSolidsKelvinVector6, Inverse)
{
auto f = [](KelvinVectorType<3> const& v) {
auto const error =
(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<KelvinVectorType<3>>(
f, 1000,
ac::make_arbitrary(kelvinVectorGenerator)
.discard_if([](KelvinVectorType<3> const& v) {
// only invertable matrices
return (std::abs(kelvinVectorToTensor(v).determinant()) == 0);
}),
gtest_reporter);
}
//
// Tests of deviatoric and spherical projection tensors.
//
TEST_F(MaterialLibSolidsKelvinVector4, DeviatoricSphericalProjections)
{
auto const& P_dev = Invariants<size>::deviatoric_projection;
auto const& P_sph = Invariants<size>::spherical_projection;
// Test product P_dev * P_sph is zero.
Eigen::Matrix<double, 4, 4> const P_dev_P_sph = P_dev * P_sph;
EXPECT_LT(P_dev_P_sph.norm(), std::numeric_limits<double>::epsilon());
EXPECT_LT(P_dev_P_sph.maxCoeff(), std::numeric_limits<double>::epsilon());
// Test product P_sph * P_dev is zero.
Eigen::Matrix<double, 4, 4> const P_sph_P_dev = P_sph * P_dev;
EXPECT_LT(P_sph_P_dev.norm(), std::numeric_limits<double>::epsilon());
EXPECT_LT(P_sph_P_dev.maxCoeff(), std::numeric_limits<double>::epsilon());
// Test sum is identity.
EXPECT_EQ(P_dev + P_sph, (Eigen::Matrix<double, size, size>::Identity()));
}
TEST_F(MaterialLibSolidsKelvinVector6, DeviatoricSphericalProjections)
{
auto const& P_dev = Invariants<size>::deviatoric_projection;
auto const& P_sph = Invariants<size>::spherical_projection;
// Test product P_dev * P_sph is zero.
Eigen::Matrix<double, 6, 6> const P_dev_P_sph = P_dev * P_sph;
EXPECT_LT(P_dev_P_sph.norm(), std::numeric_limits<double>::epsilon());
EXPECT_LT(P_dev_P_sph.maxCoeff(), std::numeric_limits<double>::epsilon());
// Test product P_sph * P_dev is zero.
Eigen::Matrix<double, 6, 6> const P_sph_P_dev = P_sph * P_dev;
EXPECT_LT(P_sph_P_dev.norm(), std::numeric_limits<double>::epsilon());
EXPECT_LT(P_sph_P_dev.maxCoeff(), std::numeric_limits<double>::epsilon());
// Test sum is identity.
EXPECT_EQ(P_dev + P_sph, (Eigen::Matrix<double, size, size>::Identity()));
}