From 2ff5c779d3d389fb80dfea343825a6974531effc Mon Sep 17 00:00:00 2001 From: Henrik Jakob <96132706+henrij22@users.noreply.github.com> Date: Mon, 21 Oct 2024 16:47:24 +0200 Subject: [PATCH] Use AutoDiffScalar Concept in materials and fix checkFESByAutoDiff test (#329) --- .../mechanics/enhancedassumedstrains.hh | 13 +++++---- .../mechanics/materials/neohooke.hh | 25 ++++++++++++---- .../mechanics/nonlinearelastic.hh | 29 +++++++------------ tests/src/checkfebyautodiff.hh | 13 +++++---- tests/src/testmaterial.cpp | 24 +++++++++++++++ tests/src/testnonlinearelasticity.hh | 4 +-- tests/src/testnonlinearelasticityneohooke.cpp | 4 +-- 7 files changed, 72 insertions(+), 40 deletions(-) diff --git a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh index 22024a19b..a9ecbe525 100644 --- a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh +++ b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh @@ -18,6 +18,7 @@ #include #include #include + #include namespace Ikarus { @@ -165,10 +166,10 @@ protected: return; decltype(auto) LMat = [this]() -> decltype(auto) { - if constexpr (std::is_same_v) - return [this]() -> Eigen::MatrixXd& { return L_; }(); - else + if constexpr (Concepts::AutodiffScalar) return Eigen::MatrixX{}; + else + return [this]() -> Eigen::MatrixXd& { return L_; }(); }(); auto calculateMatrixContribution = [&](const EAST& easFunction) { @@ -208,10 +209,10 @@ protected: calculateDAndLMatrix(easFunction, par, D, L_); decltype(auto) LMat = [this]() -> decltype(auto) { - if constexpr (std::is_same_v) - return [this]() -> Eigen::MatrixXd& { return L_; }(); - else + if constexpr (Concepts::AutodiffScalar) return Eigen::MatrixX{}; + else + return [this]() -> Eigen::MatrixXd& { return L_; }(); }(); const auto disp = Dune::viewAsFlatEigenVector(uFunction.coefficientsRef()); const auto alpha = (-D.inverse() * L_ * disp).eval(); diff --git a/ikarus/finiteelements/mechanics/materials/neohooke.hh b/ikarus/finiteelements/mechanics/materials/neohooke.hh index 588f5f39e..1dce40833 100644 --- a/ikarus/finiteelements/mechanics/materials/neohooke.hh +++ b/ikarus/finiteelements/mechanics/materials/neohooke.hh @@ -72,11 +72,13 @@ struct NeoHookeT : public Material> * \return ScalarType The stored energy. */ template - ScalarType storedEnergyImpl(const Eigen::MatrixBase& C) const noexcept { + ScalarType storedEnergyImpl(const Eigen::MatrixBase& C) const { static_assert(Concepts::EigenMatrixOrVoigtNotation3); if constexpr (!Concepts::EigenVector) { - const auto traceC = C.trace(); - const auto logdetF = log(sqrt(C.determinant())); + const auto traceC = C.trace(); + const auto detC = C.determinant(); + checkPositiveDetC(detC); + const auto logdetF = log(sqrt(detC)); return materialParameter_.mu / 2.0 * (traceC - 3 - 2 * logdetF) + materialParameter_.lambda / 2.0 * logdetF * logdetF; } else @@ -96,7 +98,9 @@ struct NeoHookeT : public Material> static_assert(Concepts::EigenMatrixOrVoigtNotation3); if constexpr (!voigt) { if constexpr (!Concepts::EigenVector) { - const auto logdetF = log(sqrt(C.determinant())); + const auto detC = C.determinant(); + checkPositiveDetC(detC); + const auto logdetF = log(sqrt(detC)); const auto invC = C.inverse().eval(); return (materialParameter_.mu * (StrainMatrix::Identity() - invC) + materialParameter_.lambda * logdetF * invC) .eval(); @@ -118,8 +122,10 @@ struct NeoHookeT : public Material> auto tangentModuliImpl(const Eigen::MatrixBase& C) const { static_assert(Concepts::EigenMatrixOrVoigtNotation3); if constexpr (!voigt) { - const auto invC = C.inverse().eval(); - const auto logdetF = log(sqrt(C.determinant())); + const auto invC = C.inverse().eval(); + const auto detC = C.determinant(); + checkPositiveDetC(detC); + const auto logdetF = log(sqrt(detC)); const auto CTinv = tensorView(invC, std::array({3, 3})); static_assert(Eigen::TensorFixedSize>::NumIndices == 2); Eigen::TensorFixedSize> moduli = @@ -144,6 +150,13 @@ struct NeoHookeT : public Material> private: MaterialParameters materialParameter_; + + void checkPositiveDetC(ScalarType detC) const { + if (Dune::FloatCmp::le(static_cast(detC), 0.0, 1e-10)) + DUNE_THROW(Dune::InvalidStateException, + "Determinant of right Cauchy Green tensor C must be greater than zero. detC = " + + std::to_string(static_cast(detC))); + } }; /** diff --git a/ikarus/finiteelements/mechanics/nonlinearelastic.hh b/ikarus/finiteelements/mechanics/nonlinearelastic.hh index 863fc6c48..8a5697919 100644 --- a/ikarus/finiteelements/mechanics/nonlinearelastic.hh +++ b/ikarus/finiteelements/mechanics/nonlinearelastic.hh @@ -152,12 +152,7 @@ public: */ template auto materialTangent(const Eigen::Vector& strain) const { - if constexpr (std::is_same_v) - return mat_.template tangentModuli(strain); - else { - decltype(auto) matAD = mat_.template rebind(); - return matAD.template tangentModuli(strain); - } + return material().template tangentModuli(strain); } /** @@ -170,12 +165,7 @@ public: */ template auto getInternalEnergy(const Eigen::Vector& strain) const { - if constexpr (std::is_same_v) - return mat_.template storedEnergy(strain); - else { - decltype(auto) matAD = mat_.template rebind(); - return matAD.template storedEnergy(strain); - } + return material().template storedEnergy(strain); } /** @@ -189,12 +179,7 @@ public: */ template auto getStress(const Eigen::Vector& strain) const { - if constexpr (std::is_same_v) - return mat_.template stresses(strain); - else { - decltype(auto) matAD = mat_.template rebind(); - return matAD.template stresses(strain); - } + return material().template stresses(strain); } const Geometry& geometry() const { return *geo_; } @@ -237,6 +222,14 @@ private: size_t numberOfNodes_{0}; int order_{}; + template + decltype(auto) material() const { + if constexpr (Concepts::AutodiffScalar) + return mat_.template rebind(); + else + return mat_; + } + protected: /** * \brief Calculate the matrix associated with the given Requirement. diff --git a/tests/src/checkfebyautodiff.hh b/tests/src/checkfebyautodiff.hh index cbb0bb055..dc71e8e5c 100644 --- a/tests/src/checkfebyautodiff.hh +++ b/tests/src/checkfebyautodiff.hh @@ -16,7 +16,8 @@ template auto checkFESByAutoDiffImpl(const GridView& gridView, const BasisHandler& basis, Skills&& skills, - AffordanceColl affordance, VectorType& d, const std::string& messageIfFailed = "") { + AffordanceColl affordance, VectorType& d, const std::string& messageIfFailed = "", + double tol = 1e-10) { double lambda = 7.3; auto fe = Ikarus::makeFE(basis, std::forward(skills)); using FE = decltype(fe); @@ -28,8 +29,7 @@ auto checkFESByAutoDiffImpl(const GridView& gridView, const BasisHandler& basis, for (auto element : elements(gridView)) { auto localView = basis.flat().localView(); localView.bind(element); - auto nDOF = localView.size(); - const double tol = 1e-10; + auto nDOF = localView.size(); fe.bind(element); @@ -80,13 +80,14 @@ auto checkFESByAutoDiffImpl(const GridView& gridView, const BasisHandler& basis, template auto checkFESByAutoDiff(const GridView& gridView, const PreBasis& pb, Skills&& skills, AffordanceColl affordance, - const std::string& testName = "") { + const std::string& testName = "", double tol = 1e-10) { Dune::TestSuite t("AutoDiff Test" + testName); auto basis = Ikarus::makeBasis(gridView, pb); Eigen::VectorXd d; d.setZero(basis.flat().dimension()); - t.subTest(checkFESByAutoDiffImpl(gridView, basis, skills, affordance, d, " Zero Displacements")); + t.subTest(checkFESByAutoDiffImpl(gridView, basis, skills, affordance, d, " Zero Displacements", tol)); d.setRandom(basis.flat().dimension()); - t.subTest(checkFESByAutoDiffImpl(gridView, basis, skills, affordance, d, " Non-zero Displacements")); + d *= 0.2; // to avoid a negative determinant of deformation gradient + t.subTest(checkFESByAutoDiffImpl(gridView, basis, skills, affordance, d, " Non-zero Displacements", tol)); return t; } diff --git a/tests/src/testmaterial.cpp b/tests/src/testmaterial.cpp index a856c3789..a67a1c2b8 100644 --- a/tests/src/testmaterial.cpp +++ b/tests/src/testmaterial.cpp @@ -228,6 +228,29 @@ auto testPlaneStrainAgainstPlaneStress(const double tol = 1e-10) { return t; } +template +auto checkThrowNeoHooke(const MAT& matNH) { + static_assert(std::is_same_v>, + "checkThrowNeoHooke is only implemented for Neo-Hooke material law."); + TestSuite t("NeoHooke Test - Checks the throw message for negative determinant of C"); + Eigen::Vector3d E; + E << 2.045327969583023, 0.05875570522766141, 0.3423966429644326; + auto reducedMat = planeStress(matNH, 1e-8); + + t.checkThrow( + [&]() { const auto moduli = (reducedMat.template tangentModuli(E)); }, + "Neo-Hooke test (tangentModuli) should have failed with negative detC for the given E"); + + t.checkThrow( + [&]() { const auto stress = (reducedMat.template stresses(E)); }, + "Neo-Hooke test (stresses) should have failed with negative detC for the given E"); + + t.checkThrow( + [&]() { const auto energy = (reducedMat.template storedEnergy(E)); }, + "Neo-Hooke test (stresses) should have failed with negative detC for the given E"); + return t; +} + int main(int argc, char** argv) { Ikarus::init(argc, argv); TestSuite t; @@ -239,6 +262,7 @@ int main(int argc, char** argv) { auto nh = NeoHooke(matPar); t.subTest(testMaterial(nh)); + t.subTest(checkThrowNeoHooke(nh)); auto le = LinearElasticity(matPar); t.subTest(testMaterial(le)); diff --git a/tests/src/testnonlinearelasticity.hh b/tests/src/testnonlinearelasticity.hh index d83c9106b..ad9768dfc 100644 --- a/tests/src/testnonlinearelasticity.hh +++ b/tests/src/testnonlinearelasticity.hh @@ -290,7 +290,7 @@ auto SingleElementTest(const Material& mat) { } template -void autoDiffTest(TestSuiteType& t, const MAT& mat, const std::string& testName = "") { +void autoDiffTest(TestSuiteType& t, const MAT& mat, const std::string& testName = "", double tol = 1e-10) { using namespace Ikarus; using namespace Dune::Functions::BasisFactory; auto vL = []([[maybe_unused]] const VectorType& globalCoord, auto& lamb) { @@ -317,5 +317,5 @@ void autoDiffTest(TestSuiteType& t, const MAT& mat, const std::string& testName t.subTest(checkFESByAutoDiff( gridView, power(lagrange<1>()), skills(Ikarus::nonLinearElastic(mat), volumeLoad(vL), neumannBoundaryLoad(&neumannBoundary, nBL)), - Ikarus::AffordanceCollections::elastoStatics, testName)); + Ikarus::AffordanceCollections::elastoStatics, testName, tol)); } \ No newline at end of file diff --git a/tests/src/testnonlinearelasticityneohooke.cpp b/tests/src/testnonlinearelasticityneohooke.cpp index 665578858..776260b7d 100644 --- a/tests/src/testnonlinearelasticityneohooke.cpp +++ b/tests/src/testnonlinearelasticityneohooke.cpp @@ -28,8 +28,8 @@ int main(int argc, char** argv) { t.subTest(NonLinearElasticityLoadControlNRandTR(matNH1)); t.subTest(NonLinearElasticityLoadControlNRandTR(matNH1)); - autoDiffTest<2>(t, planeStressMat1, " nu != 0"); - autoDiffTest<2>(t, planeStressMat2, " nu = 0"); + autoDiffTest<2>(t, planeStressMat1, " nu != 0", 1e-7); + autoDiffTest<2>(t, planeStressMat2, " nu = 0", 1e-7); autoDiffTest<3>(t, matNH1, " nu != 0"); autoDiffTest<3>(t, matNH2, " nu = 0"); return t.exit();