From 22eee422a478d682630a534b90e5a6c4fbb0cc7f Mon Sep 17 00:00:00 2001
From: Alex_Mueller <rathet@gmail.com>
Date: Wed, 10 Jan 2024 12:53:19 +0100
Subject: [PATCH] fix KL shell

---
 .../mechanics/membranestrains.hh               |  3 ++-
 tests/src/checkfebyautodiff.hh                 |  9 +++------
 tests/src/testadaptivestepsizing.cpp           | 18 +++++++++++-------
 tests/src/testhelpers.hh                       | 18 +++++++++++++-----
 4 files changed, 29 insertions(+), 19 deletions(-)

diff --git a/ikarus/finiteelements/mechanics/membranestrains.hh b/ikarus/finiteelements/mechanics/membranestrains.hh
index 66e949731..891ca4e9d 100644
--- 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/tests/src/checkfebyautodiff.hh b/tests/src/checkfebyautodiff.hh
index 6fb079498..230617f42 100644
--- a/tests/src/checkfebyautodiff.hh
+++ b/tests/src/checkfebyautodiff.hh
@@ -11,7 +11,7 @@
 
 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;
@@ -19,9 +19,8 @@ auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const Eleme
   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
index 498af1391..8e6a956a3 100644
--- 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,12 +218,14 @@ 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>");
+  const double tol          = 1e-13;
+  checkScalars(t, std::ranges::max(d), expectedResults[0][0], message1 + " <Max Displacement>", tol);
   checkScalars(t, lambda, expectedResults[0][1], message1 + "<Lambda>");
   resetNonLinearOperatorParametersToZero(nonLinOp);
 
   const auto controlInfoWoSS = crWoSS.run();
-  checkScalars(t, std::ranges::max(d), expectedResults[1][0], message2 + "<Max Displacement>");
+
+  checkScalars(t, std::ranges::max(d), expectedResults[1][0], message2 + " <Max Displacement>", tol);
   checkScalars(t, lambda, expectedResults[1][1], message2 + "<Lambda>");
   resetNonLinearOperatorParametersToZero(nonLinOp);
 
@@ -246,6 +247,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 +268,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
index 6d08fef16..ef5b819e0 100644
--- a/tests/src/testhelpers.hh
+++ b/tests/src/testhelpers.hh
@@ -31,14 +31,22 @@ 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 << messageIfFailed;
 }
 
 template <typename TestSuiteType, typename ControlInformation>