Use AutoDiffScalar Concept in materials and fix checkFESByAutoDiff te…

…st (#329)

ikarus/finiteelements/mechanics/enhancedassumedstrains.hh

@@ -18,6 +18,7 @@
   #include <ikarus/finiteelements/ferequirements.hh>
   #include <ikarus/finiteelements/mechanics/easvariants.hh>
   #include <ikarus/finiteelements/mixin.hh>
+  #include <ikarus/utils/concepts.hh>
 
 namespace Ikarus {
 
@@ -165,10 +166,10 @@ protected:
       return;
 
     decltype(auto) LMat = [this]() -> decltype(auto) {
-      if constexpr (std::is_same_v<ScalarType, double>)
-        return [this]() -> Eigen::MatrixXd& { return L_; }();
-      else
+      if constexpr (Concepts::AutodiffScalar<ScalarType>)
         return Eigen::MatrixX<ScalarType>{};
+      else
+        return [this]() -> Eigen::MatrixXd& { return L_; }();
     }();
 
     auto calculateMatrixContribution = [&]<typename EAST>(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<ScalarType, double>)
-          return [this]() -> Eigen::MatrixXd& { return L_; }();
-        else
+        if constexpr (Concepts::AutodiffScalar<ScalarType>)
           return Eigen::MatrixX<ScalarType>{};
+        else
+          return [this]() -> Eigen::MatrixXd& { return L_; }();
       }();
       const auto disp  = Dune::viewAsFlatEigenVector(uFunction.coefficientsRef());
       const auto alpha = (-D.inverse() * L_ * disp).eval();

ikarus/finiteelements/mechanics/materials/neohooke.hh

@@ -72,11 +72,13 @@ struct NeoHookeT : public Material<NeoHookeT<ST>>
    * \return ScalarType The stored energy.
    */
   template <typename Derived>
-  ScalarType storedEnergyImpl(const Eigen::MatrixBase<Derived>& C) const noexcept {
+  ScalarType storedEnergyImpl(const Eigen::MatrixBase<Derived>& C) const {
     static_assert(Concepts::EigenMatrixOrVoigtNotation3<Derived>);
     if constexpr (!Concepts::EigenVector<Derived>) {
-      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<NeoHookeT<ST>>
     static_assert(Concepts::EigenMatrixOrVoigtNotation3<Derived>);
     if constexpr (!voigt) {
       if constexpr (!Concepts::EigenVector<Derived>) {
-        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<NeoHookeT<ST>>
   auto tangentModuliImpl(const Eigen::MatrixBase<Derived>& C) const {
     static_assert(Concepts::EigenMatrixOrVoigtNotation3<Derived>);
     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<Eigen::Index, 2>({3, 3}));
       static_assert(Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3>>::NumIndices == 2);
       Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<3, 3, 3, 3>> moduli =
@@ -144,6 +150,13 @@ struct NeoHookeT : public Material<NeoHookeT<ST>>
 
 private:
   MaterialParameters materialParameter_;
+
+  void checkPositiveDetC(ScalarType detC) const {
+    if (Dune::FloatCmp::le(static_cast<double>(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<double>(detC)));
+  }
 };
 
 /**

ikarus/finiteelements/mechanics/nonlinearelastic.hh

@@ -152,12 +152,7 @@ public:
    */
   template <typename ScalarType, int strainDim, bool voigt = true>
   auto materialTangent(const Eigen::Vector<ScalarType, strainDim>& strain) const {
-    if constexpr (std::is_same_v<ScalarType, double>)
-      return mat_.template tangentModuli<strainType, voigt>(strain);
-    else {
-      decltype(auto) matAD = mat_.template rebind<ScalarType>();
-      return matAD.template tangentModuli<strainType, voigt>(strain);
-    }
+    return material<ScalarType>().template tangentModuli<strainType, voigt>(strain);
   }
 
   /**
@@ -170,12 +165,7 @@ public:
    */
   template <typename ScalarType, int strainDim>
   auto getInternalEnergy(const Eigen::Vector<ScalarType, strainDim>& strain) const {
-    if constexpr (std::is_same_v<ScalarType, double>)
-      return mat_.template storedEnergy<strainType>(strain);
-    else {
-      decltype(auto) matAD = mat_.template rebind<ScalarType>();
-      return matAD.template storedEnergy<strainType>(strain);
-    }
+    return material<ScalarType>().template storedEnergy<strainType>(strain);
   }
 
   /**
@@ -189,12 +179,7 @@ public:
    */
   template <typename ScalarType, int strainDim, bool voigt = true>
   auto getStress(const Eigen::Vector<ScalarType, strainDim>& strain) const {
-    if constexpr (std::is_same_v<ScalarType, double>)
-      return mat_.template stresses<strainType, voigt>(strain);
-    else {
-      decltype(auto) matAD = mat_.template rebind<ScalarType>();
-      return matAD.template stresses<strainType, voigt>(strain);
-    }
+    return material<ScalarType>().template stresses<strainType, voigt>(strain);
   }
 
   const Geometry& geometry() const { return *geo_; }
@@ -237,6 +222,14 @@ private:
   size_t numberOfNodes_{0};
   int order_{};
 
+  template <typename ScalarType>
+  decltype(auto) material() const {
+    if constexpr (Concepts::AutodiffScalar<ScalarType>)
+      return mat_.template rebind<ScalarType>();
+    else
+      return mat_;
+  }
+
 protected:
   /**
    * \brief Calculate the matrix associated with the given Requirement.

tests/src/checkfebyautodiff.hh

@@ -16,7 +16,8 @@
 
 template <typename GridView, typename BasisHandler, typename Skills, typename AffordanceColl, typename VectorType>
 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>(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 <typename GridView, typename PreBasis, typename Skills, typename AffordanceColl>
 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;
 }

tests/src/testmaterial.cpp

@@ -228,6 +228,29 @@ auto testPlaneStrainAgainstPlaneStress(const double tol = 1e-10) {
   return t;
 }
 
+template <typename MAT>
+auto checkThrowNeoHooke(const MAT& matNH) {
+  static_assert(std::is_same_v<MAT, NeoHookeT<double>>,
+                "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<Dune::InvalidStateException>(
+      [&]() { const auto moduli = (reducedMat.template tangentModuli<StrainTags::greenLagrangian, true>(E)); },
+      "Neo-Hooke test (tangentModuli) should have failed with negative detC for the given E");
+
+  t.checkThrow<Dune::InvalidStateException>(
+      [&]() { const auto stress = (reducedMat.template stresses<StrainTags::greenLagrangian, true>(E)); },
+      "Neo-Hooke test (stresses) should have failed with negative detC for the given E");
+
+  t.checkThrow<Dune::InvalidStateException>(
+      [&]() { const auto energy = (reducedMat.template storedEnergy<StrainTags::greenLagrangian>(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));

tests/src/testnonlinearelasticity.hh

@@ -290,7 +290,7 @@ auto SingleElementTest(const Material& mat) {
 }
 
 template <int gridDim, typename MAT, typename TestSuiteType>
-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 = []<typename VectorType>([[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<gridDim>(lagrange<1>()),
       skills(Ikarus::nonLinearElastic(mat), volumeLoad<gridDim>(vL), neumannBoundaryLoad(&neumannBoundary, nBL)),
-      Ikarus::AffordanceCollections::elastoStatics, testName));
+      Ikarus::AffordanceCollections::elastoStatics, testName, tol));
 }

tests/src/testnonlinearelasticityneohooke.cpp

@@ -28,8 +28,8 @@ int main(int argc, char** argv) {
   t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::Yasp>(matNH1));
   t.subTest(NonLinearElasticityLoadControlNRandTR<Grids::IgaSurfaceIn2D>(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();