fix KL shell

ikarus-project · Jan 10, 2024 · dc4ffaf · dc4ffaf
1 parent 6078e73
commit dc4ffaf
Show file tree

Hide file tree

Showing 7 changed files with 38 additions and 32 deletions.
diff --git a/ikarus/finiteelements/mechanics/membranestrains.hh b/ikarus/finiteelements/mechanics/membranestrains.hh
@@ -40,7 +40,8 @@ namespace Ikarus {
       const Eigen::Matrix<ScalarType, 2, 3> j = J + gradu.transpose();
 
       epsV << J.row(0).dot(gradu.col(0)) + 0.5 * gradu.col(0).squaredNorm(),
-          J.row(1).dot(gradu.col(1)) + 0.5 * gradu.col(1).squaredNorm(), j.row(0).dot(j.row(1));
+          J.row(1).dot(gradu.col(1)) + 0.5 * gradu.col(1).squaredNorm(),
+          j.row(0).dot(j.row(1)) - J.row(0).dot(J.row(1));
       return epsV;
     }
 

diff --git a/ikarus/python/finiteelements/nonlinearelastic.hh b/ikarus/python/finiteelements/nonlinearelastic.hh
@@ -39,8 +39,8 @@ namespace Ikarus::Python {
             }),
             pybind11::keep_alive<1, 2>(), pybind11::keep_alive<1, 3>());
 
-    using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(Eigen::Vector<double, Traits::worlddim>,
-                                                                               const double&)>;
+    using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(
+        Dune::FieldVector<double, Traits::worlddim>, const double&)>;
 
     cls.def(pybind11::init(
                 [](const GlobalBasis& basis, const Element& element, const Material& mat,

diff --git a/ikarus/python/finiteelements/registerelement.hh b/ikarus/python/finiteelements/registerelement.hh
@@ -59,8 +59,8 @@ namespace Ikarus::Python {
               }),
               pybind11::keep_alive<1, 2>(), pybind11::keep_alive<1, 3>());
 
-    using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(Eigen::Vector<double, Traits::worlddim>,
-                                                                               const double&)>;
+    using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(
+        Dune::FieldVector<double, Traits::worlddim>, const double&)>;
     if constexpr (defaultInitializers)
       cls.def(
           pybind11::init([](const GlobalBasis& basis, const Element& element, double emod, double nu,

diff --git a/python/ikarus/finite_elements/__init__.py b/python/ikarus/finite_elements/__init__.py
@@ -104,13 +104,7 @@ def KirchhoffLoveShell(
     includes = ["ikarus/python/finiteelements/kirchhoffloveshell.hh"]
 
     includes += ["ikarus/finiteelements/febases/autodifffe.hh"]
-    autodiffWrapper = "AutoDiffFE"
-    element_type = (
-        "Ikarus::"
-        + autodiffWrapper
-        + f"<Ikarus::KirchhoffLoveShell<{basis.cppTypeName},Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>,"
-        f"Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>"
-    )
+    element_type = f"Ikarus::KirchhoffLoveShell<{basis.cppTypeName},Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>"
 
     # else:
     #     element_type = "Ikarus::" + func.__name__ + f"<{basis.cppTypeName},  {material.cppTypeName} ,Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>"

diff --git a/tests/src/checkfebyautodiff.hh b/tests/src/checkfebyautodiff.hh
@@ -11,17 +11,16 @@
 
 template <template <typename...> class FE, typename GridView, typename PreBasis, typename... ElementArgsType>
 auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const ElementArgsType&... eleArgs) {
-  auto basis     = Ikarus::makeBasis(gridView, pb);
+  auto basis = Ikarus::makeBasis(gridView, pb);
   Eigen::VectorXd d;
   d.setRandom(basis.flat().dimension());
   double lambda = 7.3;
 
   auto req = Ikarus::FErequirements().addAffordance(Ikarus::AffordanceCollections::elastoStatics);
   req.insertGlobalSolution(Ikarus::FESolutions::displacement, d)
       .insertParameter(Ikarus::FEParameter::loadfactor, lambda);
-  Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation"
-                    );
-  for(auto element : elements(gridView)) {
+  Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation");
+  for (auto element : elements(gridView)) {
     auto localView = basis.flat().localView();
     localView.bind(element);
     auto nDOF        = localView.size();
@@ -34,8 +33,6 @@ auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const Eleme
     using AutoDiffBasedFE = Ikarus::AutoDiffFE<decltype(fe), typename decltype(fe)::FERequirementType, false, true>;
     AutoDiffBasedFE feAutoDiff(fe);
 
-
-
     Eigen::MatrixXd K, KAutoDiff;
     K.setZero(nDOF, nDOF);
     KAutoDiff.setZero(nDOF, nDOF);

diff --git a/tests/src/testadaptivestepsizing.cpp b/tests/src/testadaptivestepsizing.cpp
@@ -3,9 +3,9 @@
 
 #include <config.h>
 
+#include "checkfebyautodiff.hh"
 #include "testcommon.hh"
 #include "testhelpers.hh"
-#include "checkfebyautodiff.hh"
 
 #include <dune/common/test/testsuite.hh>
 #include <dune/functions/functionspacebases/boundarydofs.hh>
@@ -112,8 +112,8 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
   for (auto& element : elements(gridView))
     fes.emplace_back(basis, element, E, nu, thickness, utils::LoadDefault{}, &neumannBoundary, neumannBoundaryLoad);
 
-  t.subTest(checkFEByAutoDiff<KirchhoffLoveShellHelper>(gridView, power<3>(nurbs()), E, nu, thickness, utils::LoadDefault{},
-                                                      &neumannBoundary, neumannBoundaryLoad));
+  t.subTest(checkFEByAutoDiff<KirchhoffLoveShellHelper>(gridView, power<3>(nurbs()), E, nu, thickness,
+                                                        utils::LoadDefault{}, &neumannBoundary, neumannBoundaryLoad));
 
   auto basisP = std::make_shared<const decltype(basis)>(basis);
   DirichletValues dirichletValues(basisP->flat());
@@ -166,7 +166,6 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
   auto nonLinOpFull
       = Ikarus::NonLinearOperator(functions(energyFunction, residualFunction, KFunction), parameter(d, lambda));
 
-
   auto nonLinOp  = nonLinOpFull.template subOperator<1, 2>();
   auto linSolver = LinearSolver(SolverTypeTag::sd_SimplicialLDLT);
 
@@ -219,13 +218,16 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
 
   resetNonLinearOperatorParametersToZero(nonLinOp);
   const auto controlInfoWSS = crWSS.run();
-  checkScalars(t, std::ranges::max(d), expectedResults[0][0], message1 + "<Max Displacement>");
-  checkScalars(t, lambda, expectedResults[0][1], message1 + "<Lambda>");
+  const double tolDisp      = 1e-13;
+  const double tolLoad      = 1e-12;
+  checkScalars(t, std::ranges::max(d), expectedResults[0][0], message1 + " <Max Displacement>", tolDisp);
+  checkScalars(t, lambda, expectedResults[0][1], message1 + " <Lambda>", tolLoad);
   resetNonLinearOperatorParametersToZero(nonLinOp);
 
   const auto controlInfoWoSS = crWoSS.run();
-  checkScalars(t, std::ranges::max(d), expectedResults[1][0], message2 + "<Max Displacement>");
-  checkScalars(t, lambda, expectedResults[1][1], message2 + "<Lambda>");
+
+  checkScalars(t, std::ranges::max(d), expectedResults[1][0], message2 + " <Max Displacement>", tolDisp);
+  checkScalars(t, lambda, expectedResults[1][1], message2 + " <Lambda>", tolLoad);
   resetNonLinearOperatorParametersToZero(nonLinOp);
 
   const int controlInfoWSSIterations
@@ -246,6 +248,9 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
   checkSolverInfos(t, expectedIterations[0], controlInfoWSS, loadSteps, message1);
   checkSolverInfos(t, expectedIterations[1], controlInfoWoSS, loadSteps, message2);
 
+  t.check(utils::checkGradient(nonLinOpFull, {.draw = false})) << "Check gradient failed";
+  t.check(utils::checkHessian(nonLinOpFull, {.draw = false})) << "Check hessian failed";
+
   return t;
 }
 
@@ -264,7 +269,7 @@ int main(int argc, char** argv) {
   const std::vector<std::vector<double>> expectedResultsALC
       = {{0.1032139637288574, 0.0003103004514250302}, {0.162759603260405, 0.0007765975850229621}};
   const std::vector<std::vector<double>> expectedResultsLC
-      = {{0.08741028329554587, 0.0002318693543601816}, {0.144353999993308, 6e-4}};
+      = {{0.08741028329554552, 0.0002318693543601816}, {0.144353999993308, 6e-4}};
 
   t.subTest(KLShellAndAdaptiveStepSizing(alc, expectedIterationsALC, expectedResultsALC, 3, 0.4));
   t.subTest(KLShellAndAdaptiveStepSizing(lc, expectedIterationsLC, expectedResultsLC, 2, 1e-4));

diff --git a/tests/src/testhelpers.hh b/tests/src/testhelpers.hh
@@ -31,14 +31,23 @@ requires(std::convertible_to<Derived, Eigen::EigenBase<Derived> const&>and std::
 }
 
 template <typename TestSuiteType, typename ScalarType>
+requires std::is_integral_v<ScalarType>
 void checkScalars(TestSuiteType& t, const ScalarType val, const ScalarType expectedVal,
                   const std::string& messageIfFailed = "") {
   if constexpr (std::is_integral_v<ScalarType>)
-    t.check(val == expectedVal) << std::setprecision(16) << "Incorrect Scalars:\t" << expectedVal << " Actual:\t" << val
-                                << messageIfFailed;
-  else
-    t.check(Dune::FloatCmp::eq(val, expectedVal))
-        << std::setprecision(16) << "Incorrect Scalars:\t" << expectedVal << " Actual:\t" << val << messageIfFailed;
+    t.check(val == expectedVal) << std::setprecision(16) << "Incorrect Scalar integer:\t" << expectedVal << " Actual:\t"
+                                << val << messageIfFailed;
+}
+
+template <typename TestSuiteType, typename ScalarType>
+requires(not std::is_integral_v<ScalarType>) void checkScalars(TestSuiteType& t, const ScalarType val,
+                                                               const ScalarType expectedVal,
+                                                               const std::string& messageIfFailed = "",
+                                                               double tol
+                                                               = Dune::FloatCmp::DefaultEpsilon<ScalarType>::value()) {
+  t.check(Dune::FloatCmp::eq(val, expectedVal, tol))
+      << std::setprecision(16) << "Incorrect Scalar floating point:\t" << expectedVal << " Actual:\t" << val
+      << ". The used tolerance was " << tol << messageIfFailed;
 }
 
 template <typename TestSuiteType, typename ControlInformation>