diff --git a/ikarus/finiteelements/mechanics/linearelastic.hh b/ikarus/finiteelements/mechanics/linearelastic.hh
index 33f0b7338..567ce62c3 100644
--- a/ikarus/finiteelements/mechanics/linearelastic.hh
+++ b/ikarus/finiteelements/mechanics/linearelastic.hh
@@ -233,10 +233,10 @@ namespace Ikarus {
     Dune::CachedLocalBasis<
         std::remove_cvref_t<decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis())>>
         localBasis;
-    std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
+    std::function<Eigen::Vector<double, Traits::worlddim>(const Dune::FieldVector<double, Traits::worlddim>&,
                                                           const double&)>
         volumeLoad;
-    std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
+    std::function<Eigen::Vector<double, Traits::worlddim>(const Dune::FieldVector<double, Traits::worlddim>&,
                                                           const double&)>
         neumannBoundaryLoad;
     const BoundaryPatch<GridView>* neumannBoundary;
@@ -270,7 +270,7 @@ namespace Ikarus {
       if (volumeLoad) {
         for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
           const auto uVal                              = u.evaluate(gpIndex);
-          Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(toEigen(geo.global(gp.position())), lambda);
+          Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(geo.global(gp.position()), lambda);
           energy -= uVal.dot(fext) * geo.integrationElement(gp.position()) * gp.weight();
         }
       }
@@ -294,7 +294,7 @@ namespace Ikarus {
           const auto uVal = u.evaluate(quadPos);
 
           // Value of the Neumann data at the current position
-          auto neumannValue = neumannBoundaryLoad(toEigen(intersection.geometry().global(curQuad.position())), lambda);
+          auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
 
           energy -= neumannValue.dot(uVal) * curQuad.weight() * integrationElement;
         }
@@ -327,7 +327,7 @@ namespace Ikarus {
       if (volumeLoad) {
         const auto u = getDisplacementFunction(par, dx);
         for (const auto& [gpIndex, gp] : u.viewOverIntegrationPoints()) {
-          Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(toEigen(geo.global(gp.position())), lambda);
+          Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(geo.global(gp.position()), lambda);
           for (size_t i = 0; i < numberOfNodes; ++i) {
             const auto udCi = u.evaluateDerivative(gpIndex, wrt(coeff(i)));
             force.template segment<myDim>(myDim * i)
@@ -357,8 +357,7 @@ namespace Ikarus {
             const auto udCi = u.evaluateDerivative(quadPos, wrt(coeff(i)));
 
             // Value of the Neumann data at the current position
-            auto neumannValue
-                = neumannBoundaryLoad(toEigen(intersection.geometry().global(curQuad.position())), lambda);
+            auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
             force.template segment<myDim>(myDim * i) -= udCi * neumannValue * curQuad.weight() * integrationElement;
           }
         }
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/ikarus/python/finiteelements/nonlinearelastic.hh b/ikarus/python/finiteelements/nonlinearelastic.hh
index 550014728..db6b58913 100644
--- 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
index b7d1d025d..c9e6cf033 100644
--- 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
index f395d4eac..7072c1dbd 100644
--- 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
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 cd76f11fd..5b232f11d 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>
@@ -31,7 +31,7 @@
 
 using Dune::TestSuite;
 
-template <typename Basis_, typename FERequirements_ = Ikarus::FErequirements<>>
+template <typename Basis_, typename FERequirements_ = Ikarus::FERequirements<>>
 struct KirchhoffLoveShellHelper : Ikarus::KirchhoffLoveShell<Basis_, FERequirements_, false> {
   using Base = Ikarus::KirchhoffLoveShell<Basis_, FERequirements_, false>;
   using Base::Base;
@@ -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/testfeelement.hh b/tests/src/testfeelement.hh
index d82ce4f22..49c2d3ad1 100644
--- a/tests/src/testfeelement.hh
+++ b/tests/src/testfeelement.hh
@@ -43,7 +43,7 @@ auto testFEElement(const PreBasis& preBasis, const std::string& elementName, con
   auto localView = flatBasis.localView();
 
   auto volumeLoad = []<typename VectorType>([[maybe_unused]] const VectorType& globalCoord, auto& lamb) {
-    VectorType fext;
+    Eigen::Vector<typename VectorType::field_type, VectorType::dimension> fext;
     fext.setZero();
     fext[1] = 2 * lamb;
     fext[0] = lamb;
@@ -51,7 +51,7 @@ auto testFEElement(const PreBasis& preBasis, const std::string& elementName, con
   };
 
   auto neumannBoundaryLoad = []<typename VectorType>([[maybe_unused]] const VectorType& globalCoord, auto& lamb) {
-    VectorType fext;
+    Eigen::Vector<typename VectorType::field_type, VectorType::dimension> fext;
     fext.setZero();
     fext[1] = lamb / 40;
     fext[0] = 0;
diff --git a/tests/src/testhelpers.hh b/tests/src/testhelpers.hh
index 6d08fef16..bb282aa3d 100644
--- 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>
diff --git a/tests/src/testklshell.cpp b/tests/src/testklshell.cpp
index 971a2f19d..f77384eeb 100644
--- a/tests/src/testklshell.cpp
+++ b/tests/src/testklshell.cpp
@@ -153,7 +153,7 @@ static auto NonLinearKLShellLoadControlTR() {
   return t;
 }
 
-template <typename Basis_, typename FERequirements_ = Ikarus::FErequirements<>>
+template <typename Basis_, typename FERequirements_ = Ikarus::FERequirements<>>
 struct KirchhoffLoveShellHelper : Ikarus::KirchhoffLoveShell<Basis_, FERequirements_, false> {
   using Base = Ikarus::KirchhoffLoveShell<Basis_, FERequirements_, false>;
   using Base::Base;
diff --git a/tests/src/testnonlinearelasticitysvk.cpp b/tests/src/testnonlinearelasticitysvk.cpp
index c84a686a6..d7fbc9603 100644
--- a/tests/src/testnonlinearelasticitysvk.cpp
+++ b/tests/src/testnonlinearelasticitysvk.cpp
@@ -9,7 +9,7 @@
 
 using Dune::TestSuite;
 
-template <typename Basis_, typename Material, typename FERequirements_ = Ikarus::FErequirements<>>
+template <typename Basis_, typename Material, typename FERequirements_ = Ikarus::FERequirements<>>
 struct NonLinearElasticHelper : Ikarus::NonLinearElastic<Basis_, Material, FERequirements_, false> {
   using Base = Ikarus::NonLinearElastic<Basis_, Material, FERequirements_, false>;
   using Base::Base;