diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9337026e8..1664a56ef 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -84,6 +84,12 @@ SPDX-License-Identifier: LGPL-3.0-or-later
   ([#337](https://github.com/ikarus-project/ikarus/pull/337))
 - `Observers` and `Observables` are replaced with `Broadcasters` and `Listeners`. Existing loggers work almost the same.
  A noteworthy difference is that, the Broadcaster (e.g. a nonlinear solver) has to be registered to a Listener (e.g. a logger) with `logger.subscribeTo(solver)` ([#349](https://github.com/ikarus-project/ikarus/pull/349))
+- Refactor EAS to handle NonLinearElastic ([#325](https://github.com/ikarus-project/ikarus/pull/325))
+    - This also includes a new interface method in the `Mixin` called `updateState` with the respective `impl()` function.
+    - For EAS, `updateStateImpl()` is used to update the internal variable `alpha_` in a nonlinear analysis.
+    - Missing functions like `getStress` and `materialTangentFunction` are added to `LinearElastic` and `NonLinearElastic`, respectively.
+    - `EnhancedAssumedStrains` class now takes a struct denoting the enhanced strain type as a template argument, which implements the
+      enhanced strain value, its first and second derivatives w.r.t the displacements and the internal variable `alpha_`.
 
 ## Release v0.4 (Ganymede)
 
diff --git a/docs/literature.bib b/docs/literature.bib
index 982ab13ba..1bb1dbb49 100644
--- a/docs/literature.bib
+++ b/docs/literature.bib
@@ -107,6 +107,17 @@ @incollection{wunderlich_strategies_1981
 	pages = {63--89},
 }
 
+@article{bischoffShearDeformableShell1997b,
+  title = {{Shear deformable shell elements for large strains and rotations}},
+  author = {Bischoff, M. and Ramm, E.},
+  year = {1997},
+  journal = {International Journal for Numerical Methods in Engineering},
+  volume = {40},
+  number = {23},
+  pages = {4427-4449},
+  doi = {10.1002/(SICI)1097-0207(19971215)40:23<4427::AID-NME268>3.0.CO;2-9},
+}
+
 @article{riks_application_1972,
 	title = {The {Application} of {Newton}’s {Method} to the {Problem} of {Elastic} {Stability}},
 	volume = {39},
@@ -210,6 +221,14 @@ @book{bonet2008nonlinear
 	doi  = {https://doi.org/10.1017/CBO9780511755446}
 }
 
+@phdthesis{bieberLockingHourglassingNonlinear2024,
+  title = {Locking and Hourglassing in Nonlinear Finite Element Technology},
+  author = {Bieber, Simon},
+  year = {2024},
+  address = {Stuttgart},
+  langid = {english},
+  school = {Institut f{\"u}r Baustatik und Baudynamik, Universit{\"a}t Stuttgart}
+}
 
 @book{trustregion,
 	author = {Conn, Andrew R. and Gould, Nicholas I. M. and Toint, Philippe L.},
diff --git a/ikarus/finiteelements/ferequirements.hh b/ikarus/finiteelements/ferequirements.hh
index b314eccf2..edb5f97d1 100644
--- a/ikarus/finiteelements/ferequirements.hh
+++ b/ikarus/finiteelements/ferequirements.hh
@@ -204,7 +204,7 @@ auto scalarAffordance(VectorAffordance affordanceV) {
 namespace AffordanceCollections {
   inline constexpr AffordanceCollection elastoStatics(ScalarAffordance::mechanicalPotentialEnergy,
                                                       VectorAffordance::forces, MatrixAffordance::stiffness);
-}
+} // namespace AffordanceCollections
 
 /**
  * \class FERequirements
diff --git a/ikarus/finiteelements/mechanics/CMakeLists.txt b/ikarus/finiteelements/mechanics/CMakeLists.txt
index 68c915320..d13ba1faf 100644
--- a/ikarus/finiteelements/mechanics/CMakeLists.txt
+++ b/ikarus/finiteelements/mechanics/CMakeLists.txt
@@ -1,6 +1,6 @@
 # SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
 # SPDX-License-Identifier: LGPL-3.0-or-later
-add_subdirectory(eas)
+add_subdirectory(strainenhancements)
 add_subdirectory(loads)
 add_subdirectory(materials)
 # install headers
@@ -12,7 +12,6 @@ install(
         kirchhoffloveshell.hh
         membranestrains.hh
         loads.hh
-        easvariants.hh
         truss.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements/mechanics
 )
diff --git a/ikarus/finiteelements/mechanics/eas/eas2d.hh b/ikarus/finiteelements/mechanics/eas/eas2d.hh
deleted file mode 100644
index 54cf1bf29..000000000
--- a/ikarus/finiteelements/mechanics/eas/eas2d.hh
+++ /dev/null
@@ -1,157 +0,0 @@
-// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
-// SPDX-License-Identifier: LGPL-3.0-or-later
-
-/**
- * \file eas2d.hh
- * \brief Definition of the types of EAS formulations for 2D elements.
- * \ingroup  eas
- */
-
-#pragma once
-
-#include <ikarus/utils/tensorutils.hh>
-
-namespace Ikarus::EAS {
-
-/**
- * \brief Dummy struct for displacement-based EAS elements, i.e. 0 enhanced modes
- */
-template <typename GEO>
-struct E0
-{
-  static constexpr int strainSize         = GEO::mydimension * (GEO::mydimension + 1) / 2;
-  static constexpr int enhancedStrainSize = 0;
-  using MType                             = Eigen::Matrix<double, strainSize, enhancedStrainSize>;
-  using DType                             = Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>;
-
-  E0() = default;
-  explicit E0(const GEO& /*geometry*/) {}
-
-  // returns an Eigen zero expression for optimization purposes
-  static auto calcM(const Dune::FieldVector<double, 2>& /*quadPos*/) { return MType::Zero(); }
-};
-
-/**
- * \brief Q1E4 structure for EAS with linear strains and 4 enhanced modes.
- *
- * \details The Q1E4 structure represents an implementation of EAS for a
- * specific case where linear strains are considered, and the method includes
- * enhancement with 4 additional modes to improve the accuracy of finite element solutions.
- *
- * \tparam GEO The geometry type.
- */
-template <typename GEO>
-struct Q1E4
-{
-  static constexpr int strainSize         = 3;
-  static constexpr int enhancedStrainSize = 4;
-  using MType                             = Eigen::Matrix<double, strainSize, enhancedStrainSize>;
-  using DType                             = Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>;
-
-  Q1E4() = default;
-  explicit Q1E4(const GEO& geometry)
-      : geometry_{std::make_optional<GEO>(geometry)},
-        T0InverseTransformed_{calcTransformationMatrix2D(geometry)} {}
-
-  auto calcM(const Dune::FieldVector<double, 2>& quadPos) const {
-    MType M;
-    M.setZero(strainSize, enhancedStrainSize);
-    const double xi   = quadPos[0];
-    const double eta  = quadPos[1];
-    M(0, 0)           = 2 * xi - 1.0;
-    M(1, 1)           = 2 * eta - 1.0;
-    M(2, 2)           = 2 * xi - 1.0;
-    M(2, 3)           = 2 * eta - 1.0;
-    const double detJ = geometry_->integrationElement(quadPos);
-    M                 = T0InverseTransformed_ / detJ * M;
-    return M;
-  }
-
-private:
-  std::optional<GEO> geometry_;
-  Eigen::Matrix3d T0InverseTransformed_;
-};
-
-/**
- * \brief Structure representing EAS for Q1 with 5 enhanced strains.
- *
- * This structure defines the EAS for Q1 elements with 5 enhanced strains.
- *
- * \tparam GEO The geometry type.
- */
-template <typename GEO>
-struct Q1E5
-{
-  static constexpr int strainSize         = 3;
-  static constexpr int enhancedStrainSize = 5;
-  using MType                             = Eigen::Matrix<double, strainSize, enhancedStrainSize>;
-  using DType                             = Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>;
-
-  Q1E5() = default;
-  explicit Q1E5(const GEO& geometry)
-      : geometry_{std::make_optional<GEO>(geometry)},
-        T0InverseTransformed_{calcTransformationMatrix2D(geometry)} {}
-
-  auto calcM(const Dune::FieldVector<double, 2>& quadPos) const {
-    MType M;
-    M.setZero();
-    const double xi   = quadPos[0];
-    const double eta  = quadPos[1];
-    M(0, 0)           = 2 * xi - 1.0;
-    M(1, 1)           = 2 * eta - 1.0;
-    M(2, 2)           = 2 * xi - 1.0;
-    M(2, 3)           = 2 * eta - 1.0;
-    M(2, 4)           = (2 * xi - 1.0) * (2 * eta - 1.0);
-    const double detJ = geometry_->integrationElement(quadPos);
-    M                 = T0InverseTransformed_ / detJ * M;
-    return M;
-  }
-
-private:
-  std::optional<GEO> geometry_;
-  Eigen::Matrix3d T0InverseTransformed_;
-};
-
-/**
- * \brief Structure representing EAS for Q1 with 7 enhanced strains.
- *
- * This structure defines the EAS for Q1 elements with 7 enhanced strains.
- *
- * \tparam GEO The geometry type.
- */
-template <typename GEO>
-struct Q1E7
-{
-  static constexpr int strainSize         = 3;
-  static constexpr int enhancedStrainSize = 7;
-  using MType                             = Eigen::Matrix<double, strainSize, enhancedStrainSize>;
-  using DType                             = Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>;
-
-  Q1E7() = default;
-  explicit Q1E7(const GEO& geometry)
-      : geometry_{std::make_optional<GEO>(geometry)},
-        T0InverseTransformed_{calcTransformationMatrix2D(geometry)} {}
-
-  auto calcM(const Dune::FieldVector<double, 2>& quadPos) const {
-    MType M;
-    M.setZero();
-    const double xi   = quadPos[0];
-    const double eta  = quadPos[1];
-    M(0, 0)           = 2 * xi - 1.0;
-    M(1, 1)           = 2 * eta - 1.0;
-    M(2, 2)           = 2 * xi - 1.0;
-    M(2, 3)           = 2 * eta - 1.0;
-    M(0, 4)           = (2 * xi - 1.0) * (2 * eta - 1.0);
-    M(1, 5)           = (2 * xi - 1.0) * (2 * eta - 1.0);
-    M(2, 6)           = (2 * xi - 1.0) * (2 * eta - 1.0);
-    const double detJ = geometry_->integrationElement(quadPos);
-    M                 = T0InverseTransformed_ / detJ * M;
-    return M;
-  }
-
-private:
-  std::optional<GEO> geometry_;
-  Eigen::Matrix3d T0InverseTransformed_;
-};
-
-} // namespace Ikarus::EAS
diff --git a/ikarus/finiteelements/mechanics/eas/eas3d.hh b/ikarus/finiteelements/mechanics/eas/eas3d.hh
deleted file mode 100644
index 62440c13d..000000000
--- a/ikarus/finiteelements/mechanics/eas/eas3d.hh
+++ /dev/null
@@ -1,120 +0,0 @@
-// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
-// SPDX-License-Identifier: LGPL-3.0-or-later
-
-/**
- * \file eas3d.hh
- * \brief Definition of the types of EAS formulations for 3D elements.
- * \ingroup  mechanics
- */
-
-#pragma once
-
-#include <ikarus/utils/tensorutils.hh>
-
-namespace Ikarus::EAS {
-
-/**
- * \brief Structure representing EAS for H1 with 9 enhanced strains.
- *
- * This structure defines the EAS for H1 elements with 9 enhanced strains.
- *
- * \tparam GEO The geometry type.
- */
-template <typename GEO>
-struct H1E9
-{
-  static constexpr int strainSize         = 6;
-  static constexpr int enhancedStrainSize = 9;
-  using MType                             = Eigen::Matrix<double, strainSize, enhancedStrainSize>;
-  using DType                             = Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>;
-
-  H1E9() = default;
-  explicit H1E9(const GEO& geometry)
-      : geometry_{std::make_optional<GEO>(geometry)},
-        T0InverseTransformed_{calcTransformationMatrix3D(geometry)} {}
-
-  auto calcM(const Dune::FieldVector<double, 3>& quadPos) const {
-    MType M;
-    M.setZero();
-    const double xi   = quadPos[0];
-    const double eta  = quadPos[1];
-    const double zeta = quadPos[2];
-    M(0, 0)           = 2 * xi - 1.0;
-    M(1, 1)           = 2 * eta - 1.0;
-    M(2, 2)           = 2 * zeta - 1.0;
-    M(3, 3)           = 2 * xi - 1.0;
-    M(3, 4)           = 2 * eta - 1.0;
-    M(4, 5)           = 2 * xi - 1.0;
-    M(4, 6)           = 2 * zeta - 1.0;
-    M(5, 7)           = 2 * eta - 1.0;
-    M(5, 8)           = 2 * zeta - 1.0;
-    const double detJ = geometry_->integrationElement(quadPos);
-    M                 = T0InverseTransformed_ / detJ * M;
-    return M;
-  }
-
-private:
-  std::optional<GEO> geometry_;
-  Eigen::Matrix<double, strainSize, strainSize> T0InverseTransformed_;
-};
-
-/**
- * \brief Structure representing EAS for H1 with 21 enhanced strains.
- *
- * This structure defines the EAS for H1 elements with 21 enhanced strains.
- *
- * \tparam GEO The geometry type.
- */
-template <typename GEO>
-struct H1E21
-{
-  static constexpr int strainSize         = 6;
-  static constexpr int enhancedStrainSize = 21;
-  using MType                             = Eigen::Matrix<double, strainSize, enhancedStrainSize>;
-  using DType                             = Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>;
-
-  H1E21() = default;
-  explicit H1E21(const GEO& geometry)
-      : geometry_{std::make_optional<GEO>(geometry)},
-        T0InverseTransformed{calcTransformationMatrix3D(geometry)} {}
-
-  auto calcM(const Dune::FieldVector<double, 3>& quadPos) const {
-    MType M;
-    M.setZero();
-    const double xi   = quadPos[0];
-    const double eta  = quadPos[1];
-    const double zeta = quadPos[2];
-    M(0, 0)           = 2 * xi - 1.0;
-    M(1, 1)           = 2 * eta - 1.0;
-    M(2, 2)           = 2 * zeta - 1.0;
-    M(3, 3)           = 2 * xi - 1.0;
-    M(3, 4)           = 2 * eta - 1.0;
-    M(4, 5)           = 2 * xi - 1.0;
-    M(4, 6)           = 2 * zeta - 1.0;
-    M(5, 7)           = 2 * eta - 1.0;
-    M(5, 8)           = 2 * zeta - 1.0;
-
-    M(3, 9)  = (2 * xi - 1.0) * (2 * zeta - 1.0);
-    M(3, 10) = (2 * eta - 1.0) * (2 * zeta - 1.0);
-    M(4, 11) = (2 * xi - 1.0) * (2 * eta - 1.0);
-    M(4, 12) = (2 * eta - 1.0) * (2 * zeta - 1.0);
-    M(5, 13) = (2 * xi - 1.0) * (2 * eta - 1.0);
-    M(5, 14) = (2 * xi - 1.0) * (2 * zeta - 1.0);
-
-    M(0, 15) = (2 * xi - 1.0) * (2 * eta - 1.0);
-    M(0, 16) = (2 * xi - 1.0) * (2 * zeta - 1.0);
-    M(1, 17) = (2 * xi - 1.0) * (2 * eta - 1.0);
-    M(1, 18) = (2 * eta - 1.0) * (2 * zeta - 1.0);
-    M(2, 19) = (2 * xi - 1.0) * (2 * zeta - 1.0);
-    M(2, 20) = (2 * eta - 1.0) * (2 * zeta - 1.0);
-
-    const double detJ = geometry_->integrationElement(quadPos);
-    M                 = T0InverseTransformed / detJ * M;
-    return M;
-  }
-
-private:
-  std::optional<GEO> geometry_;
-  Eigen::Matrix<double, strainSize, strainSize> T0InverseTransformed;
-};
-} // namespace Ikarus::EAS
diff --git a/ikarus/finiteelements/mechanics/easvariants.hh b/ikarus/finiteelements/mechanics/easvariants.hh
deleted file mode 100644
index 92e340f0a..000000000
--- a/ikarus/finiteelements/mechanics/easvariants.hh
+++ /dev/null
@@ -1,119 +0,0 @@
-// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
-// SPDX-License-Identifier: LGPL-3.0-or-later
-
-/**
- * \file easvariants.hh
- * \brief Definition of the EAS variants.
- * \ingroup  mechanics
- */
-
-#pragma once
-
-#include <ikarus/finiteelements/mechanics/eas/eas2d.hh>
-#include <ikarus/finiteelements/mechanics/eas/eas3d.hh>
-#include <ikarus/utils/tensorutils.hh>
-
-namespace Ikarus::EAS {
-namespace Impl {
-  /**
-   * \brief EASVariants structure holding variants of different Enhanced Assumed Strains (EAS).
-   *
-   * \details EASVariants holds the different variants of EAS formulations for
-   * both 2D and 3D elements.
-   *
-   * \tparam GEO The geometry type.
-   */
-  template <typename GEO>
-  struct Variants
-  {
-    static constexpr int worldDim = GEO::coorddimension;
-    static_assert((worldDim == 2) or (worldDim == 3), "EAS variants are only available 2D and 3D elements.");
-
-    /** \brief Variant type holding different EAS formulations for 2D elements. */
-    using EAS2D = std::variant<E0<GEO>, Q1E4<GEO>, Q1E5<GEO>, Q1E7<GEO>>;
-
-    /** \brief Variant type holding different EAS formulations for 3D elements. */
-    using EAS3D = std::variant<E0<GEO>, H1E9<GEO>, H1E21<GEO>>;
-
-    /** \brief Type of the EAS variant depending on the grid dimension. */
-    using type = std::conditional_t<worldDim == 2, EAS2D, EAS3D>;
-  };
-
-  /**
-   * \brief Wrapper around the EAS variant, contains helper functions
-   * \tparam Geometry the Geometry type of the udderlying grid element
-   */
-  template <typename Geometry>
-  struct EASVariant
-  {
-    using Variant = Impl::Variants<Geometry>::type;
-
-    /**
-     * \brief Executes a function F on the variant, if `enhancedStrainSize` is not zero
-     * \tparam F the Type of the function F
-     * \param f the function, which takes the eas element as an argument
-     */
-    template <typename F>
-    void operator()(F&& f) const {
-      std::visit(
-          [&]<typename EAST>(const EAST& easFunction) {
-            if constexpr (EAST::enhancedStrainSize != 0)
-              f(easFunction);
-          },
-          var_);
-    }
-
-    auto numberOfEASParameters() const {
-      return std::visit([]<typename EAST>(const EAST&) -> int { return EAST::enhancedStrainSize; }, var_);
-    }
-    bool isDisplacmentBased() const { return numberOfEASParameters() == 0; }
-    void setEASType(int numberOfEASParameters) {
-      numberOfEASParameters_ = numberOfEASParameters;
-      if (geometry_)
-        createEASType();
-    }
-    void bind(const Geometry& geometry) {
-      geometry_ = std::make_optional<Geometry>(geometry);
-      createEASType();
-    }
-
-  private:
-    void createEASType() {
-      if (numberOfEASParameters_ == 0) {
-        var_ = E0(geometry_.value());
-        return;
-      }
-
-      if constexpr (Geometry::mydimension == 2) {
-        switch (numberOfEASParameters_) {
-          case 4:
-            var_ = Q1E4(geometry_.value());
-            break;
-          case 5:
-            var_ = Q1E5(geometry_.value());
-            break;
-          case 7:
-            var_ = Q1E7(geometry_.value());
-            break;
-          default:
-            DUNE_THROW(Dune::NotImplemented, "The given EAS parameters are not available for the 2D case.");
-        }
-      } else if constexpr (Geometry::mydimension == 3) {
-        switch (numberOfEASParameters_) {
-          case 9:
-            var_ = H1E9(geometry_.value());
-            break;
-          case 21:
-            var_ = H1E21(geometry_.value());
-            break;
-          default:
-            DUNE_THROW(Dune::NotImplemented, "The given EAS parameters are not available for the 3D case.");
-        }
-      }
-    }
-    std::optional<Geometry> geometry_;
-    Variant var_;
-    int numberOfEASParameters_;
-  };
-} // namespace Impl
-} // namespace Ikarus::EAS
diff --git a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh
index 1e69a1ff5..01e01cbe0 100644
--- a/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh
+++ b/ikarus/finiteelements/mechanics/enhancedassumedstrains.hh
@@ -4,7 +4,7 @@
 /**
  * \file enhancedassumedstrains.hh
  * \brief Definition of the EAS class.
- * \ingroup  mechanics
+ * \ingroup mechanics
  */
 
 #pragma once
@@ -16,24 +16,27 @@
 
   #include <ikarus/finiteelements/fehelper.hh>
   #include <ikarus/finiteelements/ferequirements.hh>
-  #include <ikarus/finiteelements/mechanics/easvariants.hh>
-  #include <ikarus/finiteelements/mixin.hh>
+  #include <ikarus/finiteelements/mechanics/materials/tags.hh>
+  #include <ikarus/finiteelements/mechanics/strainenhancements/easvariants.hh>
+  #include <ikarus/utils/broadcaster/broadcastermessages.hh>
   #include <ikarus/utils/concepts.hh>
 
 namespace Ikarus {
 
-template <typename PreFE, typename FE>
+template <typename PreFE, typename FE, typename ES>
 class EnhancedAssumedStrains;
 
 /**
  * \brief A PreFE struct for Enhanced Assumed Strains.
+ * \tparam ES The enhanced strain type.
  */
+template <typename ES>
 struct EnhancedAssumedStrainsPre
 {
   int numberOfParameters{};
 
   template <typename PreFE, typename FE>
-  using Skill = EnhancedAssumedStrains<PreFE, FE>;
+  using Skill = EnhancedAssumedStrains<PreFE, FE, ES>;
 };
 
 /**
@@ -45,23 +48,42 @@ struct EnhancedAssumedStrainsPre
  *
  * \tparam PreFE Type of the pre finite element.
  * \tparam FE Type of the finite element.
+ * \tparam ES The enhanced strain type.
  */
-template <typename PreFE, typename FE>
+template <typename PreFE, typename FE, typename ES>
 class EnhancedAssumedStrains
+    : public std::conditional_t<std::same_as<ES, EAS::LinearStrain>,
+                                ResultTypeBase<ResultTypes::linearStress, ResultTypes::linearStressFull>,
+                                ResultTypeBase<ResultTypes::PK2Stress, ResultTypes::PK2StressFull>>
 {
 public:
-  using Traits      = PreFE::Traits;
-  using Requirement = FERequirements<FESolutions::displacement, FEParameter::loadfactor>;
-  using LocalView   = typename Traits::LocalView;
-  using Geometry    = typename Traits::Geometry;
-  using GridView    = typename Traits::GridView;
-  using Pre         = EnhancedAssumedStrainsPre;
+  using Traits         = PreFE::Traits;
+  using Requirement    = FERequirements<FESolutions::displacement, FEParameter::loadfactor>;
+  using LocalView      = typename Traits::LocalView;
+  using Geometry       = typename Traits::Geometry;
+  using GridView       = typename Traits::GridView;
+  using Pre            = EnhancedAssumedStrainsPre<ES>;
+  using EnhancedStrain = ES;
+
+  static constexpr int myDim = Traits::mydim;
+
+  template <typename ST>
+  using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
+
+  template <template <typename, int, int> class RT>
+  using RTWrapperType = ResultWrapper<RT<typename Traits::ctype, Traits::mydim, Traits::worlddim>, ResultShape::Vector>;
 
   /**
    * \brief Constructor for Enhanced Assumed Strains elements.
    * \param pre The pre finite element
    */
-  explicit EnhancedAssumedStrains(const Pre& pre) { this->setEASType(pre.numberOfParameters); }
+  explicit EnhancedAssumedStrains(const Pre& pre) {
+    static_assert(Concepts::Formulations::TotalLagrangian<FE::strainType, FE::stressType>,
+                  "EAS method is only implemented for the total Lagrangian setting.");
+    static_assert(not(FE::strainType == StrainTags::linear) or (std::same_as<EnhancedStrain, EAS::LinearStrain>),
+                  "If FE::strainType is linear, then enhancedStrain must also be linear.");
+    this->setEASType(pre.numberOfParameters);
+  }
 
   /**
    * \brief Checks if the element is displacement-based and the EAS is turned off.
@@ -99,29 +121,28 @@ public:
    * \tparam RT The type representing the requested result.
    */
   template <template <typename, int, int> class RT>
+  requires(EnhancedAssumedStrains::template supportsResultType<RT>())
   auto calculateAtImpl(const Requirement& req, const Dune::FieldVector<double, Traits::mydim>& local,
                        Dune::PriorityTag<2>) const {
-    if (isDisplacementBased())
-      return underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());
-
-    using namespace Dune::Indices;
-    using namespace Dune::DerivativeDirections;
-    auto resultVector = underlying().template calculateAtImpl<RT>(req, local, Dune::PriorityTag<1>());
-    if constexpr (isSameResultType<RT, ResultTypes::linearStress>) {
-      using RTWrapper = ResultWrapper<RT<typename Traits::ctype, Traits::mydim, Traits::worlddim>, ResultShape::Vector>;
-      RTWrapper resultWrapper;
-
+    if constexpr (isSameResultType<RT, ResultTypes::linearStress> or isSameResultType<RT, ResultTypes::PK2Stress> or
+                  isSameResultType<RT, ResultTypes::linearStressFull> or
+                  isSameResultType<RT, ResultTypes::PK2StressFull>) {
+      const auto ufunc     = underlying().displacementFunction(req);
+      const auto rFunction = underlying().template resultFunction<RT>();
+      auto disp            = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
+      const auto geo       = underlying().localView().element().geometry();
+
+      RTWrapperType<RT> resultWrapper{};
       auto calculateAtContribution = [&]<typename EAST>(const EAST& easFunction) {
-        typename EAST::DType D;
-        calculateDAndLMatrix(easFunction, req, D, L_);
-        const auto ufunc = underlying().displacementFunction(req);
-        const auto disp  = Dune::viewAsFlatEigenVector(ufunc.coefficientsRef());
-        const auto C     = underlying().materialTangentFunction(req);
-        const auto alpha = (-D.inverse() * L_ * disp).eval();
-        const auto M     = easFunction.calcM(local);
-        const auto CEval = C(local);
-        auto easStress   = (CEval * M * alpha).eval();
-        resultWrapper    = resultVector.asVec() + easStress;
+        Eigen::VectorXd alpha;
+        alpha.setZero(numberOfEASParameters());
+        if constexpr (EAST::enhancedStrainSize != 0) {
+          typename EAST::DType D;
+          calculateDAndLMatrix(easFunction, req, D, L_);
+          alpha = -D.inverse() * L_ * disp;
+        }
+        const auto enhancedStrain = EnhancedStrain::value(geo, ufunc, local, easFunction, alpha);
+        resultWrapper             = rFunction(enhancedStrain);
       };
       easVariant_(calculateAtContribution);
       return resultWrapper;
@@ -138,144 +159,245 @@ public:
     if (numberOfEASParameters != 0)
       easApplicabilityCheck();
     easVariant_.setEASType(numberOfEASParameters);
+    initializeState();
   }
 
+  /**
+   * \brief Gets the internal state variable alpha for the EAS element.
+   *
+   * \return Internal state variable (alpha).
+   */
+  const auto& alpha() const { return alpha_; }
+
 protected:
   void bindImpl() {
     assert(underlying().localView().bound());
     easVariant_.bind(underlying().localView().element().geometry());
+    initializeState();
+  }
+
+  /**
+   * \brief Updates the internal state variable alpha_ at the end of an iteration before the update of the displacements
+   * done by the non-linear solver. See \cite bieberLockingHourglassingNonlinear2024 for implementation details and
+   * further references.
+   *
+   * \param par The Requirement object.
+   * \param correction The correction in displacement (DeltaD) vector passed based on which the internal state variable
+   * alpha is to be updated.
+   */
+  void updateStateImpl(const Requirement& par,
+                       const std::remove_reference_t<typename Traits::template VectorType<>>& correction) {
+    using ScalarType = Traits::ctype;
+    easApplicabilityCheck();
+    auto correctAlpha = [&]<typename EAST>(const EAST& easFunction) {
+      if constexpr (EAST::enhancedStrainSize != 0) {
+        const auto& Rtilde      = calculateRtilde<ScalarType>(par);
+        const auto localdxBlock = Ikarus::FEHelper::localSolutionBlockVector<Traits, Eigen::VectorXd, double>(
+            correction, underlying().localView());
+        const auto localdx               = Dune::viewAsFlatEigenVector(localdxBlock);
+        constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
+        Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize> D;
+        calculateDAndLMatrix(easFunction, par, D, L_);
+        this->alpha_ -= D.inverse() * (Rtilde + (L_ * localdx));
+      }
+    };
+
+    easVariant_(correctAlpha);
   }
 
-public:
-protected:
   inline void easApplicabilityCheck() const {
     const auto& numberOfNodes = underlying().numberOfNodes();
-    assert(not(not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 8 and Traits::mydim == 3)) and
-               (not isDisplacementBased())) &&
-           "EAS only supported for Q1 or H1 elements");
+    if (not((numberOfNodes == 4 and Traits::mydim == 2) or (numberOfNodes == 9 and Traits::mydim == 2) or
+            (numberOfNodes == 8 and Traits::mydim == 3)) and
+        (not isDisplacementBased()))
+      DUNE_THROW(Dune::NotImplemented, "EAS is only supported for Q1, Q2 and H1 elements");
   }
 
   template <typename ScalarType>
-  void calculateMatrixImpl(
-      const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<> K,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
-    using namespace Dune::DerivativeDirections;
-    using namespace Dune;
+  void calculateMatrixImpl(const Requirement& par, const MatrixAffordance& affordance,
+                           typename Traits::template MatrixType<> K,
+                           const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
+    if (affordance != MatrixAffordance::stiffness)
+      DUNE_THROW(Dune::NotImplemented, "MatrixAffordance not implemented: " + toString(affordance));
     easApplicabilityCheck();
-    if (isDisplacementBased())
-      return;
 
-    decltype(auto) LMat = [this]() -> decltype(auto) {
-      if constexpr (Concepts::AutodiffScalar<ScalarType>)
-        return Eigen::MatrixX<ScalarType>{};
-      else
-        return [this]() -> Eigen::MatrixXd& { return L_; }();
-    }();
+    const auto geo           = underlying().localView().element().geometry();
+    const auto& uFunction    = underlying().displacementFunction(par, dx);
+    const auto& kFunction    = underlying().template stiffnessMatrixFunction<ScalarType>(par, K, dx);
+    const auto numberOfNodes = underlying().numberOfNodes();
+    const auto& localBasis   = underlying().localBasis();
 
     auto calculateMatrixContribution = [&]<typename EAST>(const EAST& easFunction) {
-      typename EAST::DType D;
-      calculateDAndLMatrix(easFunction, par, D, LMat, dx);
+      for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
+        const auto enhancedStrain = EnhancedStrain::value(geo, uFunction, gp.position(), easFunction, alpha_);
+        const auto stresses       = underlying().stress(enhancedStrain);
+        for (size_t i = 0; i < numberOfNodes; ++i) {
+          const auto E_dI = EnhancedStrain::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
+                                                                        gp.position(), easFunction, alpha_, i);
+          for (size_t j = 0; j < numberOfNodes; ++j) {
+            const auto E_dJ = EnhancedStrain::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
+                                                                          gp.position(), easFunction, alpha_, j);
+            const auto E_dd = EnhancedStrain::template secondDerivative<0>(
+                geo, uFunction, localBasis, gpIndex, gp.position(), stresses, easFunction, alpha_, i, j);
+            kFunction(enhancedStrain, E_dI, E_dJ, E_dd, i, j, gp);
+          }
+        }
+      }
 
-      K.template triangularView<Eigen::Upper>() -= LMat.transpose() * D.inverse() * LMat;
-      K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
+      if constexpr (EAST::enhancedStrainSize != 0) {
+        typename EAST::DType D;
+        calculateDAndLMatrix(easFunction, par, D, L_);
+        K.template triangularView<Eigen::Upper>() -= L_.transpose() * D.inverse() * L_;
+        K.template triangularView<Eigen::StrictlyLower>() = K.transpose();
+      }
     };
     easVariant_(calculateMatrixContribution);
   }
 
   template <typename ScalarType>
-  inline ScalarType calculateScalarImpl(
-      const Requirement& par, ScalarAffordance affordance,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  inline ScalarType calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
+                                        const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
     easApplicabilityCheck();
     if (isDisplacementBased())
-      return 0.0;
+      return underlying().template energyFunction<ScalarType>(par, dx)();
+
     DUNE_THROW(Dune::NotImplemented,
                "EAS element do not support any scalar calculations, i.e. they are not derivable from a potential");
   }
 
   template <typename ScalarType>
-  void calculateVectorImpl(
-      const Requirement& par, VectorAffordance affordance, typename Traits::template VectorType<ScalarType> force,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  void calculateVectorImpl(const Requirement& par, VectorAffordance affordance,
+                           typename Traits::template VectorType<ScalarType> force,
+                           const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
+    if (affordance != VectorAffordance::forces)
+      DUNE_THROW(Dune::NotImplemented, "VectorAffordance not implemented: " + toString(affordance));
     easApplicabilityCheck();
-    using namespace Dune;
-    using namespace Dune::DerivativeDirections;
-    const auto uFunction      = underlying().displacementFunction(par, dx);
-    auto strainFunction       = underlying().strainFunction(par, dx);
-    const auto& numberOfNodes = underlying().numberOfNodes();
+
+    const auto geo           = underlying().localView().element().geometry();
+    const auto& uFunction    = underlying().displacementFunction(par, dx);
+    const auto numberOfNodes = underlying().numberOfNodes();
+    const auto& localBasis   = underlying().localBasis();
+    const auto fIntFunction  = underlying().template internalForcesFunction<ScalarType>(par, force, dx);
 
     auto calculateForceContribution = [&]<typename EAST>(const EAST& easFunction) {
-      typename EAST::DType D;
-      calculateDAndLMatrix(easFunction, par, D, L_);
-
-      decltype(auto) LMat = [this]() -> decltype(auto) {
-        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();
-      const auto geo   = underlying().localView().element().geometry();
-      auto C           = underlying().materialTangentFunction(par);
-
-      for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
-        const auto M            = easFunction.calcM(gp.position());
-        const double intElement = geo.integrationElement(gp.position()) * gp.weight();
-        const auto CEval        = C(gpIndex);
-        auto stresses           = (CEval * M * alpha).eval();
+      for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
+        const auto enhancedStrain = EnhancedStrain::value(geo, uFunction, gp.position(), easFunction, alpha_);
         for (size_t i = 0; i < numberOfNodes; ++i) {
-          const auto bopI = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
-          force.template segment<Traits::worlddim>(Traits::worlddim * i) += bopI.transpose() * stresses * intElement;
+          const auto E_dI = EnhancedStrain::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
+                                                                        gp.position(), easFunction, alpha_, i);
+          fIntFunction(enhancedStrain, E_dI, i, gp);
         }
       }
+      if constexpr (EAST::enhancedStrainSize != 0) {
+        // Here, L_ and D doesn't see the ScalarType, while Rtilde sees it. This is needed because only then when we use
+        // automatic differentiation w.r.t to the displacements, we get the correct static-condensed stiffness matrix.
+        // If L_ and D sees the ScalarType, then their derivatives w.r.t. the displacements is also performed which is
+        // then not equivalent to the static-condensed stiffness matrix in general.
+        typename EAST::DType D;
+        calculateDAndLMatrix(easFunction, par, D, L_);
+        const auto& Rtilde = calculateRtilde(par, dx);
+        force -= L_.transpose() * D.inverse() * Rtilde;
+      }
     };
     easVariant_(calculateForceContribution);
   }
 
+  template <typename MT, typename BC>
+  void subscribeToImpl(BC& bc) {
+    if constexpr (std::same_as<MT, NonLinearSolverMessages>) {
+      using NLSState = typename BC::State;
+      underlying().subscribe(bc, [&](NonLinearSolverMessages message, const NLSState& state) {
+        if (message == NonLinearSolverMessages::CORRECTION_UPDATED)
+          this->updateStateImpl(state.domain, state.correction);
+      });
+    }
+  }
+
 private:
-  EAS::Impl::EASVariant<Geometry> easVariant_;
+  EAS::EASVariant<ES, Geometry> easVariant_;
   mutable Eigen::MatrixXd L_;
+  Eigen::VectorXd alpha_;
 
   //> CRTP
   const auto& underlying() const { return static_cast<const FE&>(*this); }
   auto& underlying() { return static_cast<FE&>(*this); }
 
-  template <typename ScalarType, int enhancedStrainSize>
-  void calculateDAndLMatrix(
-      const auto& easFunction, const auto& par, Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>& DMat,
-      Eigen::MatrixX<ScalarType>& LMat,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  /**
+   * \brief Initializes the internal state variable alpha_ based on the number of EAS parameters.
+   */
+  void initializeState() {
+    alpha_.resize(numberOfEASParameters());
+    alpha_.setZero();
+  }
+
+  template <int enhancedStrainSize>
+  void calculateDAndLMatrix(const auto& easFunction, const auto& par,
+                            Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>& DMat,
+                            Eigen::MatrixX<double>& LMat) const {
     using namespace Dune;
     using namespace Dune::DerivativeDirections;
 
-    auto strainFunction      = underlying().strainFunction(par, dx);
-    const auto C             = underlying().materialTangentFunction(par);
+    const auto& uFunction    = underlying().displacementFunction(par);
     const auto geo           = underlying().localView().element().geometry();
     const auto numberOfNodes = underlying().numberOfNodes();
+    const auto& localBasis   = underlying().localBasis();
     DMat.setZero();
     LMat.setZero(enhancedStrainSize, underlying().localView().size());
-    for (const auto& [gpIndex, gp] : strainFunction.viewOverIntegrationPoints()) {
-      const auto M            = easFunction.calcM(gp.position());
-      const auto CEval        = C(gpIndex);
-      const double detJTimesW = geo.integrationElement(gp.position()) * gp.weight();
-      DMat += M.transpose() * CEval * M * detJTimesW;
+    for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
+      const auto enhancedStrain = EnhancedStrain::value(geo, uFunction, gp.position(), easFunction, alpha_);
+      const auto C              = underlying().materialTangent(enhancedStrain);
+      const auto stresses       = underlying().stress(enhancedStrain);
+      const double intElement   = geo.integrationElement(gp.position()) * gp.weight();
+      const auto E_a  = EnhancedStrain::template firstDerivative<1>(geo, uFunction, localBasis, gpIndex, gp.position(),
+                                                                    easFunction, alpha_);
+      const auto E_aa = EnhancedStrain::template secondDerivative<1>(geo, uFunction, localBasis, gpIndex, gp.position(),
+                                                                     stresses, easFunction, alpha_);
+      DMat += (E_a.transpose() * C * E_a + E_aa) * intElement;
       for (size_t i = 0U; i < numberOfNodes; ++i) {
-        const size_t I = Traits::worlddim * i;
-        const auto Bi  = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
-        LMat.template block<enhancedStrainSize, Traits::worlddim>(0, I) += M.transpose() * CEval * Bi * detJTimesW;
+        const auto E_dI = EnhancedStrain::template firstDerivative<0>(geo, uFunction, localBasis, gpIndex,
+                                                                      gp.position(), easFunction, alpha_, i);
+        const auto E_ad = EnhancedStrain::template secondDerivative<2>(geo, uFunction, localBasis, gpIndex,
+                                                                       gp.position(), stresses, easFunction, alpha_, i);
+        LMat.template block<enhancedStrainSize, myDim>(0, myDim * i) +=
+            (E_a.transpose() * C * E_dI + E_ad) * intElement;
       }
     }
   }
+
+  template <typename ScalarType>
+  Eigen::VectorX<ScalarType> calculateRtilde(const Requirement& par,
+                                             const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
+    const auto geo         = underlying().localView().element().geometry();
+    const auto& uFunction  = underlying().displacementFunction(par, dx);
+    const auto& localBasis = underlying().localBasis();
+    Eigen::VectorX<ScalarType> Rtilde;
+    Rtilde.setZero(numberOfEASParameters());
+
+    auto calculateRtildeContribution = [&]<typename EAST>(const EAST& easFunction) {
+      for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
+        const auto enhancedStrain = EnhancedStrain::value(geo, uFunction, gp.position(), easFunction, alpha_);
+        const auto E_a = EnhancedStrain::template firstDerivative<1>(geo, uFunction, localBasis, gpIndex, gp.position(),
+                                                                     easFunction, alpha_);
+        const double intElement = geo.integrationElement(gp.position()) * gp.weight();
+        auto stresses           = underlying().stress(enhancedStrain);
+        Rtilde += (E_a.transpose() * stresses).eval() * intElement;
+      }
+    };
+
+    easVariant_(calculateRtildeContribution);
+    return Rtilde;
+  }
 };
 
 /**
  * \brief A helper function to create an enhanced assumed strain pre finite element.
+ * \tparam ES The enhanced strain type.
  * \param numberOfEASParameters Number of EAS parameters
  * \return An enhanced assumed strain pre finite element.
  */
+template <typename ES = EAS::LinearStrain>
 auto eas(int numberOfEASParameters = 0) {
-  EnhancedAssumedStrainsPre pre(numberOfEASParameters);
+  EnhancedAssumedStrainsPre<ES> pre(numberOfEASParameters);
 
   return pre;
 }
diff --git a/ikarus/finiteelements/mechanics/linearelastic.hh b/ikarus/finiteelements/mechanics/linearelastic.hh
index a1b503bcb..0b4e334a4 100644
--- a/ikarus/finiteelements/mechanics/linearelastic.hh
+++ b/ikarus/finiteelements/mechanics/linearelastic.hh
@@ -22,6 +22,8 @@
   #include <ikarus/finiteelements/ferequirements.hh>
   #include <ikarus/finiteelements/feresulttypes.hh>
   #include <ikarus/finiteelements/mechanics/materials.hh>
+  #include <ikarus/finiteelements/mechanics/materials/tags.hh>
+  #include <ikarus/finiteelements/mechanics/strainenhancements/easfunctions.hh>
   #include <ikarus/finiteelements/physicshelper.hh>
 
 namespace Ikarus {
@@ -66,8 +68,28 @@ public:
   using Material     = PRE::Material;
   using Pre          = PRE;
 
-  static constexpr int myDim = Traits::mydim;
-  using LocalBasisType       = decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis());
+  using LocalBasisType = decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis());
+
+  template <typename ST>
+  using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
+
+  static constexpr int myDim       = Traits::mydim;
+  static constexpr int strainDim   = myDim * (myDim + 1) / 2;
+  static constexpr auto strainType = StrainTags::linear;
+  static constexpr auto stressType = StressTags::linear;
+  static constexpr bool hasEAS     = FE::template hasEAS<EAS::LinearStrain>;
+
+  template <template <typename, int, int> class RT>
+  using RTWrapperType = ResultWrapper<RT<typename Traits::ctype, myDim, Traits::worlddim>, ResultShape::Vector>;
+
+  template <typename ST>
+  using StrainType = Eigen::Vector<ST, strainDim>;
+
+  template <typename ST>
+  using BopType = Eigen::Matrix<ST, strainDim, myDim>;
+
+  template <typename ST>
+  using KgType = Eigen::Matrix<ST, myDim, myDim>;
 
   /**
    * \brief Constructor for the LinearElastic class.
@@ -109,9 +131,7 @@ public:
    * \return The displacement function.
    */
   template <typename ScalarType = double>
-  auto displacementFunction(
-      const Requirement& par,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  auto displacementFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
     const auto& d = par.globalSolution();
     auto disp     = Ikarus::FEHelper::localSolutionBlockVector<Traits>(d, underlying().localView(), dx);
     Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_);
@@ -127,35 +147,59 @@ public:
    * \return The strain function.
    */
   template <class ScalarType = double>
-  auto strainFunction(
-      const Requirement& par,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  auto strainFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
     return Dune::linearStrains(displacementFunction(par, dx));
   }
 
   /**
-   * \brief Gets the material tangent matrix for the linear elastic material.
+   * \brief Get the internal energy for the given strain.
    *
-   * \return The material tangent matrix.
+   * \tparam ScalarType The scalar type for the material and strain.
+   * \tparam strainDim The dimension of the strain vector.
+   * \param strain The strain vector in Voigt notation.
+   * \return The internal energy calculated using the material's storedEnergy function.
    */
-  auto materialTangent() const {
-    // Since that material is independent of the strains, a zero strain is passed here
-    return mat_.template tangentModuli<StrainTags::linear, true>(Eigen::Vector<double, 6>::Zero());
+  template <typename ScalarType, int strainDim>
+  auto internalEnergy(const Eigen::Vector<ScalarType, strainDim>& strain) const {
+    const auto C = materialTangent<ScalarType, strainDim>(strain);
+    return (0.5 * strain.dot(C * strain));
   }
 
   /**
-   * \brief Gets the material tangent function for the given Requirement.
+   * \brief Get the stress for the given strain.
    *
-   * \param par The Requirement object.
-   * \return The material tangent function.
+   * \tparam ScalarType The scalar type for the material and strain.
+   * \tparam strainDim The dimension of the strain vector.
+   * \tparam voigt A boolean indicating whether to use the Voigt notation for stress.
+   * \param strain The strain vector in Voigt notation.
+   * \return The stress vector calculated using the material's stresses function.
+   */
+  template <typename ScalarType, int strainDim, bool voigt = true>
+  auto stress(const Eigen::Vector<ScalarType, strainDim>& strain) const {
+    const auto C = materialTangent<ScalarType, strainDim, voigt>(strain);
+    return (C * strain).eval();
+  }
+
+  /**
+   * \brief Get the linear elastic material tangent for the given strain for the given Requirement.
+   *
+   * \tparam ScalarType The scalar type for the material and strain.
+   * \tparam strainDim The dimension of the strain vector.
+   * \tparam voigt Flag indicating whether to use Voigt notation.
+   * \param strain The strain vector.
+   * \return The material tangent calculated using the material's tangentModuli function.
    */
-  auto materialTangentFunction([[maybe_unused]] const Requirement& par) const {
-    return [&]([[maybe_unused]] auto gp) { return materialTangent(); };
+  template <typename ScalarType, int strainDim, bool voigt = true>
+  auto materialTangent(
+      const Eigen::Vector<ScalarType, strainDim>& strain = Eigen::Vector<ScalarType, strainDim>::Zero()) const {
+    // Since that material is independent of the strains, a zero strain is passed here
+    return material<ScalarType>().template tangentModuli<strainType, voigt>(strain);
   }
 
   const Geometry& geometry() const { return *geo_; }
   [[nodiscard]] size_t numberOfNodes() const { return numberOfNodes_; }
   [[nodiscard]] int order() const { return order_; }
+  const Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>>& localBasis() const { return localBasis_; }
 
   template <typename ScalarType = double>
   decltype(auto) material() const {
@@ -166,33 +210,48 @@ public:
   }
 
 public:
+  /**
+   * \brief Get a lambda function that evaluates the requested result type for a given strain (in Voigt notation).
+   * \tparam RT The type representing the requested result.
+   * \return A lambda function that evaluates the requested result type for a given strain (in Voigt notation).
+   */
+  template <template <typename, int, int> class RT>
+  requires(LinearElastic::template supportsResultType<RT>())
+  auto resultFunction() const {
+    return [&](const Eigen::Vector<double, strainDim>& strainInVoigt) {
+      if constexpr (isSameResultType<RT, ResultTypes::linearStress> or
+                    isSameResultType<RT, ResultTypes::linearStressFull>) {
+        decltype(auto) mat = [&]() {
+          if constexpr (isSameResultType<RT, ResultTypes::linearStressFull> and requires { mat_.underlying(); })
+            return mat_.underlying();
+          else
+            return mat_;
+        }();
+        return RTWrapperType<RT>{mat.template stresses<StrainTags::linear>(enlargeIfReduced<Material>(strainInVoigt))};
+      }
+    };
+  }
+
   /**
    * \brief Calculates a requested result at a specific local position.
-   *
    * \param req The Requirement object holding the global solution.
    * \param local Local position vector.
-   * \tparam RT The requested result type
-   * \return calculated result
+   * \return Calculated result.
    *
    * \tparam RT The type representing the requested result.
    */
   template <template <typename, int, int> class RT>
-  requires(supportsResultType<RT>())
+  requires(LinearElastic::template supportsResultType<RT>())
   auto calculateAtImpl(const Requirement& req, const Dune::FieldVector<double, Traits::mydim>& local,
                        Dune::PriorityTag<1>) const {
-    using RTWrapper = ResultWrapper<RT<typename Traits::ctype, myDim, Traits::worlddim>, ResultShape::Vector>;
-
+    if constexpr (hasEAS)
+      return RTWrapperType<RT>{};
     if constexpr (isSameResultType<RT, ResultTypes::linearStress> or
                   isSameResultType<RT, ResultTypes::linearStressFull>) {
-      const auto eps     = strainFunction(req);
-      auto epsVoigt      = eps.evaluate(local, Dune::on(Dune::DerivativeDirections::gridElement));
-      decltype(auto) mat = [&]() -> decltype(auto) {
-        if constexpr (isSameResultType<RT, ResultTypes::linearStressFull> and requires { mat_.underlying(); })
-          return mat_.underlying();
-        else
-          return mat_;
-      }();
-      return RTWrapper{mat.template stresses<StrainTags::linear>(enlargeIfReduced<Material>(epsVoigt))};
+      const auto rFunction = resultFunction<RT>();
+      const auto eps       = strainFunction(req);
+      auto epsVoigt        = eps.evaluate(local, Dune::on(Dune::DerivativeDirections::gridElement));
+      return rFunction(epsVoigt);
     }
   }
 
@@ -207,23 +266,92 @@ private:
   size_t numberOfNodes_{0};
   int order_{};
 
+public:
+  /**
+   * \brief Get a lambda function that evaluates the stiffness matrix for a given strain, integration point and its
+   * index.
+   *
+   * \tparam ST The scalar type for the material and strain.
+   * \param par The Requirement object.
+   * \param dx Optional displacement vector.
+   * \param K The matrix to store the calculated result.
+   * \return A lambda function that evaluates the stiffness matrix for a given strain, integration point and its index.
+   */
+  template <typename ST>
+  auto stiffnessMatrixFunction(const Requirement& par, typename Traits::template MatrixType<ST>& K,
+                               const VectorXOptRef<ST>& dx = std::nullopt) const {
+    return [&](const StrainType<ST>& strain, const BopType<ST>& bopI, const BopType<ST>& bopJ, const KgType<ST>& kgIJ,
+               const int I, const int J, const auto& gp) {
+      const auto C            = materialTangent(strain);
+      const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
+      K.template block<myDim, myDim>(I * myDim, J * myDim) += (bopI.transpose() * C * bopJ) * intElement;
+    };
+  }
+
+  /**
+   * \brief Get a lambda function that evaluates the internal force vector for a given strain, integration point and its
+   * index.
+   *
+   * \tparam ST The scalar type for the material and strain.
+   * \param par The Requirement object.
+   * \param dx Optional displacement vector.
+   * \param force The vector to store the calculated result.
+   * \return A lambda function that evaluates the intenral force vector for a given strain, integration point and its
+   * index.
+   */
+  template <typename ST>
+  auto internalForcesFunction(const Requirement& par, typename Traits::template VectorType<ST>& force,
+                              const VectorXOptRef<ST>& dx = std::nullopt) const {
+    return [&](const StrainType<ST>& strain, const BopType<ST>& bopI, const int I, const auto& gp) {
+      const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
+      auto stresses           = stress(strain);
+      force.template segment<myDim>(myDim * I) += bopI.transpose() * stresses * intElement;
+    };
+  }
+
+  /**
+   * \brief Get a lambda function that evaluates the internal energy at a given integration point and its index.
+   *
+   * \tparam ScalarType The scalar type for the material and strain.
+   * \param par The Requirement object.
+   * \param dx Optional displacement vector.
+   * \return A lambda function that returns the intenral energy at a given integration point and its index.
+   */
+  template <typename ScalarType>
+  auto energyFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
+    return [&]() -> ScalarType {
+      using namespace Dune::DerivativeDirections;
+      using namespace Dune;
+      ScalarType energy = 0.0;
+      const auto eps    = strainFunction(par, dx);
+      for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
+        const auto epsVoigt = eps.evaluate(gpIndex, on(gridElement));
+        energy += internalEnergy(epsVoigt) * geo_->integrationElement(gp.position()) * gp.weight();
+      }
+      return energy;
+    };
+  }
+
 protected:
   template <typename ScalarType>
-  void calculateMatrixImpl(
-      const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<> K,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
-    const auto eps = strainFunction(par, dx);
+  void calculateMatrixImpl(const Requirement& par, const MatrixAffordance& affordance,
+                           typename Traits::template MatrixType<> K,
+                           const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
+    if constexpr (hasEAS)
+      return;
     using namespace Dune::DerivativeDirections;
     using namespace Dune;
+    const auto eps       = strainFunction(par, dx);
+    const auto kFunction = stiffnessMatrixFunction<ScalarType>(par, K, dx);
+    const auto& kgIJ     = KgType<ScalarType>::Zero().eval();
 
-    const auto C = materialTangent();
     for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
-      const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
+      const auto epsVoigt = eps.evaluate(gpIndex, on(gridElement));
       for (size_t i = 0; i < numberOfNodes_; ++i) {
         const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
         for (size_t j = 0; j < numberOfNodes_; ++j) {
           const auto bopJ = eps.evaluateDerivative(gpIndex, wrt(coeff(j)), on(gridElement));
-          K.template block<myDim, myDim>(i * myDim, j * myDim) += bopI.transpose() * C * bopJ * intElement;
+          kFunction(epsVoigt, bopI, bopJ, kgIJ, i, j, gp);
         }
       }
     }
@@ -231,41 +359,29 @@ protected:
 
   template <typename ScalarType>
   auto calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
-                           const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx =
-                               std::nullopt) const -> ScalarType {
-    const auto uFunction = displacementFunction(par, dx);
-    const auto eps       = strainFunction(par, dx);
-    const auto& lambda   = par.parameter();
-    using namespace Dune::DerivativeDirections;
-    using namespace Dune;
-
-    const auto C = materialTangent();
-
-    ScalarType energy = 0.0;
-    for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
-      const auto EVoigt = eps.evaluate(gpIndex, on(gridElement));
-      energy += (0.5 * EVoigt.dot(C * EVoigt)) * geo_->integrationElement(gp.position()) * gp.weight();
-    }
-    return energy;
+                           const VectorXOptRef<ScalarType>& dx = std::nullopt) const -> ScalarType {
+    if constexpr (hasEAS)
+      return ScalarType{0.0};
+    return energyFunction(par, dx)();
   }
 
   template <typename ScalarType>
-  void calculateVectorImpl(
-      const Requirement& par, VectorAffordance affordance, typename Traits::template VectorType<ScalarType> force,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
-    const auto eps = strainFunction(par, dx);
+  void calculateVectorImpl(const Requirement& par, VectorAffordance affordance,
+                           typename Traits::template VectorType<ScalarType> force,
+                           const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
+    if constexpr (hasEAS)
+      return;
     using namespace Dune::DerivativeDirections;
     using namespace Dune;
-
-    const auto C = materialTangent();
+    const auto eps          = strainFunction(par, dx);
+    const auto fIntFunction = internalForcesFunction<ScalarType>(par, force, dx);
 
     // Internal forces
     for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
-      const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
-      auto stresses           = (C * eps.evaluate(gpIndex, on(gridElement))).eval();
+      const auto epsVoigt = eps.evaluate(gpIndex, on(gridElement));
       for (size_t i = 0; i < numberOfNodes_; ++i) {
         const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
-        force.template segment<myDim>(myDim * i) += bopI.transpose() * stresses * intElement;
+        fIntFunction(epsVoigt, bopI, i, gp);
       }
     }
   }
diff --git a/ikarus/finiteelements/mechanics/materials/hyperelastic/deviatoric/blatzko.hh b/ikarus/finiteelements/mechanics/materials/hyperelastic/deviatoric/blatzko.hh
index 063bc1401..463150468 100644
--- a/ikarus/finiteelements/mechanics/materials/hyperelastic/deviatoric/blatzko.hh
+++ b/ikarus/finiteelements/mechanics/materials/hyperelastic/deviatoric/blatzko.hh
@@ -82,8 +82,8 @@ struct BlatzKoT
    */
   template <typename ST>
   FirstDerivative<ST> firstDerivativeImpl(const PrincipalStretches<ST>& lambda) const {
-    auto dWdLambda     = FirstDerivative<ST>::Zero().eval();
-    const ScalarType J = Impl::determinantFromPrincipalValues(lambda);
+    auto dWdLambda = FirstDerivative<ST>::Zero().eval();
+    const ST J     = Impl::determinantFromPrincipalValues(lambda);
 
     return mu_ * (-lambda.cwisePow(3).cwiseInverse() + (J * lambda.cwiseInverse()));
   }
@@ -97,8 +97,8 @@ struct BlatzKoT
    */
   template <typename ST>
   SecondDerivative<ST> secondDerivativeImpl(const PrincipalStretches<ST>& lambda) const {
-    auto dS            = SecondDerivative<ST>::Zero().eval();
-    const ScalarType J = Impl::determinantFromPrincipalValues(lambda);
+    auto dS    = SecondDerivative<ST>::Zero().eval();
+    const ST J = Impl::determinantFromPrincipalValues(lambda);
 
     auto lam      = lambda.array();
     dS            = J / (lambda * lambda.transpose()).array();
diff --git a/ikarus/finiteelements/mechanics/materials/linearelasticity.hh b/ikarus/finiteelements/mechanics/materials/linearelasticity.hh
index 4cbb8c122..a64659d2a 100644
--- a/ikarus/finiteelements/mechanics/materials/linearelasticity.hh
+++ b/ikarus/finiteelements/mechanics/materials/linearelasticity.hh
@@ -115,8 +115,7 @@ struct LinearElasticityT : Material<LinearElasticityT<ST>>
    */
   template <typename STO>
   auto rebind() const {
-    LinearElasticityT<STO> leRebound;
-    leRebound.svk = svk_.template rebind<STO>();
+    LinearElasticityT<STO> leRebound{svk_.template rebind<STO>().materialParametersImpl()};
     return leRebound;
   }
 
diff --git a/ikarus/finiteelements/mechanics/membranestrains.hh b/ikarus/finiteelements/mechanics/membranestrains.hh
index c9f1dfcdd..1c463ef13 100644
--- a/ikarus/finiteelements/mechanics/membranestrains.hh
+++ b/ikarus/finiteelements/mechanics/membranestrains.hh
@@ -7,7 +7,6 @@
  */
 
 #pragma once
-#include <ranges>
 
 #include <dune/common/fvector.hh>
 
diff --git a/ikarus/finiteelements/mechanics/nonlinearelastic.hh b/ikarus/finiteelements/mechanics/nonlinearelastic.hh
index 64c1de3ba..e7bf14d0a 100644
--- a/ikarus/finiteelements/mechanics/nonlinearelastic.hh
+++ b/ikarus/finiteelements/mechanics/nonlinearelastic.hh
@@ -23,6 +23,7 @@
   #include <ikarus/finiteelements/ferequirements.hh>
   #include <ikarus/finiteelements/mechanics/loads.hh>
   #include <ikarus/finiteelements/mechanics/materials/tags.hh>
+  #include <ikarus/finiteelements/mechanics/strainenhancements/easfunctions.hh>
   #include <ikarus/finiteelements/physicshelper.hh>
   #include <ikarus/utils/defaultfunctions.hh>
   #include <ikarus/utils/eigendunetransformations.hh>
@@ -73,8 +74,26 @@ public:
 
   using LocalBasisType = decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis());
 
+  template <typename ST>
+  using VectorXOptRef = std::optional<std::reference_wrapper<const Eigen::VectorX<ST>>>;
+
   static constexpr int myDim       = Traits::mydim;
+  static constexpr int strainDim   = myDim * (myDim + 1) / 2;
   static constexpr auto strainType = StrainTags::greenLagrangian;
+  static constexpr auto stressType = StressTags::PK2;
+  static constexpr bool hasEAS     = FE::template hasEAS<EAS::GreenLagrangeStrain>;
+
+  template <template <typename, int, int> class RT>
+  using RTWrapperType = ResultWrapper<RT<typename Traits::ctype, myDim, Traits::worlddim>, ResultShape::Vector>;
+
+  template <typename ST>
+  using StrainType = Eigen::Vector<ST, strainDim>;
+
+  template <typename ST>
+  using BopType = Eigen::Matrix<ST, strainDim, myDim>;
+
+  template <typename ST>
+  using KgType = Eigen::Matrix<ST, myDim, myDim>;
 
   /**
    * \brief Constructor for the NonLinearElastic class.
@@ -116,9 +135,7 @@ public:
    * \return A StandardLocalFunction representing the displacement function.
    */
   template <typename ScalarType = double>
-  auto displacementFunction(
-      const Requirement& par,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  auto displacementFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
     const auto& d = par.globalSolution();
     auto disp     = Ikarus::FEHelper::localSolutionBlockVector<Traits>(d, underlying().localView(), dx);
     Dune::StandardLocalFunction uFunction(localBasis_, disp, geo_);
@@ -134,26 +151,10 @@ public:
    * \return The strain function calculated using greenLagrangeStrains.
    */
   template <typename ScalarType = double>
-  inline auto strainFunction(
-      const Requirement& par,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  inline auto strainFunction(const Requirement& par, const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
     return Dune::greenLagrangeStrains(displacementFunction(par, dx));
   }
 
-  /**
-   * \brief Get the material tangent for the given strain.
-   *
-   * \tparam ScalarType The scalar type for the material and strain.
-   * \tparam strainDim The dimension of the strain vector.
-   * \tparam voigt Flag indicating whether to use Voigt notation.
-   * \param strain The strain vector.
-   * \return The material tangent calculated using the material's tangentModuli function.
-   */
-  template <typename ScalarType, int strainDim, bool voigt = true>
-  auto materialTangent(const Eigen::Vector<ScalarType, strainDim>& strain) const {
-    return material<ScalarType>().template tangentModuli<strainType, voigt>(strain);
-  }
-
   /**
    * \brief Get the internal energy for the given strain.
    *
@@ -163,7 +164,7 @@ public:
    * \return The internal energy calculated using the material's storedEnergy function.
    */
   template <typename ScalarType, int strainDim>
-  auto getInternalEnergy(const Eigen::Vector<ScalarType, strainDim>& strain) const {
+  auto internalEnergy(const Eigen::Vector<ScalarType, strainDim>& strain) const {
     return material<ScalarType>().template storedEnergy<strainType>(strain);
   }
 
@@ -177,13 +178,28 @@ public:
    * \return The stress vector calculated using the material's stresses function.
    */
   template <typename ScalarType, int strainDim, bool voigt = true>
-  auto getStress(const Eigen::Vector<ScalarType, strainDim>& strain) const {
+  auto stress(const Eigen::Vector<ScalarType, strainDim>& strain) const {
     return material<ScalarType>().template stresses<strainType, voigt>(strain);
   }
 
+  /**
+   * \brief Get the material tangent for the given strain for the given Requirement.
+   *
+   * \tparam ScalarType The scalar type for the material and strain.
+   * \tparam strainDim The dimension of the strain vector.
+   * \tparam voigt Flag indicating whether to use Voigt notation.
+   * \param strain The strain vector.
+   * \return The material tangent calculated using the material's tangentModuli function.
+   */
+  template <typename ScalarType, int strainDim, bool voigt = true>
+  auto materialTangent(const Eigen::Vector<ScalarType, strainDim>& strain) const {
+    return material<ScalarType>().template tangentModuli<strainType, voigt>(strain);
+  }
+
   const Geometry& geometry() const { return *geo_; }
   [[nodiscard]] size_t numberOfNodes() const { return numberOfNodes_; }
   [[nodiscard]] int order() const { return order_; }
+  const Dune::CachedLocalBasis<std::remove_cvref_t<LocalBasisType>>& localBasis() const { return localBasis_; }
 
   template <typename ScalarType = double>
   decltype(auto) material() const {
@@ -194,6 +210,27 @@ public:
   }
 
 public:
+  /**
+   * \brief Get a lambda function that evaluates the requested result type for a given strain (in Voigt notation).
+   * \tparam RT The type representing the requested result.
+   * \return A lambda function that evaluates the requested result type for a given strain (in Voigt notation).
+   */
+  template <template <typename, int, int> class RT>
+  requires(supportsResultType<RT>())
+  auto resultFunction() const {
+    return [&](const Eigen::Vector<double, strainDim>& strainInVoigt) {
+      if constexpr (isSameResultType<RT, ResultTypes::PK2Stress> or isSameResultType<RT, ResultTypes::PK2StressFull>) {
+        decltype(auto) mat = [&]() {
+          if constexpr (isSameResultType<RT, ResultTypes::PK2StressFull> and requires { mat_.underlying(); })
+            return mat_.underlying();
+          else
+            return mat_;
+        }();
+        return RTWrapperType<RT>{mat.template stresses<strainType>(enlargeIfReduced<Material>(strainInVoigt))};
+      }
+    };
+  }
+
   /**
    * \brief Calculates a requested result at a specific local position.
    *
@@ -209,19 +246,15 @@ public:
                        Dune::PriorityTag<1>) const {
     using namespace Dune::DerivativeDirections;
 
-    using RTWrapper = ResultWrapper<RT<typename Traits::ctype, myDim, Traits::worlddim>, ResultShape::Vector>;
+    if constexpr (hasEAS)
+      return RTWrapperType<RT>{};
     if constexpr (isSameResultType<RT, ResultTypes::PK2Stress> or isSameResultType<RT, ResultTypes::PK2StressFull>) {
       const auto uFunction = displacementFunction(req);
+      const auto rFunction = resultFunction<RT>();
       const auto H         = uFunction.evaluateDerivative(local, Dune::wrt(spatialAll), Dune::on(gridElement));
       const auto E         = (0.5 * (H.transpose() + H + H.transpose() * H)).eval();
 
-      decltype(auto) mat = [&]() {
-        if constexpr (isSameResultType<RT, ResultTypes::PK2StressFull> and requires { mat_.underlying(); })
-          return mat_.underlying();
-        else
-          return mat_;
-      }();
-      return RTWrapper{mat.template stresses<StrainTags::greenLagrangian>(enlargeIfReduced<Material>(toVoigt(E)))};
+      return rFunction(toVoigt(E));
     }
   }
 
@@ -235,6 +268,72 @@ private:
   size_t numberOfNodes_{0};
   int order_{};
 
+public:
+  /**
+   * \brief Get a lambda function that evaluates the stiffness matrix for a given strain, integration point and its
+   * index.
+   *
+   * \tparam ST The scalar type for the material and strain.
+   * \param par The Requirement object.
+   * \param dx Optional displacement vector.
+   * \param K The matrix to store the calculated result.
+   * \return A lambda function that evaluates the stiffness matrix for a given strain, integration point and its index.
+   */
+  template <typename ST>
+  auto stiffnessMatrixFunction(const Requirement& par, typename Traits::template MatrixType<ST>& K,
+                               const VectorXOptRef<ST>& dx = std::nullopt) const {
+    return [&](const StrainType<ST>& strain, const BopType<ST>& bopI, const BopType<ST>& bopJ, const KgType<ST>& kgIJ,
+               const int I, const int J, const auto& gp) {
+      const auto C            = materialTangent(strain);
+      const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
+      K.template block<myDim, myDim>(I * myDim, J * myDim) += (bopI.transpose() * C * bopJ + kgIJ) * intElement;
+    };
+  }
+
+  /**
+   * \brief Get a lambda function that evaluates the internal force vector for a given strain, integration point and its
+   * index.
+   *
+   * \tparam ST The scalar type for the material and strain.
+   * \param par The Requirement object.
+   * \param dx Optional displacement vector.
+   * \param force The vector to store the calculated result.
+   * \return A lambda function that evaluates the intenral force vector for a given strain, integration point and its
+   * index.
+   */
+  template <typename ST>
+  auto internalForcesFunction(const Requirement& par, typename Traits::template VectorType<ST>& force,
+                              const VectorXOptRef<ST>& dx = std::nullopt) const {
+    return [&](const StrainType<ST>& strain, const BopType<ST>& bopI, const int I, const auto& gp) {
+      const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
+      auto stresses           = stress(strain);
+      force.template segment<myDim>(myDim * I) += bopI.transpose() * stresses * intElement;
+    };
+  }
+
+  /**
+   * \brief Get a lambda function that evaluates the internal energy at a given integration point and its index.
+   *
+   * \tparam ST The scalar type for the material and strain.
+   * \param par The Requirement object.
+   * \param dx Optional displacement vector.
+   * \return A lambda function that returns the intenral energy at a given integration point and its index.
+   */
+  template <typename ST>
+  auto energyFunction(const Requirement& par, const VectorXOptRef<ST>& dx = std::nullopt) const {
+    return [&]() -> ST {
+      using namespace Dune::DerivativeDirections;
+      using namespace Dune;
+      ST energy      = 0.0;
+      const auto eps = strainFunction(par, dx);
+      for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
+        const auto EVoigt = eps.evaluate(gpIndex, on(gridElement));
+        energy += internalEnergy(EVoigt) * geo_->integrationElement(gp.position()) * gp.weight();
+      }
+      return energy;
+    };
+  }
+
 protected:
   /**
    * \brief Calculate the matrix associated with the given Requirement.
@@ -244,26 +343,25 @@ protected:
    * \param K The matrix to store the calculated result.
    */
   template <typename ScalarType>
-  void calculateMatrixImpl(
-      const Requirement& par, const MatrixAffordance& affordance, typename Traits::template MatrixType<> K,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  void calculateMatrixImpl(const Requirement& par, const MatrixAffordance& affordance,
+                           typename Traits::template MatrixType<> K,
+                           const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
+    if constexpr (hasEAS)
+      return;
     using namespace Dune::DerivativeDirections;
     using namespace Dune;
     const auto uFunction = displacementFunction(par, dx);
     const auto eps       = strainFunction(par, dx);
+    const auto kFunction = stiffnessMatrixFunction<ScalarType>(par, K, dx);
     for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
-      const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
-      const auto EVoigt       = (eps.evaluate(gpIndex, on(gridElement))).eval();
-      const auto u            = (uFunction.evaluate(gpIndex, on(gridElement))).eval();
-      const auto C            = materialTangent(EVoigt);
-
-      const auto stresses = getStress(EVoigt);
+      const auto EVoigt   = (eps.evaluate(gpIndex, on(gridElement))).eval();
+      const auto stresses = stress(EVoigt);
       for (size_t i = 0; i < numberOfNodes_; ++i) {
         const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
         for (size_t j = 0; j < numberOfNodes_; ++j) {
           const auto bopJ = eps.evaluateDerivative(gpIndex, wrt(coeff(j)), on(gridElement));
           const auto kgIJ = eps.evaluateDerivative(gpIndex, wrt(coeff(i, j)), along(stresses), on(gridElement));
-          K.template block<myDim, myDim>(i * myDim, j * myDim) += (bopI.transpose() * C * bopJ + kgIJ) * intElement;
+          kFunction(EVoigt, bopI, bopJ, kgIJ, i, j, gp);
         }
       }
     }
@@ -271,40 +369,29 @@ protected:
 
   template <typename ScalarType>
   auto calculateScalarImpl(const Requirement& par, ScalarAffordance affordance,
-                           const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx =
-                               std::nullopt) const -> ScalarType {
-    using namespace Dune::DerivativeDirections;
-    using namespace Dune;
-
-    const auto eps     = strainFunction(par, dx);
-    const auto& lambda = par.parameter();
-    ScalarType energy  = 0.0;
-
-    for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
-      const auto EVoigt         = (eps.evaluate(gpIndex, on(gridElement))).eval();
-      const auto internalEnergy = getInternalEnergy(EVoigt);
-      energy += internalEnergy * geo_->integrationElement(gp.position()) * gp.weight();
-    }
-
-    return energy;
+                           const VectorXOptRef<ScalarType>& dx = std::nullopt) const -> ScalarType {
+    if constexpr (hasEAS)
+      return ScalarType{0.0};
+    return energyFunction(par, dx)();
   }
 
   template <typename ScalarType>
-  void calculateVectorImpl(
-      const Requirement& par, VectorAffordance affordance, typename Traits::template VectorType<ScalarType> force,
-      const std::optional<std::reference_wrapper<const Eigen::VectorX<ScalarType>>>& dx = std::nullopt) const {
+  void calculateVectorImpl(const Requirement& par, VectorAffordance affordance,
+                           typename Traits::template VectorType<ScalarType> force,
+                           const VectorXOptRef<ScalarType>& dx = std::nullopt) const {
+    if constexpr (hasEAS)
+      return;
     using namespace Dune::DerivativeDirections;
     using namespace Dune;
-    const auto eps = strainFunction(par, dx);
+    const auto eps          = strainFunction(par, dx);
+    const auto fIntFunction = internalForcesFunction<ScalarType>(par, force, dx);
 
     // Internal forces
     for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
-      const double intElement = geo_->integrationElement(gp.position()) * gp.weight();
-      const auto EVoigt       = (eps.evaluate(gpIndex, on(gridElement))).eval();
-      const auto stresses     = getStress(EVoigt);
+      const auto EVoigt = (eps.evaluate(gpIndex, on(gridElement))).eval();
       for (size_t i = 0; i < numberOfNodes_; ++i) {
         const auto bopI = eps.evaluateDerivative(gpIndex, wrt(coeff(i)), on(gridElement));
-        force.template segment<myDim>(myDim * i) += bopI.transpose() * stresses * intElement;
+        fIntFunction(EVoigt, bopI, i, gp);
       }
     }
   }
diff --git a/ikarus/finiteelements/mechanics/eas/CMakeLists.txt b/ikarus/finiteelements/mechanics/strainenhancements/CMakeLists.txt
similarity index 60%
rename from ikarus/finiteelements/mechanics/eas/CMakeLists.txt
rename to ikarus/finiteelements/mechanics/strainenhancements/CMakeLists.txt
index 22cfc2324..eb56a0aba 100644
--- a/ikarus/finiteelements/mechanics/eas/CMakeLists.txt
+++ b/ikarus/finiteelements/mechanics/strainenhancements/CMakeLists.txt
@@ -1,7 +1,8 @@
 # SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
 # SPDX-License-Identifier: LGPL-3.0-or-later
-
+add_subdirectory(easfunctions)
+add_subdirectory(easvariants)
 # install headers
-install(FILES eas2d.hh eas3d.hh
-        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements/mechanics/eas
+install(FILES easvariants.hh easfunctions.hh
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements/mechanics/strainenhancements
 )
diff --git a/ikarus/finiteelements/mechanics/strainenhancements/easfunctions.hh b/ikarus/finiteelements/mechanics/strainenhancements/easfunctions.hh
new file mode 100644
index 000000000..64fb26885
--- /dev/null
+++ b/ikarus/finiteelements/mechanics/strainenhancements/easfunctions.hh
@@ -0,0 +1,13 @@
+// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+/**
+ * \file easfunctions.hh
+ * \brief Header file for various EAS functions.
+ * \ingroup eas
+ */
+
+#pragma once
+
+#include <ikarus/finiteelements/mechanics/strainenhancements/easfunctions/greenlagrangestrain.hh>
+#include <ikarus/finiteelements/mechanics/strainenhancements/easfunctions/linearstrain.hh>
diff --git a/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/CMakeLists.txt b/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/CMakeLists.txt
new file mode 100644
index 000000000..a17dfe9eb
--- /dev/null
+++ b/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/CMakeLists.txt
@@ -0,0 +1,9 @@
+# SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+# SPDX-License-Identifier: LGPL-3.0-or-later
+
+# install headers
+install(
+  FILES greenlagrangestrain.hh linearstrain.hh
+  DESTINATION
+    ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements/mechanics/strainenhancements/easfunctions
+)
diff --git a/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/greenlagrangestrain.hh b/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/greenlagrangestrain.hh
new file mode 100644
index 000000000..14753a7de
--- /dev/null
+++ b/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/greenlagrangestrain.hh
@@ -0,0 +1,152 @@
+// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+/**
+ * \file greenlagrangestrain.hh
+ * \brief Implementation of the EAS function, where the Green-Lagrange strain is enhanced.
+ *
+ * \ingroup eas
+ */
+
+#pragma once
+
+#include <dune/common/fvector.hh>
+#include <dune/localfefunctions/eigenDuneTransformations.hh>
+#include <dune/localfefunctions/expressions/greenLagrangeStrains.hh>
+
+#include <Eigen/Core>
+
+namespace Ikarus::EAS {
+/**
+ * \brief A struct computing the value, first and second derivatives of the Green-Lagrange strain tensor (in Voigt
+ * notation), where the Green-Lagrange strain itself is enhanced.
+ */
+struct GreenLagrangeStrain
+{
+  /**
+   * \brief Compute the strain vector at a given integration point or its index.
+   *
+   * \param geo The geometry object providing the transposed Jacobian.
+   * \param uFunction The function representing the displacement field.
+   * \param gpPos The position of the integration point.
+   * \param easFunction The EAS function.
+   * \param alpha The coefficients of the EAS function.
+   *
+   * \tparam GEO The type of the geometry object.
+   * \tparam EAST The type of the EAS function.
+   *
+   * \return The Green-Lagrange strain vector at the given integration point.
+   */
+  template <typename GEO, typename EAST>
+  static auto value(const GEO& geo, const auto& uFunction, const Dune::FieldVector<double, GEO::mydimension>& gpPos,
+                    const EAST& easFunction, const auto& alpha) {
+    using namespace Dune::DerivativeDirections;
+    using namespace Dune;
+    constexpr int myDim        = GEO::mydimension;
+    const auto& strainFunction = Dune::greenLagrangeStrains(uFunction);
+    const auto E               = (strainFunction.evaluate(gpPos, on(gridElement))).eval();
+    const auto Etilde          = (easFunction(gpPos) * alpha).eval();
+    return (E + Etilde).eval();
+  }
+
+  /**
+   * \brief Compute the first derivative of the Green-Lagrange strain w.r.t d or alpha for a given node and integration
+   * point.
+   *
+   * \param geo The geometry object of the finite element.
+   * \param uFunction The function representing the displacement field.
+   * \param localBasis The local basis object.
+   * \param gpPos The position of the integration point.
+   * \param gpIndex The index of the integration point.
+   * \param node The FE node index (defaults to sNaN).
+   * \param easFunction The EAS function.
+   * \param alpha The coefficients of the EAS function.
+   *
+   * \tparam wrtCoeff An integer indicating the coefficient w.r.t which the first derivative of Green-Lagrange strain is
+   * to be returned (0 = E,d; 1 = E,a).
+   * \tparam GEO The type of the geometry object.
+   * \tparam EAST The type of the EAS
+   * function.
+   *
+   * \return The first derivative of the Green-Lagrange strain w.r.t d or alpha.
+   */
+  template <int wrtCoeff, typename GEO, typename EAST>
+  static auto firstDerivative(const GEO& geo, const auto& uFunction, const auto& localBasis, const auto& gpIndex,
+                              const Dune::FieldVector<double, GEO::mydimension>& gpPos, const EAST& easFunction,
+                              const auto& alpha, const int node = sNaN) {
+    if constexpr (wrtCoeff == 0) {
+      using namespace Dune::DerivativeDirections;
+      using namespace Dune;
+      constexpr int myDim        = GEO::mydimension;
+      const auto& strainFunction = Dune::greenLagrangeStrains(uFunction);
+      const auto bopI            = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(node)), on(gridElement));
+      return bopI.eval();
+    } else if constexpr (wrtCoeff == 1) {
+      return easFunction(gpPos);
+    } else
+      static_assert(Dune::AlwaysFalse<GEO>::value,
+                    "firstDerivative can only be called with wrtCoeff as 0 and 1 indicating first derivative of the "
+                    "Green-Lagrange strain w.r.t d and "
+                    "alpha, respectively.");
+  }
+
+  /**
+   * \brief Compute the second derivative of the Green-Lagrange strain w.r.t d or alpha for a given node and integration
+   * point.
+   *
+   * \param geo The geometry object of the finite element.
+   * \param uFunction The function representing the displacement field.
+   * \param localBasis The local basis object.
+   * \param gpPos The position of the integration point.
+   * \param gpIndex The index of the integration point.
+   * \param I The FE node index I (defaults to sNaN).
+   * \param J The FE node index J (defaults to sNaN).
+   * \param S The PK2 stress (in Voigt notation).
+   * \param easFunction The EAS function.
+   * \param alpha The coefficients of the EAS function.
+   *
+   * \tparam wrtCoeff An integer indicating the coefficient w.r.t which the second derivative of Green-Lagrange strain
+   * is to be returned (0 = E,dd; 1 = E,aa; 2 = E,ad = E,da).
+   * \tparam GEO The type of the geometry object.
+   * \tparam ST
+   * The underlying scalar type.
+   * \tparam EAST The type of the EAS function.
+   *
+   * \return The second derivative of the Green-Lagrange strain w.r.t d, alpha or the mixed derivative
+   * (0 = E,dd; 1 = E,aa; 2 = E,ad = E,da).
+   */
+  template <int wrtCoeff, typename GEO, typename ST, typename EAST>
+  static auto secondDerivative(const GEO& geo, const auto& uFunction, const auto& localBasis, const auto& gpIndex,
+                               const Dune::FieldVector<double, GEO::mydimension>& gpPos,
+                               const Eigen::Vector<ST, GEO::mydimension*(GEO::mydimension + 1) / 2>& S,
+                               const EAST& easFunction, const auto& alpha, const int I = sNaN, const int J = sNaN) {
+    constexpr int myDim              = GEO::mydimension;
+    constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
+    if constexpr (wrtCoeff == 0) {
+      using namespace Dune::DerivativeDirections;
+      using namespace Dune;
+      const auto& strainFunction = Dune::greenLagrangeStrains(uFunction);
+      const auto kg = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(I, J)), along(S), on(gridElement));
+      return kg.eval(); // E_dd
+    } else if constexpr (wrtCoeff == 1) {
+      const auto kg = Eigen::Matrix<ST, enhancedStrainSize, enhancedStrainSize>::Zero().eval();
+      return kg; // E_aa
+    } else if constexpr (wrtCoeff == 2) {
+      const auto kg = Eigen::Matrix<ST, enhancedStrainSize, myDim>::Zero().eval();
+      return kg; // E_ad
+    } else
+      static_assert(
+          Dune::AlwaysFalse<GEO>::value,
+          "secondDerivative can only be called with wrtCoeff as 0, 1 and 2 indicating second derivative of the "
+          "Green-Lagrange strain w.r.t d and "
+          "alpha and the mixed derivative, respectively.");
+  }
+
+  /** \brief The name of the strain measure enhanced w.r.t EAS method. */
+  static constexpr auto name() { return std::string("Green-Lagrange Strain"); }
+
+private:
+  static constexpr int sNaN = std::numeric_limits<int>::signaling_NaN();
+};
+
+} // namespace Ikarus::EAS
diff --git a/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/linearstrain.hh b/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/linearstrain.hh
new file mode 100644
index 000000000..a9681ec7c
--- /dev/null
+++ b/ikarus/finiteelements/mechanics/strainenhancements/easfunctions/linearstrain.hh
@@ -0,0 +1,147 @@
+// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+/**
+ * \file linearstrain.hh
+ * \brief Implementation of the EAS function, where the linear strain is enhanced.
+ *
+ * \ingroup eas
+ */
+
+#pragma once
+
+#include <dune/common/fvector.hh>
+#include <dune/localfefunctions/eigenDuneTransformations.hh>
+#include <dune/localfefunctions/expressions/linearStrainsExpr.hh>
+
+#include <Eigen/Core>
+
+namespace Ikarus::EAS {
+/**
+ * \brief A struct computing the value, first and second derivatives of the linear strain tensor (in Voigt
+ * notation), where the linear strain itself is enhanced.
+ */
+struct LinearStrain
+{
+  /**
+   * \brief Compute the strain vector at a given integration point or its index.
+   *
+   * \param geo The geometry object providing the transposed Jacobian.
+   * \param uFunction The function representing the displacement field.
+   * \param gpPos The position of the integration point.
+   * \param easFunction The EAS function.
+   * \param alpha The coefficients of the EAS function.
+   *
+   * \tparam GEO The type of the geometry object.
+   * \tparam EAST The type of the EAS function.
+   *
+   * \return The linear strain vector at the given integration point.
+   */
+  template <typename GEO, typename EAST>
+  static auto value(const GEO& geo, const auto& uFunction, const Dune::FieldVector<double, GEO::mydimension>& gpPos,
+                    const EAST& easFunction, const auto& alpha) {
+    using namespace Dune::DerivativeDirections;
+    using namespace Dune;
+    constexpr int myDim        = GEO::mydimension;
+    const auto& strainFunction = Dune::linearStrains(uFunction);
+    const auto E               = (strainFunction.evaluate(gpPos, on(gridElement))).eval();
+    const auto Etilde          = (easFunction(gpPos) * alpha).eval();
+    return (E + Etilde).eval();
+  }
+
+  /**
+   * \brief Compute the first derivative of the linear strain w.r.t d or alpha for a given node and integration
+   * point.
+   *
+   * \param geo The geometry object of the finite element.
+   * \param uFunction The function representing the displacement field.
+   * \param localBasis The local basis object.
+   * \param gpPos The position of the integration point.
+   * \param gpIndex The index of the integration point.
+   * \param node The FE node index (defaults to sNaN).
+   * \param easFunction The EAS function.
+   * \param alpha The coefficients of the EAS function.
+   *
+   * \tparam wrtCoeff An integer indicating the coefficient w.r.t which the first derivative of linear strain is
+   * to be returned (0 = E,d; 1 = E,a).
+   * \tparam GEO The type of the geometry object.
+   * \tparam EAST The type of the EAS function.
+   *
+   * \return The first derivative of the linear strain w.r.t d or alpha.
+   */
+  template <int wrtCoeff, typename GEO, typename EAST>
+  static auto firstDerivative(const GEO& geo, const auto& uFunction, const auto& localBasis, const auto& gpIndex,
+                              const Dune::FieldVector<double, GEO::mydimension>& gpPos, const EAST& easFunction,
+                              const auto& alpha, const int node = sNaN) {
+    if constexpr (wrtCoeff == 0) {
+      using namespace Dune::DerivativeDirections;
+      using namespace Dune;
+      constexpr int myDim        = GEO::mydimension;
+      const auto& strainFunction = Dune::linearStrains(uFunction);
+      const auto bopI            = strainFunction.evaluateDerivative(gpIndex, wrt(coeff(node)), on(gridElement));
+      return bopI.eval();
+    } else if constexpr (wrtCoeff == 1) {
+      return easFunction(gpPos);
+    } else
+      static_assert(Dune::AlwaysFalse<GEO>::value,
+                    "firstDerivative can only be called with wrtCoeff as 0 and 1 indicating first derivative of the "
+                    "linear strain w.r.t d and "
+                    "alpha, respectively.");
+  }
+
+  /**
+   * \brief Compute the second derivative of the linear strain w.r.t d or alpha for a given node and integration
+   * point.
+   *
+   * \param geo The geometry object of the finite element.
+   * \param uFunction The function representing the displacement field.
+   * \param localBasis The local basis object.
+   * \param gpPos The position of the integration point.
+   * \param gpIndex The index of the integration point.
+   * \param I The FE node index I (defaults to sNaN).
+   * \param J The FE node index J (defaults to sNaN).
+   * \param S The PK2 stress (in Voigt notation).
+   * \param easFunction The EAS function.
+   * \param alpha The coefficients of the EAS function.
+   *
+   * \tparam wrtCoeff An integer indicating the coefficient w.r.t which the second derivative of linear strain
+   * is to be returned (0 = E,dd; 1 = E,aa; 2 = E,ad = E,da).
+   * \tparam GEO The type of the geometry object.
+   * \tparam ST The underlying scalar type.
+   * \tparam EAST The type of the EAS function.
+   *
+   * \return A tuple containing second derivative of the linear strain w.r.t d, alpha and the mixed derivative
+   * (0 = E,dd; 1 = E,aa; 2 = E,ad = E,da).
+   */
+  template <int wrtCoeff, typename GEO, typename ST, typename EAST>
+  static auto secondDerivative(const GEO& geo, const auto& uFunction, const auto& localBasis, const auto& gpIndex,
+                               const Dune::FieldVector<double, GEO::mydimension>& gpPos,
+                               const Eigen::Vector<ST, GEO::mydimension*(GEO::mydimension + 1) / 2>& S,
+                               const EAST& easFunction, const auto& alpha, const int I = sNaN, const int J = sNaN) {
+    constexpr int myDim              = GEO::mydimension;
+    constexpr int enhancedStrainSize = EAST::enhancedStrainSize;
+    if constexpr (wrtCoeff == 0) {
+      const auto kg = Eigen::Matrix<ST, myDim, myDim>::Zero().eval();
+      return kg; // E_dd
+    } else if constexpr (wrtCoeff == 1) {
+      const auto kg = Eigen::Matrix<ST, enhancedStrainSize, enhancedStrainSize>::Zero().eval();
+      return kg; // E_aa
+    } else if constexpr (wrtCoeff == 2) {
+      const auto kg = Eigen::Matrix<ST, enhancedStrainSize, myDim>::Zero().eval();
+      return kg; // E_ad
+    } else
+      static_assert(
+          Dune::AlwaysFalse<GEO>::value,
+          "secondDerivative can only be called with wrtCoeff as 0, 1 and 2 indicating second derivative of the "
+          "linear strain w.r.t d and "
+          "alpha and the mixed derivative, respectively.");
+  }
+
+  /** \brief The name of the strain measure enhanced w.r.t EAS method. */
+  static constexpr auto name() { return std::string("Linear Strain"); }
+
+private:
+  static constexpr int sNaN = std::numeric_limits<int>::signaling_NaN();
+};
+
+} // namespace Ikarus::EAS
diff --git a/ikarus/finiteelements/mechanics/strainenhancements/easvariants.hh b/ikarus/finiteelements/mechanics/strainenhancements/easvariants.hh
new file mode 100644
index 000000000..d4abca030
--- /dev/null
+++ b/ikarus/finiteelements/mechanics/strainenhancements/easvariants.hh
@@ -0,0 +1,128 @@
+// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+/**
+ * \file easvariants.hh
+ * \brief Definition of the EAS variants.
+ * \ingroup eas
+ */
+
+#pragma once
+
+#include <ikarus/finiteelements/mechanics/strainenhancements/easfunctions.hh>
+#include <ikarus/finiteelements/mechanics/strainenhancements/easvariants/linearandglstrains.hh>
+#include <ikarus/utils/tensorutils.hh>
+
+namespace Ikarus::EAS {
+namespace Impl {
+  template <typename GEO, int myDim, typename ESS>
+  struct EASType;
+  template <typename GEO>
+  struct EASType<GEO, 3, EAS::LinearStrain>
+  {
+    using type = std::variant<E0<GEO>, E9<GEO>, E21<GEO>>;
+  };
+  template <typename GEO>
+  struct EASType<GEO, 2, EAS::LinearStrain>
+  {
+    using type = std::variant<E0<GEO>, E4<GEO>, E5<GEO>, E7<GEO>, E11<GEO>>;
+  };
+  template <typename GEO, int myDim>
+  struct EASType<GEO, myDim, EAS::GreenLagrangeStrain>
+  {
+    using type = EASType<GEO, myDim, EAS::LinearStrain>::type;
+  };
+} // namespace Impl
+
+/**
+ * \brief Wrapper around the EAS variant, contains helper functions
+ * \tparam GEO the Geometry type of the udderlying grid element
+ */
+template <typename ES, typename GEO>
+struct EASVariant
+{
+  static constexpr int myDim = GEO::mydimension;
+  using Variant              = Impl::EASType<GEO, myDim, ES>::type;
+  static_assert((myDim == 2) or (myDim == 3), "EAS variants are only available for 2D and 3D elements.");
+
+  /**
+   * \brief Executes a function F on the variant, if `enhancedStrainSize` is not zero
+   * \tparam F the Type of the function F
+   * \param f the function, which takes the eas element as an argument
+   */
+  template <typename F>
+  void operator()(F&& f) const {
+    std::visit([&]<typename EAST>(const EAST& easFunction) { f(easFunction); }, var_);
+  }
+
+  /**
+   * \brief A helper function to get the number of EAS parameters.
+   * \return Number of EAS parameters.
+   */
+  auto numberOfEASParameters() const {
+    return std::visit([]<typename EAST>(const EAST&) -> int { return EAST::enhancedStrainSize; }, var_);
+  }
+
+  /**
+   * \brief A helper function to identify if the formulation is purely displacement-based, i.e., number of EAS
+   * parameters is zero.
+   * \return A bool indicating if the formulation is displacement-based or not.
+   */
+  bool isDisplacmentBased() const { return numberOfEASParameters() == 0; }
+
+  void setEASType(int numberOfEASParameters) {
+    numberOfEASParameters_ = numberOfEASParameters;
+    if (geometry_)
+      createEASType();
+  }
+  void bind(const GEO& geometry) {
+    geometry_ = std::make_optional<GEO>(geometry);
+    createEASType();
+  }
+
+private:
+  void createEASType() {
+    const std::string& errorMessage = "The given EAS parameters are not available for enhancing " + ES::name() +
+                                      " strain measure for the " + std::to_string(myDim) + "D case.";
+
+    if constexpr (std::same_as<ES, EAS::LinearStrain> or std::same_as<ES, EAS::GreenLagrangeStrain>) {
+      if (numberOfEASParameters_ == 0) {
+        var_ = E0(geometry_.value());
+        return;
+      }
+      if constexpr (myDim == 2) {
+        switch (numberOfEASParameters_) {
+          case 4:
+            var_ = E4(geometry_.value());
+            break;
+          case 5:
+            var_ = E5(geometry_.value());
+            break;
+          case 7:
+            var_ = E7(geometry_.value());
+            break;
+          case 11:
+            var_ = E11(geometry_.value());
+            break;
+          default:
+            DUNE_THROW(Dune::NotImplemented, errorMessage);
+        }
+      } else if constexpr (myDim == 3) {
+        switch (numberOfEASParameters_) {
+          case 9:
+            var_ = E9(geometry_.value());
+            break;
+          case 21:
+            var_ = E21(geometry_.value());
+            break;
+          default:
+            DUNE_THROW(Dune::NotImplemented, errorMessage);
+        }
+      }
+    }
+  }
+  std::optional<GEO> geometry_;
+  Variant var_;
+  int numberOfEASParameters_;
+};
+} // namespace Ikarus::EAS
diff --git a/ikarus/finiteelements/mechanics/strainenhancements/easvariants/CMakeLists.txt b/ikarus/finiteelements/mechanics/strainenhancements/easvariants/CMakeLists.txt
new file mode 100644
index 000000000..94b5831bc
--- /dev/null
+++ b/ikarus/finiteelements/mechanics/strainenhancements/easvariants/CMakeLists.txt
@@ -0,0 +1,9 @@
+# SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+# SPDX-License-Identifier: LGPL-3.0-or-later
+
+# install headers
+install(
+  FILES linearandglstrains.hh helperfunctions.hh
+  DESTINATION
+    ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements/mechanics/strainenhancements/easvariants
+)
diff --git a/ikarus/finiteelements/mechanics/strainenhancements/easvariants/helperfunctions.hh b/ikarus/finiteelements/mechanics/strainenhancements/easvariants/helperfunctions.hh
new file mode 100644
index 000000000..c30e83f8a
--- /dev/null
+++ b/ikarus/finiteelements/mechanics/strainenhancements/easvariants/helperfunctions.hh
@@ -0,0 +1,26 @@
+// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+/**
+ * \file helperfunctions.hh
+ * \brief Definition of the certain helper functions used to transform strain measures.
+ *
+ * \ingroup eas
+ */
+
+#pragma once
+
+#include <dune/localfefunctions/eigenDuneTransformations.hh>
+
+#include <ikarus/utils/tensorutils.hh>
+
+namespace Ikarus::EAS::Impl {
+template <typename GEO>
+auto transformationMatrixAtCenterWithDetJ(const GEO& geometry) {
+  const auto& referenceElement = Dune::ReferenceElements<double, GEO::mydimension>::general(geometry.type());
+  const auto quadPos0          = referenceElement.position(0, 0);
+  const auto detJ0             = geometry.integrationElement(quadPos0);
+
+  return (transformationMatrix(geometry, quadPos0) * detJ0).eval();
+}
+} // namespace Ikarus::EAS::Impl
\ No newline at end of file
diff --git a/ikarus/finiteelements/mechanics/strainenhancements/easvariants/linearandglstrains.hh b/ikarus/finiteelements/mechanics/strainenhancements/easvariants/linearandglstrains.hh
new file mode 100644
index 000000000..bfe3f7c37
--- /dev/null
+++ b/ikarus/finiteelements/mechanics/strainenhancements/easvariants/linearandglstrains.hh
@@ -0,0 +1,330 @@
+// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+/**
+ * \file linearandglstrains.hh
+ * \brief Definition of the types of EAS formulations for 2D and 3D linear solid elements, where linear and
+ * Green-Lagrange strains are enhanced.
+ *
+ * \ingroup eas
+ */
+
+#pragma once
+
+#include <ikarus/finiteelements/mechanics/strainenhancements/easvariants/helperfunctions.hh>
+#include <ikarus/utils/tensorutils.hh>
+
+namespace Ikarus::EAS {
+
+/**
+ * \brief Interface for displacement-based EAS elements, where linear or Green-Lagrange strains are enhanced.
+ */
+template <typename GEO, int ess>
+struct EX
+{
+  static constexpr int myDim              = GEO::mydimension;
+  static constexpr int strainSize         = myDim * (myDim + 1) / 2;
+  static constexpr int enhancedStrainSize = ess;
+  using AnsatzType                        = Eigen::Matrix<double, strainSize, enhancedStrainSize>;
+  using DType                             = Eigen::Matrix<double, enhancedStrainSize, enhancedStrainSize>;
+
+  EX() = default;
+  explicit EX(const GEO& geometry)
+      : geometry_{std::make_optional<GEO>(geometry)},
+        T0InverseTransformed_{Impl::transformationMatrixAtCenterWithDetJ(geometry).inverse()} {}
+
+protected:
+  std::optional<GEO> geometry_;
+  Eigen::Matrix<double, strainSize, strainSize> T0InverseTransformed_;
+};
+
+/**
+ * \brief Dummy struct for displacement-based EAS elements, i.e. 0 enhanced modes
+ */
+template <typename GEO>
+struct E0 : EX<GEO, 0>
+{
+  using Base                              = EX<GEO, 0>;
+  static constexpr int myDim              = Base::myDim;
+  static constexpr int strainSize         = Base::strainSize;
+  static constexpr int enhancedStrainSize = Base::enhancedStrainSize;
+  using AnsatzType                        = typename Base::AnsatzType;
+  using DType                             = typename Base::DType;
+
+  E0() = default;
+  explicit E0(const GEO& geo)
+      : Base(geo) {}
+
+  // returns an Eigen zero expression for optimization purposes
+  auto operator()(const Dune::FieldVector<double, myDim>& /*quadPos*/) const { return AnsatzType::Zero(); }
+};
+
+/**
+ * \brief E4 structure for EAS with linear strains and 4 enhanced modes.
+ *
+ * \details The E4 structure represents an implementation of EAS for a
+ * specific case where linear strains are considered, and the method includes
+ * enhancement with 4 additional modes to improve the accuracy of finite element solutions.
+ *
+ * \tparam GEO The geometry type.
+ */
+template <typename GEO>
+struct E4 : EX<GEO, 4>
+{
+  using Base                              = EX<GEO, 4>;
+  static constexpr int myDim              = Base::myDim;
+  static constexpr int strainSize         = Base::strainSize;
+  static constexpr int enhancedStrainSize = Base::enhancedStrainSize;
+  using AnsatzType                        = typename Base::AnsatzType;
+  using DType                             = typename Base::DType;
+
+  E4() = default;
+  explicit E4(const GEO& geo)
+      : Base(geo) {}
+
+  AnsatzType operator()(const Dune::FieldVector<double, myDim>& quadPos) const {
+    AnsatzType M;
+    M.setZero(strainSize, enhancedStrainSize);
+    const double xi   = quadPos[0];
+    const double eta  = quadPos[1];
+    M(0, 0)           = 2 * xi - 1.0;
+    M(1, 1)           = 2 * eta - 1.0;
+    M(2, 2)           = 2 * xi - 1.0;
+    M(2, 3)           = 2 * eta - 1.0;
+    const double detJ = this->geometry_->integrationElement(quadPos);
+    M                 = this->T0InverseTransformed_ / detJ * M;
+    return M;
+  }
+};
+
+/**
+ * \brief Structure representing EAS for Q1 with 5 enhanced strains.
+ *
+ * This structure defines the EAS for Q1 elements with 5 enhanced strains.
+ *
+ * \tparam GEO The geometry type.
+ */
+template <typename GEO>
+struct E5 : EX<GEO, 5>
+{
+  using Base                              = EX<GEO, 5>;
+  static constexpr int myDim              = Base::myDim;
+  static constexpr int strainSize         = Base::strainSize;
+  static constexpr int enhancedStrainSize = Base::enhancedStrainSize;
+  using AnsatzType                        = typename Base::AnsatzType;
+  using DType                             = typename Base::DType;
+
+  E5() = default;
+  explicit E5(const GEO& geo)
+      : Base(geo) {}
+
+  AnsatzType operator()(const Dune::FieldVector<double, myDim>& quadPos) const {
+    AnsatzType M;
+    M.setZero();
+    const double xi   = quadPos[0];
+    const double eta  = quadPos[1];
+    M(0, 0)           = 2 * xi - 1.0;
+    M(1, 1)           = 2 * eta - 1.0;
+    M(2, 2)           = 2 * xi - 1.0;
+    M(2, 3)           = 2 * eta - 1.0;
+    M(2, 4)           = (2 * xi - 1.0) * (2 * eta - 1.0);
+    const double detJ = this->geometry_->integrationElement(quadPos);
+    M                 = this->T0InverseTransformed_ / detJ * M;
+    return M;
+  }
+};
+
+/**
+ * \brief Structure representing EAS for Q1 with 7 enhanced strains.
+ *
+ * This structure defines the EAS for Q1 elements with 7 enhanced strains.
+ *
+ * \tparam GEO The geometry type.
+ */
+template <typename GEO>
+struct E7 : EX<GEO, 7>
+{
+  using Base                              = EX<GEO, 7>;
+  static constexpr int myDim              = Base::myDim;
+  static constexpr int strainSize         = Base::strainSize;
+  static constexpr int enhancedStrainSize = Base::enhancedStrainSize;
+  using AnsatzType                        = typename Base::AnsatzType;
+  using DType                             = typename Base::DType;
+
+  E7() = default;
+  explicit E7(const GEO& geo)
+      : Base(geo) {}
+
+  AnsatzType operator()(const Dune::FieldVector<double, myDim>& quadPos) const {
+    AnsatzType M;
+    M.setZero();
+    const double xi   = quadPos[0];
+    const double eta  = quadPos[1];
+    M(0, 0)           = 2 * xi - 1.0;
+    M(1, 1)           = 2 * eta - 1.0;
+    M(2, 2)           = 2 * xi - 1.0;
+    M(2, 3)           = 2 * eta - 1.0;
+    M(0, 4)           = (2 * xi - 1.0) * (2 * eta - 1.0);
+    M(1, 5)           = (2 * xi - 1.0) * (2 * eta - 1.0);
+    M(2, 6)           = (2 * xi - 1.0) * (2 * eta - 1.0);
+    const double detJ = this->geometry_->integrationElement(quadPos);
+    M                 = this->T0InverseTransformed_ / detJ * M;
+    return M;
+  }
+};
+
+/**
+ * \brief Structure representing EAS for Q2 with 11 enhanced strains.
+ *
+ * \details See \cite bischoffShearDeformableShell1997b for details.
+ *
+ * This structure defines the EAS for Q2 elements with 11 enhanced strains.
+ *
+ * \tparam GEO The geometry type.
+ */
+template <typename GEO>
+struct E11 : EX<GEO, 11>
+{
+  using Base                              = EX<GEO, 11>;
+  static constexpr int myDim              = Base::myDim;
+  static constexpr int strainSize         = Base::strainSize;
+  static constexpr int enhancedStrainSize = Base::enhancedStrainSize;
+  using AnsatzType                        = typename Base::AnsatzType;
+  using DType                             = typename Base::DType;
+
+  E11() = default;
+  explicit E11(const GEO& geo)
+      : Base(geo) {}
+
+  AnsatzType operator()(const Dune::FieldVector<double, myDim>& quadPos) const {
+    AnsatzType M;
+    M.setZero();
+    const double xi    = quadPos[0];
+    const double eta   = quadPos[1];
+    const double xi_t  = 2 * xi - 1;
+    const double eta_t = 2 * eta - 1;
+
+    const double xi2    = xi_t * xi_t;
+    const double eta2   = eta_t * eta_t;
+    const double xi_eta = xi_t * eta_t;
+
+    M(0, 0) = 1 - 3 * xi2;
+    M(0, 1) = eta_t - 3 * xi2 * eta_t;
+    M(0, 2) = eta2 - 3 * xi2 * eta2;
+
+    M(1, 3) = 1 - 3 * eta2;
+    M(1, 4) = xi_t - 3 * xi_t * eta2;
+    M(1, 5) = xi2 - 3 * xi2 * eta2;
+
+    M(2, 6)  = 1 - 3 * xi2;
+    M(2, 7)  = 1 - 3 * eta2;
+    M(2, 8)  = eta_t - 3 * xi2 * eta_t;
+    M(2, 9)  = xi_t - 3 * xi_t * eta2;
+    M(2, 10) = 1 - 3 * (xi2 + eta2) + 9 * xi2 * eta2;
+
+    const double detJ = this->geometry_->integrationElement(quadPos);
+    M                 = this->T0InverseTransformed_ / detJ * M;
+    return M;
+  }
+};
+
+/**
+ * \brief Structure representing EAS for H1 with 9 enhanced strains.
+ *
+ * This structure defines the EAS for H1 elements with 9 enhanced strains.
+ *
+ * \tparam GEO The geometry type.
+ */
+template <typename GEO>
+struct E9 : EX<GEO, 9>
+{
+  using Base                              = EX<GEO, 9>;
+  static constexpr int myDim              = Base::myDim;
+  static constexpr int strainSize         = Base::strainSize;
+  static constexpr int enhancedStrainSize = Base::enhancedStrainSize;
+  using AnsatzType                        = typename Base::AnsatzType;
+  using DType                             = typename Base::DType;
+
+  E9() = default;
+  explicit E9(const GEO& geo)
+      : Base(geo) {}
+
+  AnsatzType operator()(const Dune::FieldVector<double, myDim>& quadPos) const {
+    AnsatzType M;
+    M.setZero();
+    const double xi   = quadPos[0];
+    const double eta  = quadPos[1];
+    const double zeta = quadPos[2];
+    M(0, 0)           = 2 * xi - 1.0;
+    M(1, 1)           = 2 * eta - 1.0;
+    M(2, 2)           = 2 * zeta - 1.0;
+    M(3, 3)           = 2 * xi - 1.0;
+    M(3, 4)           = 2 * eta - 1.0;
+    M(4, 5)           = 2 * xi - 1.0;
+    M(4, 6)           = 2 * zeta - 1.0;
+    M(5, 7)           = 2 * eta - 1.0;
+    M(5, 8)           = 2 * zeta - 1.0;
+    const double detJ = this->geometry_->integrationElement(quadPos);
+    M                 = this->T0InverseTransformed_ / detJ * M;
+    return M;
+  }
+};
+
+/**
+ * \brief Structure representing EAS for H1 with 21 enhanced strains.
+ *
+ * This structure defines the EAS for H1 elements with 21 enhanced strains.
+ *
+ * \tparam GEO The geometry type.
+ */
+template <typename GEO>
+struct E21 : EX<GEO, 21>
+{
+  using Base                              = EX<GEO, 21>;
+  static constexpr int myDim              = Base::myDim;
+  static constexpr int strainSize         = Base::strainSize;
+  static constexpr int enhancedStrainSize = Base::enhancedStrainSize;
+  using AnsatzType                        = typename Base::AnsatzType;
+  using DType                             = typename Base::DType;
+
+  E21() = default;
+  explicit E21(const GEO& geo)
+      : Base(geo) {}
+
+  AnsatzType operator()(const Dune::FieldVector<double, myDim>& quadPos) const {
+    AnsatzType M;
+    M.setZero();
+    const double xi   = quadPos[0];
+    const double eta  = quadPos[1];
+    const double zeta = quadPos[2];
+    M(0, 0)           = 2 * xi - 1.0;
+    M(1, 1)           = 2 * eta - 1.0;
+    M(2, 2)           = 2 * zeta - 1.0;
+    M(3, 3)           = 2 * xi - 1.0;
+    M(3, 4)           = 2 * eta - 1.0;
+    M(4, 5)           = 2 * xi - 1.0;
+    M(4, 6)           = 2 * zeta - 1.0;
+    M(5, 7)           = 2 * eta - 1.0;
+    M(5, 8)           = 2 * zeta - 1.0;
+
+    M(3, 9)  = (2 * xi - 1.0) * (2 * zeta - 1.0);
+    M(3, 10) = (2 * eta - 1.0) * (2 * zeta - 1.0);
+    M(4, 11) = (2 * xi - 1.0) * (2 * eta - 1.0);
+    M(4, 12) = (2 * eta - 1.0) * (2 * zeta - 1.0);
+    M(5, 13) = (2 * xi - 1.0) * (2 * eta - 1.0);
+    M(5, 14) = (2 * xi - 1.0) * (2 * zeta - 1.0);
+
+    M(0, 15) = (2 * xi - 1.0) * (2 * eta - 1.0);
+    M(0, 16) = (2 * xi - 1.0) * (2 * zeta - 1.0);
+    M(1, 17) = (2 * xi - 1.0) * (2 * eta - 1.0);
+    M(1, 18) = (2 * eta - 1.0) * (2 * zeta - 1.0);
+    M(2, 19) = (2 * xi - 1.0) * (2 * zeta - 1.0);
+    M(2, 20) = (2 * eta - 1.0) * (2 * zeta - 1.0);
+
+    const double detJ = this->geometry_->integrationElement(quadPos);
+    M                 = this->T0InverseTransformed_ / detJ * M;
+    return M;
+  }
+};
+
+} // namespace Ikarus::EAS
diff --git a/ikarus/finiteelements/mixin.hh b/ikarus/finiteelements/mixin.hh
index 670a0b1b0..923e763bc 100644
--- a/ikarus/finiteelements/mixin.hh
+++ b/ikarus/finiteelements/mixin.hh
@@ -94,6 +94,9 @@ public:
   using LocalView               = typename Traits::LocalView;
   static constexpr int worldDim = Traits::worlddim;
 
+  template <typename ES>
+  static constexpr bool hasEAS = hasSkill<EnhancedAssumedStrainsPre<ES>::template Skill>();
+
   /**
    * \brief Create a Requirement object.
    *
@@ -248,11 +251,34 @@ public:
      ...);
   }
 
+private:
+  template <typename Sk>
+  auto invokeUpdateState(const Requirement& par,
+                         const std::remove_reference_t<typename Traits::template VectorType<>>& correction) {
+    if constexpr (requires { Sk::updateStateImpl(par, correction); })
+      Sk::updateStateImpl(par, correction);
+  }
+
+public:
+  /**
+   * \brief  Call all updateStateImpl functions if the skill implements it.
+   *
+   * \details Update the state variables related to a particular skill.
+   *
+   * \param req The Requirement object specifying the requirements for the update itself.
+   * \param force A correction vector (for example, the displacement increment) based on which the state variables are
+   * to be updated.
+   */
+  void updateState(const Requirement& par,
+                   const std::remove_reference_t<typename Traits::template VectorType<>>& correction) {
+    (invokeUpdateState<Skills<PreFE, typename PreFE::template FE<Skills...>>>(par, correction), ...);
+  }
+
 private:
   template <typename Sk, typename BC, typename MT>
   auto invokeSubscribeTo(BC& bc) {
-    // For Clang-16: we need the this-> otherwise the code in the if clause will never be called. For Gcc-12.2: with the
-    // this-> it throws a compiler error in certain cases
+    // For Clang-16: we need the this-> otherwise the code in the if clause will never be called. For Gcc-12.2: with
+    // the this-> it throws a compiler error in certain cases
 #if defined(__clang__)
     if constexpr (requires { this->Sk::template subscribeToImpl<MT>(bc); }) {
 #else
diff --git a/ikarus/python/finiteelements/fe.hh b/ikarus/python/finiteelements/fe.hh
index 248e91d19..9a7365761 100644
--- a/ikarus/python/finiteelements/fe.hh
+++ b/ikarus/python/finiteelements/fe.hh
@@ -100,6 +100,8 @@ void registerFE(pybind11::handle scope, pybind11::class_<FE, options...> cls) {
           pybind11::keep_alive<1, 2>());
 
   cls.def("bind", [](FE& self, const GridElement& e) { self.bind(e); });
+  cls.def("updateState",
+          [](FE& self, const Requirement& req, Eigen::Ref<Eigen::VectorXd> dx) { self.updateState(req, dx); });
   cls.def("calculateScalar", [](FE& self, const Requirement& req, Ikarus::ScalarAffordance affordance) {
     return calculateScalar(self, req, affordance);
   });
diff --git a/ikarus/python/finiteelements/registerpreelement.hh b/ikarus/python/finiteelements/registerpreelement.hh
index 50e675a27..1e662475e 100644
--- a/ikarus/python/finiteelements/registerpreelement.hh
+++ b/ikarus/python/finiteelements/registerpreelement.hh
@@ -3,8 +3,11 @@
 
 #pragma once
 
+#include <dune/common/fvector.hh>
 #include <dune/python/pybind11/pybind11.h>
 
+#include <Eigen/Core>
+
 namespace Ikarus::Python {
 
 /**
diff --git a/ikarus/python/test/helperfunctions.py b/ikarus/python/test/helperfunctions.py
new file mode 100644
index 000000000..013c37753
--- /dev/null
+++ b/ikarus/python/test/helperfunctions.py
@@ -0,0 +1,67 @@
+import numpy as np
+import scipy as sp
+from scipy.optimize import minimize
+
+
+def updateAllElements(fes, req, dx):
+    for fe in fes:
+        fe.updateState(req, dx)
+
+
+def solveWithSciPyMinimize(energyFun, x0, jacFun=None, hessFun=None, callBackFun=None):
+    if jacFun == None:
+        return minimize(
+            energyFun,
+            x0=x0,
+            options={"disp": True, "gtol": 1e-3},
+            tol=1e-14,
+            callback=callBackFun,
+        )
+    elif hessFun == None:
+        return minimize(
+            energyFun,
+            x0=x0,
+            jac=jacFun,
+            options={"disp": True, "gtol": 1e-3},
+            tol=1e-14,
+            callback=callBackFun,
+        )
+    else:
+        return minimize(
+            energyFun,
+            method="trust-constr",
+            x0=x0,
+            jac=jacFun,
+            hess=hessFun,
+            options={"disp": True},
+            callback=callBackFun,
+        )
+
+
+# Using derivative free version root, as the others do not accept a callback
+def solveWithSciPyRoot(fun, x0, callBackFun=None):
+    if callBackFun != None:
+        return sp.optimize.root(
+            fun, x0=x0, tol=1e-10, method="krylov", callback=callBackFun
+        )
+
+
+def solveWithNewtonRaphson(gradFun, hessFun, assembler, fes, req):
+    maxiter = 100
+    abs_tolerance = 1e-14
+    d = np.zeros(assembler.reducedSize())
+    for k in range(maxiter):
+        R = gradFun(d)
+        K = hessFun(d)
+        r_norm = np.linalg.norm(R)
+
+        deltad = sp.linalg.solve(K, R)
+
+        updateAllElements(fes, req, assembler.createFullVector(deltad))
+        d -= deltad
+
+        if r_norm < abs_tolerance or k > maxiter:
+            break
+
+    success = r_norm < abs_tolerance
+    return success, d
diff --git a/ikarus/python/test/nonlinearelastictest.py b/ikarus/python/test/nonlinearelastictest.py
index 8ac7832d8..aeb64f7d8 100644
--- a/ikarus/python/test/nonlinearelastictest.py
+++ b/ikarus/python/test/nonlinearelastictest.py
@@ -6,15 +6,22 @@
 debug_info.setDebugFlags()
 
 import ikarus as iks
-from ikarus import finite_elements, utils, assembler
+from ikarus import finite_elements, utils, assembler, materials
 import numpy as np
 import scipy as sp
 from scipy.optimize import minimize
+from helperfunctions import (
+    updateAllElements,
+    solveWithSciPyMinimize,
+    solveWithSciPyRoot,
+    solveWithNewtonRaphson,
+)
 
 import dune.grid
 import dune.functions
 
-if __name__ == "__main__":
+
+def nonlinElasticTest(easBool):
     lowerLeft = []
     upperRight = []
     elements = []
@@ -59,12 +66,20 @@ def loadTopEdgePredicate(x):
     svkPS = svk.asPlaneStress()
 
     nonLinElastic = iks.finite_elements.nonLinearElastic(svkPS)
+    easF = iks.finite_elements.eas(4, "GreenLagrangeStrain")
 
     fes = []
     for e in grid.elements:
-        fes.append(
-            iks.finite_elements.makeFE(basisLagrange1, nonLinElastic, vLoad, nBLoad)
-        )
+        if easBool:
+            fes.append(
+                iks.finite_elements.makeFE(
+                    basisLagrange1, nonLinElastic, easF, vLoad, nBLoad
+                )
+            )
+        else:
+            fes.append(
+                iks.finite_elements.makeFE(basisLagrange1, nonLinElastic, vLoad, nBLoad)
+            )
         fes[-1].bind(e)
 
     dirichletValues = iks.dirichletValues(flatBasis)
@@ -84,61 +99,97 @@ def fixLeftHandEdge(vec, localIndex, localView, intersection):
     dirichletValues.fixBoundaryDOFs(fixAnotherVertex)
     dirichletValues.fixBoundaryDOFs(fixLeftHandEdge)
 
-    assembler = iks.assembler.sparseFlatAssembler(fes, dirichletValues)
+    sparseAssembler = iks.assembler.sparseFlatAssembler(fes, dirichletValues)
 
-    dRed = np.zeros(assembler.reducedSize())
+    dRed = np.zeros(sparseAssembler.reducedSize())
 
     lambdaLoad = iks.Scalar(3.0)
 
     feReq = fes[0].createRequirement()
+    feReq.insertGlobalSolution(np.zeros(sparseAssembler.size))
 
     def energy(dRedInput):
         feReq = fes[0].createRequirement()
         feReq.insertParameter(lambdaLoad)
-        dBig = assembler.createFullVector(dRedInput)
+        dBig = sparseAssembler.createFullVector(dRedInput)
         feReq.insertGlobalSolution(dBig)
         feReq.globalSolution()
-        return assembler.scalar(feReq, iks.ScalarAffordance.mechanicalPotentialEnergy)
+        return sparseAssembler.scalar(
+            feReq, iks.ScalarAffordance.mechanicalPotentialEnergy
+        )
 
     def gradient(dRedInput):
         feReq = fes[0].createRequirement()
         feReq.insertParameter(lambdaLoad)
-        dBig = assembler.createFullVector(dRedInput)
+        dBig = sparseAssembler.createFullVector(dRedInput)
         feReq.insertGlobalSolution(dBig)
-        return assembler.vector(
+        return sparseAssembler.vector(
             feReq, iks.VectorAffordance.forces, iks.DBCOption.Reduced
         )
 
     def hess(dRedInput):
         feReq = fes[0].createRequirement()
         feReq.insertParameter(lambdaLoad)
-        dBig = assembler.createFullVector(dRedInput)
+        dBig = sparseAssembler.createFullVector(dRedInput)
         feReq.insertGlobalSolution(dBig)
-        return assembler.matrix(
+        return sparseAssembler.matrix(
             feReq, iks.MatrixAffordance.stiffness, iks.DBCOption.Reduced
         ).todense()  # this is slow, but for this test we don't care
 
-    resultd = minimize(energy, x0=dRed, options={"disp": True}, tol=1e-14)
-    resultd2 = minimize(
-        energy, x0=dRed, jac=gradient, options={"disp": True}, tol=1e-14
-    )
-    resultd3 = minimize(
-        energy,
-        method="trust-constr",
-        x0=dRed,
-        jac=gradient,
-        hess=hess,
-        options={"disp": True},
+    def updateElements(intermediate_result: sp.optimize.OptimizeResult):
+        dx = (
+            sparseAssembler.createFullVector(intermediate_result.x)
+            - feReq.globalSolution()
+        )
+        updateAllElements(fes, feReq, dx)
+
+    def updateElementsForRoot(x, _):
+        dx = sparseAssembler.createFullVector(x) - feReq.globalSolution()
+        updateAllElements(fes, feReq, dx)
+
+    # The callback/update function has no effect for non-eas Elements (for now)
+    if not easBool:
+        resultd1 = solveWithSciPyMinimize(energy, dRed, callBackFun=updateElements)
+        resultd2 = solveWithSciPyMinimize(
+            energy, dRed, jacFun=gradient, callBackFun=updateElements
+        )
+        resultd3 = solveWithSciPyMinimize(
+            energy, dRed, jacFun=gradient, hessFun=hess, callBackFun=updateElements
+        )
+        resultd4 = solveWithSciPyRoot(gradient, dRed, callBackFun=updateElementsForRoot)
+
+        assert resultd1.success
+        assert resultd2.success
+        assert resultd3.success
+        assert resultd4.success
+
+        assert np.allclose(resultd1.x, resultd2.x, atol=1e-6)
+        assert np.allclose(resultd3.x, resultd4.x)
+        assert np.all(abs(resultd3.grad) < 1e-8)
+        assert np.all(abs(resultd4.fun) < 1e-8)
+
+        feReq.insertGlobalSolution(sparseAssembler.createFullVector(resultd4.x))
+        fes[0].calculateAt(feReq, np.array([0.5, 0.5]), "PK2Stress")
+
+    else:
+        resultd4 = solveWithSciPyRoot(gradient, dRed, callBackFun=updateElementsForRoot)
+        assert np.all(abs(resultd4.fun) < 1e-8)
+        assert resultd4.success
+
+        feReq.insertGlobalSolution(sparseAssembler.createFullVector(resultd4.x))
+
+    [successNR, resultNR] = solveWithNewtonRaphson(
+        gradient, hess, sparseAssembler, fes, feReq
     )
-    resultd4 = sp.optimize.root(gradient, jac=hess, x0=dRed, tol=1e-10)
+    assert successNR
+    assert np.allclose(resultd4.x, resultNR)
 
-    assert np.allclose(resultd.x, resultd2.x, atol=1e-6)
-    assert np.allclose(resultd3.x, resultd4.x)
-    assert np.all(abs(resultd3.grad) < 1e-8)
-    assert np.all(abs(resultd4.fun) < 1e-8)
+    return resultNR
 
-    feReq = fes[0].createRequirement()
-    fullD = assembler.createFullVector(resultd2.x)
-    feReq.insertGlobalSolution(fullD)
 
-    res1 = fes[0].calculateAt(feReq, np.array([0.5, 0.5]), "PK2Stress")
+if __name__ == "__main__":
+    dNonLin = nonlinElasticTest(easBool=False)
+    dEAS = nonlinElasticTest(easBool=True)
+
+    # Sanity check -> EAS shouldn't have the same result
+    assert not np.allclose(dNonLin, dEAS)
diff --git a/ikarus/solver/nonlinearsolver/newtonraphson.hh b/ikarus/solver/nonlinearsolver/newtonraphson.hh
index 9c390c427..c519f0d7c 100644
--- a/ikarus/solver/nonlinearsolver/newtonraphson.hh
+++ b/ikarus/solver/nonlinearsolver/newtonraphson.hh
@@ -184,17 +184,16 @@ public:
         linearSolver_.factorize(Ax);
         linearSolver_.solve(correction_, -rx);
         dNorm = correction_.norm();
-        updateFunction_(x, correction_);
       } else {
         correction_ = -linearSolver_(rx, Ax);
         dNorm       = norm(correction_);
-        updateFunction_(x, correction_);
       }
+      this->notify(CORRECTION_UPDATED, solverState);
+      updateFunction_(x, correction_);
       this->notify(CORRECTIONNORM_UPDATED, static_cast<double>(dNorm));
       this->notify(SOLUTION_CHANGED);
-      rx = residualFunction_(x);
-      Ax = jacobianFunction_(x);
-      this->notify(CORRECTION_UPDATED, solverState);
+      rx    = residualFunction_(x);
+      Ax    = jacobianFunction_(x);
       rNorm = norm(rx);
       this->notify(RESIDUALNORM_UPDATED, static_cast<double>(rNorm));
       this->notify(ITERATION_ENDED);
diff --git a/ikarus/solver/nonlinearsolver/solverinfos.hh b/ikarus/solver/nonlinearsolver/solverinfos.hh
index 8013d31b6..85e887907 100644
--- a/ikarus/solver/nonlinearsolver/solverinfos.hh
+++ b/ikarus/solver/nonlinearsolver/solverinfos.hh
@@ -3,7 +3,7 @@
 
 /**
  * \file solverinfos.hh
- * \brief Implementation of the Newton-Raphson method for solving nonlinear equations.
+ * \brief Implementation of the solver information returned by the nonlinear solvers.
  */
 
 #pragma once
diff --git a/ikarus/utils/concepts.hh b/ikarus/utils/concepts.hh
index 74af0578e..960a590d6 100644
--- a/ikarus/utils/concepts.hh
+++ b/ikarus/utils/concepts.hh
@@ -639,6 +639,28 @@ namespace Concepts {
     { s.stepSize } -> std::convertible_to<double>;
   };
 
+  namespace Formulations {
+    /**
+      \concept TotalLagrangian
+      \brief Concept to check if the underlying strain and stress tag correspond to a total Lagrangian formulation.
+      \tparam T1 The underlying strain tag.
+      \tparam T2 The underlying stress tag.
+      */
+    template <StrainTags T1, StressTags T2>
+    concept TotalLagrangian = (T1 == StrainTags::linear and T2 == StressTags::linear) or
+                              (T1 == StrainTags::greenLagrangian and T2 == StressTags::PK2);
+
+    /**
+      \concept TwoPoint
+      \brief Concept to check if the underlying strain and stress tag correspond to a two point formulation.
+      \tparam T1 The underlying strain tag.
+      \tparam T2 The underlying stress tag.
+      */
+    template <StrainTags T1, StressTags T2>
+    concept TwoPoint =
+        (T1 == StrainTags::deformationGradient or T1 == StrainTags::displacementGradient) and T2 == StressTags::PK1;
+  } // namespace Formulations
+
 } // namespace Concepts
 
 namespace traits {
diff --git a/ikarus/utils/tensorutils.hh b/ikarus/utils/tensorutils.hh
index 55f6ff853..406192ead 100644
--- a/ikarus/utils/tensorutils.hh
+++ b/ikarus/utils/tensorutils.hh
@@ -166,6 +166,7 @@ auto symTwoSlots(const Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim,
 /**
  * \brief Converts 2D indices to Voigt notation index.
  *  \ingroup tensor
+ * \tparam dim dimension (either 2d or 3d), defaults to 3
  * \param i Row index.
  * \param j Column index.
  * \return Eigen::Index Voigt notation index.
@@ -176,17 +177,28 @@ auto symTwoSlots(const Eigen::TensorFixedSize<ScalarType, Eigen::Sizes<dim, dim,
  *
  * If the input indices are not within the valid range (0, 1, 2), an assertion failure is triggered.
  */
+template <int dim = 3>
+requires(dim == 2 or dim == 3)
 constexpr Eigen::Index toVoigt(Eigen::Index i, Eigen::Index j) noexcept {
-  if (i == j) // _00 -> 0, _11 -> 1,  _22 -> 2
-    return i;
-  if ((i == 1 and j == 2) or (i == 2 and j == 1)) // _12 and _21 --> 3
-    return 3;
-  if ((i == 0 and j == 2) or (i == 2 and j == 0)) // _02 and _20 --> 4
-    return 4;
-  if ((i == 0 and j == 1) or (i == 1 and j == 0)) // _01 and _10 --> 5
-    return 5;
-  assert(i < 3 and j < 3 && "For Voigt notation the indices need to be 0,1 or 2.");
-  __builtin_unreachable();
+  if constexpr (dim == 2) {
+    if (i == j) // _00 -> 0, _11 -> 1
+      return i;
+    if ((i == 0 and j == 1) or (i == 1 and j == 0)) // _01 and _10 --> 2
+      return 2;
+    assert(i < 2 and j < 2 && "For Voigt notation the indices need to be 0 or 1.");
+    __builtin_unreachable();
+  } else {
+    if (i == j) // _00 -> 0, _11 -> 1,  _22 -> 2
+      return i;
+    if ((i == 1 and j == 2) or (i == 2 and j == 1)) // _12 and _21 --> 3
+      return 3;
+    if ((i == 0 and j == 2) or (i == 2 and j == 0)) // _02 and _20 --> 4
+      return 4;
+    if ((i == 0 and j == 1) or (i == 1 and j == 0)) // _01 and _10 --> 5
+      return 5;
+    assert(i < 3 and j < 3 && "For Voigt notation the indices need to be 0,1 or 2.");
+    __builtin_unreachable();
+  }
 }
 
 /**
@@ -374,29 +386,26 @@ auto fromVoigt(const Eigen::Matrix<ScalarType, 6, 6>& CVoigt) {
  *
  * \tparam Geometry The geometry type.
  * \param geometry Reference to the geometry object.
+ * \param pos The position where the transformation matrix is to be evaluated.
  * \return The transformation matrix for 2D elements.
  */
-template <typename Geometry>
-Eigen::Matrix3d calcTransformationMatrix2D(const Geometry& geometry) {
-  const auto& referenceElement = Dune::ReferenceElements<double, 2>::general(geometry.type());
-  const auto quadPos0          = referenceElement.position(0, 0);
+template <typename GEO>
+requires(GEO::mydimension == 2)
+Eigen::Matrix3d transformationMatrix(const GEO& geometry, const Dune::FieldVector<double, 2>& pos) {
+  const auto jacobian = toEigen(geometry.jacobianTransposed(pos)).eval();
 
-  const auto jacobianinvT0 = toEigen(geometry.jacobianInverseTransposed(quadPos0));
-  const auto detJ0         = geometry.integrationElement(quadPos0);
+  const auto J11 = jacobian(0, 0);
+  const auto J12 = jacobian(0, 1);
+  const auto J21 = jacobian(1, 0);
+  const auto J22 = jacobian(1, 1);
 
-  auto jaco = (jacobianinvT0).inverse().eval();
-  auto J11  = jaco(0, 0);
-  auto J12  = jaco(0, 1);
-  auto J21  = jaco(1, 0);
-  auto J22  = jaco(1, 1);
-
-  Eigen::Matrix3d T0{
+  Eigen::Matrix3d T{
       {      J11 * J11,       J12 * J12,             J11 * J12},
       {      J21 * J21,       J22 * J22,             J21 * J22},
       {2.0 * J11 * J21, 2.0 * J12 * J22, J21 * J12 + J11 * J22}
   };
 
-  return T0.inverse() * detJ0;
+  return T;
 }
 
 /**
@@ -407,39 +416,36 @@ Eigen::Matrix3d calcTransformationMatrix2D(const Geometry& geometry) {
  *
  * \tparam Geometry The geometry type.
  * \param geometry Reference to the geometry object.
+ * \param pos The position where the transformation matrix is to be evaluated.
  * \return The transformation matrix for 3D elements.
  */
-template <typename Geometry>
-Eigen::Matrix<double, 6, 6> calcTransformationMatrix3D(const Geometry& geometry) {
-  const auto& referenceElement = Dune::ReferenceElements<double, 3>::general(geometry.type());
-  const auto quadPos0          = referenceElement.position(0, 0);
-
-  const auto jacobianinvT0 = toEigen(geometry.jacobianInverseTransposed(quadPos0));
-  const auto detJ0         = geometry.integrationElement(quadPos0);
-
-  auto jaco = (jacobianinvT0).inverse().eval();
-  auto J11  = jaco(0, 0);
-  auto J12  = jaco(0, 1);
-  auto J13  = jaco(0, 2);
-  auto J21  = jaco(1, 0);
-  auto J22  = jaco(1, 1);
-  auto J23  = jaco(1, 2);
-  auto J31  = jaco(2, 0);
-  auto J32  = jaco(2, 1);
-  auto J33  = jaco(2, 2);
+template <typename GEO>
+requires(GEO::mydimension == 3)
+Eigen::Matrix<double, 6, 6> transformationMatrix(const GEO& geometry, const Dune::FieldVector<double, 3>& pos) {
+  const auto jacobian = toEigen(geometry.jacobianTransposed(pos)).eval();
+
+  const auto J11 = jacobian(0, 0);
+  const auto J12 = jacobian(0, 1);
+  const auto J13 = jacobian(0, 2);
+  const auto J21 = jacobian(1, 0);
+  const auto J22 = jacobian(1, 1);
+  const auto J23 = jacobian(1, 2);
+  const auto J31 = jacobian(2, 0);
+  const auto J32 = jacobian(2, 1);
+  const auto J33 = jacobian(2, 2);
 
   // clang-format off
-  Eigen::Matrix<double, 6, 6> T0  {
-    {J11 * J11,       J12 * J12,       J13 * J13,             J11 * J12,             J11 * J13,             J12 * J13},
-    {J21 * J21,       J22 * J22,       J23 * J23,             J21 * J22,             J21 * J23,             J22 * J23},
-    {J31 * J31,       J32 * J32,       J33 * J33,             J31 * J32,             J31 * J33,             J32 * J33},
-    {2.0 * J11 * J21, 2.0 * J12 * J22, 2.0 * J13 * J23, J11 * J22 + J21 * J12, J11 * J23 + J21 * J13, J12 * J23 + J22 * J13},
-    {2.0 * J11 * J31, 2.0 * J12 * J32, 2.0 * J13 * J33, J11 * J32 + J31 * J12, J11 * J33 + J31 * J13, J12 * J33 + J32 * J13},
-    {2.0 * J31 * J21, 2.0 * J32 * J22, 2.0 * J33 * J23, J31 * J22 + J21 * J32, J31 * J23 + J21 * J33, J32 * J23 + J22 * J33}
-  };
+   Eigen::Matrix<double, 6, 6> T  {
+     {J11 * J11,       J12 * J12,       J13 * J13,             J11 * J12,             J11 * J13,             J12 * J13},
+     {J21 * J21,       J22 * J22,       J23 * J23,             J21 * J22,             J21 * J23,             J22 * J23},
+     {J31 * J31,       J32 * J32,       J33 * J33,             J31 * J32,             J31 * J33,             J32 * J33},
+     {2.0 * J11 * J21, 2.0 * J12 * J22, 2.0 * J13 * J23, J11 * J22 + J21 * J12, J11 * J23 + J21 * J13, J12 * J23 + J22 * J13},
+     {2.0 * J11 * J31, 2.0 * J12 * J32, 2.0 * J13 * J33, J11 * J32 + J31 * J12, J11 * J33 + J31 * J13, J12 * J33 + J32 * J13},
+     {2.0 * J31 * J21, 2.0 * J32 * J22, 2.0 * J33 * J23, J31 * J22 + J21 * J32, J31 * J23 + J21 * J33, J32 * J23 + J22 * J33}
+   };
   // clang-format on
 
-  return T0.inverse() * detJ0;
+  return T;
 }
 
 } // namespace Ikarus
diff --git a/python/ikarus/finite_elements/__init__.py b/python/ikarus/finite_elements/__init__.py
index d1318e56c..e603ac7b0 100644
--- a/python/ikarus/finite_elements/__init__.py
+++ b/python/ikarus/finite_elements/__init__.py
@@ -3,6 +3,7 @@
 
 from dune.common.hashit import hashIt
 from ikarus.generator import MySimpleGenerator
+from ikarus import materials
 
 
 def registerPreElement(name, includes, element_type, *args):
@@ -146,7 +147,7 @@ def truss(youngs_modulus, cross_section):
     )
 
 
-def eas(numberofparameters):
+def eas(numberofparameters, strainTag="LinearStrain"):
     """
     @brief Creates an enhanced assumed strains pre-element.
 
@@ -155,7 +156,8 @@ def eas(numberofparameters):
     @return: The registered enhanced assumed strains pre-element function.
     """
     includes = ["ikarus/finiteelements/mechanics/enhancedassumedstrains.hh"]
-    element_type = "Ikarus::EnhancedAssumedStrainsPre"
+    formatted_strainTag = "Ikarus::EAS::" + strainTag
+    element_type = f"Ikarus::EnhancedAssumedStrainsPre<{formatted_strainTag}>"
     return registerPreElement(
         "EnhancedAssumedStrainsPre", includes, element_type, numberofparameters
     )
diff --git a/tests/src/CMakeLists.txt b/tests/src/CMakeLists.txt
index 2f8650757..34f2dfc7c 100644
--- a/tests/src/CMakeLists.txt
+++ b/tests/src/CMakeLists.txt
@@ -27,6 +27,7 @@ set(programSourceFiles
     testlinearelasticity.cpp
     testmanifolds.cpp
     testmaterial.cpp
+    testnonlineareas.cpp
     testnonlinearelasticityneohooke.cpp
     testnonlinearelasticitysvk.cpp
     testdifffunction.cpp
@@ -42,6 +43,7 @@ set(programSourceFiles
     testhyperelasticity.cpp
     testmaterialswithad.cpp
     testmueslimaterials.cpp
+    testcantileverbeam.cpp
 )
 
 set(TEST_DEPENDING_ON_LOCALFEFUNCTIONS
diff --git a/tests/src/checkfebyautodiff.hh b/tests/src/checkfebyautodiff.hh
index 006430ea3..0a5bbfd1a 100644
--- a/tests/src/checkfebyautodiff.hh
+++ b/tests/src/checkfebyautodiff.hh
@@ -16,16 +16,18 @@
 
 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& testName = "",
                             double tol = 1e-10) {
   double lambda = 7.3;
   auto fe       = Ikarus::makeFE(basis, std::forward<Skills>(skills));
   using FE      = decltype(fe);
 
   auto req = typename FE::Requirement();
-  req.insertGlobalSolution(d).insertParameter(lambda);
-  Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation" +
-                    messageIfFailed);
+  VectorType dZero;
+  dZero.setZero(d.size());
+
+  req.insertGlobalSolution(dZero).insertParameter(lambda);
+  Dune::TestSuite t("Check calculateMatrixImpl() and calculateVectorImpl() by Automatic Differentiation" + testName);
   for (auto element : elements(gridView)) {
     auto localView = basis.flat().localView();
     localView.bind(element);
@@ -33,6 +35,9 @@ auto checkFESByAutoDiffImpl(const GridView& gridView, const BasisHandler& basis,
 
     fe.bind(element);
 
+    fe.updateState(req, d); // here d = correction vector (DeltaD)
+    req.insertGlobalSolution(d).insertParameter(lambda);
+
     const std::string feClassName = Dune::className(fe);
 
     using AutoDiffBasedFE = Ikarus::AutoDiffFE<decltype(fe), true>;
@@ -49,14 +54,35 @@ auto checkFESByAutoDiffImpl(const GridView& gridView, const BasisHandler& basis,
     calculateMatrix(fe, req, affordance.matrixAffordance(), K);
     calculateMatrix(feAutoDiff, req, affordance.matrixAffordance(), KAutoDiff);
 
-    calculateVector(fe, req, affordance.vectorAffordance(), R);
-    calculateVector(feAutoDiff, req, affordance.vectorAffordance(), RAutoDiff);
-
-    checkScalars(t, calculateScalar(fe, req, affordance.scalarAffordance()),
-                 calculateScalar(feAutoDiff, req, affordance.scalarAffordance()),
-                 testLocation() + "Incorrect Energies." + feClassName, tol);
-    checkApproxVectors(t, R, RAutoDiff, testLocation() + "Incorrect residual vectors." + feClassName, tol);
-    checkApproxMatrices(t, K, KAutoDiff, testLocation() + "Incorrect stiffness matrices." + feClassName, tol);
+    checkApproxMatrices(t, K, KAutoDiff, testLocation() + "\nIncorrect stiffness matrices." + feClassName, tol);
+    checkSymmetricMatrix(t, K, tol, "K");
+    checkSymmetricMatrix(t, KAutoDiff, tol, "KAutoDiff");
+
+    if constexpr (requires { fe.numberOfEASParameters(); }) {
+      checkScalars(t, fe.numberOfEASParameters(), feAutoDiff.realFE().numberOfEASParameters(),
+                   "Incorrect number of EAS parameters");
+      checkApproxVectors(t, fe.alpha(), feAutoDiff.alpha(),
+                         testLocation() + "\nIncorrect internal variable alpha." + feClassName, tol);
+
+      if (fe.numberOfEASParameters() != 0) {
+        t.checkThrow<Dune::NotImplemented>(
+            [&]() { calculateVector(feAutoDiff, req, affordance.vectorAffordance(), RAutoDiff); },
+            "calculateVector with AutoDiff should throw a Dune::NotImplemented");
+
+        t.checkThrow<Dune::NotImplemented>(
+            [&]() { auto energyAutoDiff = calculateScalar(feAutoDiff, req, affordance.scalarAffordance()); },
+            "calculateScalar with AutoDiff should throw a Dune::NotImplemented");
+      }
+
+    } else {
+      calculateVector(fe, req, affordance.vectorAffordance(), R);
+      calculateVector(feAutoDiff, req, affordance.vectorAffordance(), RAutoDiff);
+
+      checkApproxVectors(t, R, RAutoDiff, testLocation() + "\nIncorrect residual vectors." + feClassName, tol);
+      checkScalars(t, calculateScalar(fe, req, affordance.scalarAffordance()),
+                   calculateScalar(feAutoDiff, req, affordance.scalarAffordance()),
+                   testLocation() + "\nIncorrect Energies." + feClassName, tol);
+    }
   }
 
   return t;
@@ -71,7 +97,7 @@ auto checkFESByAutoDiff(const GridView& gridView, const PreBasis& pb, Skills&& s
   d.setZero(basis.flat().dimension());
   t.subTest(checkFESByAutoDiffImpl(gridView, basis, skills, affordance, d, " Zero Displacements", tol));
   d.setRandom(basis.flat().dimension());
-  d *= 0.2; // to avoid a negative determinant of deformation gradient
+  d *= 0.01; // 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/testcantileverbeam.cpp b/tests/src/testcantileverbeam.cpp
new file mode 100644
index 000000000..dff4f919b
--- /dev/null
+++ b/tests/src/testcantileverbeam.cpp
@@ -0,0 +1,185 @@
+// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+#include <config.h>
+
+#include "testcommon.hh"
+#include "testhelpers.hh"
+
+#include <dune/common/test/testsuite.hh>
+#include <dune/functions/functionspacebases/basistags.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <Eigen/Core>
+
+#include <ikarus/assembler/assemblermanipulatorfuser.hh>
+#include <ikarus/assembler/simpleassemblers.hh>
+#include <ikarus/controlroutines/loadcontrol.hh>
+#include <ikarus/finiteelements/fefactory.hh>
+#include <ikarus/finiteelements/mechanics/enhancedassumedstrains.hh>
+#include <ikarus/finiteelements/mechanics/materials.hh>
+#include <ikarus/finiteelements/mechanics/nonlinearelastic.hh>
+#include <ikarus/solver/linearsolver/linearsolver.hh>
+#include <ikarus/solver/nonlinearsolver/newtonraphson.hh>
+#include <ikarus/solver/nonlinearsolver/nonlinearsolverfactory.hh>
+#include <ikarus/utils/basis.hh>
+#include <ikarus/utils/dirichletvalues.hh>
+#include <ikarus/utils/functionhelper.hh>
+#include <ikarus/utils/init.hh>
+#include <ikarus/utils/listener/controllogger.hh>
+#include <ikarus/utils/listener/controlvtkwriter.hh>
+#include <ikarus/utils/listener/nonlinearsolverlogger.hh>
+
+using namespace Ikarus;
+using Dune::TestSuite;
+
+template <typename ES, typename MAT>
+auto cantileverBeamResults() {
+  static_assert(MAT::isReduced,
+                "cantileverBeamResults are available only for a reduced material (planeStress or planeStrain).");
+  using namespace Materials;
+  using MATU = typename MAT::Underlying;
+  if constexpr (std::same_as<ES, EAS::GreenLagrangeStrain>) {
+    if constexpr (std::is_same_v<MATU, StVenantKirchhoff>) {
+      constexpr int expectedIterations = 80;
+      constexpr double expectedMaxDisp = 4.459851990257645;
+      return std::make_pair(expectedIterations, expectedMaxDisp);
+    } else if constexpr (std::is_same_v<MATU, NeoHooke>) {
+      constexpr int expectedIterations = 80;
+      constexpr double expectedMaxDisp = 4.479930218997457;
+      return std::make_pair(expectedIterations, expectedMaxDisp);
+    } else if constexpr (std::is_same_v<MATU, Hyperelastic<Deviatoric<BlatzKo>, Volumetric<VF0>>>) {
+      constexpr int expectedIterations = 80;
+      constexpr double expectedMaxDisp = 4.6877262928164365;
+      return std::make_pair(expectedIterations, expectedMaxDisp);
+    } else
+      static_assert(Dune::AlwaysFalse<MATU>::value,
+                    "Expected results are not available for the given material and enhancedStrain type.");
+  } else
+    static_assert(Dune::AlwaysFalse<MATU>::value,
+                  "Expected results are not available for the provided enhanced strain type.");
+}
+
+template <typename ES, typename MAT>
+auto cantileverBeamTest(const MAT& reducedMat) {
+  static_assert(MAT::isReduced,
+                "cantileverBeamTest is only valid for a reduced material (planeStress or planeStrain).");
+  TestSuite t("Cantilever Beam for Nonlinear EAS element with enhanced strain type as " + ES::name() +
+              " and material type as " + Dune::className<MAT>());
+  constexpr int gridDim = 2;
+  using Grid            = Dune::YaspGrid<gridDim>;
+  const double L        = 10;
+  const double h        = 2;
+
+  Dune::FieldVector<double, gridDim> bbox       = {L, h};
+  std::array<int, gridDim> elementsPerDirection = {10, 1};
+  auto grid                                     = std::make_shared<Grid>(bbox, elementsPerDirection);
+  auto gridView                                 = grid->leafGridView();
+  using namespace Dune::Functions::BasisFactory;
+  auto basis = Ikarus::makeBasis(gridView, power<2>(lagrange<1>(), FlatInterleaved()));
+
+  auto sk      = skills(nonLinearElastic(reducedMat), eas<ES>(4));
+  using FEType = decltype(makeFE(basis, sk));
+  std::vector<FEType> fes;
+
+  for (auto&& ge : elements(gridView)) {
+    fes.emplace_back(makeFE(basis, sk));
+    fes.back().bind(ge);
+  }
+
+  Ikarus::DirichletValues dirichletValues(basis.flat());
+
+  // fix left edge (x=0)
+  dirichletValues.fixBoundaryDOFs([&](auto& dirichletFlags, auto&& localIndex, auto&& localView, auto&& intersection) {
+    if (std::abs(intersection.geometry().center()[0]) < 1e-8)
+      dirichletFlags[localView.index(localIndex)] = true;
+  });
+
+  auto sparseFlatAssembler = makeSparseFlatAssembler(fes, dirichletValues);
+  auto sparseAssemblerAM   = makeAssemblerManipulator(*sparseFlatAssembler);
+
+  Eigen::VectorXd d;
+  d.setZero(basis.flat().size());
+  double lambda = 0.0;
+
+  auto req = typename FEType::Requirement(d, lambda);
+
+  sparseAssemblerAM->bind(req, Ikarus::AffordanceCollections::elastoStatics, DBCOption::Full);
+
+  // Apply constant point load at the top right corner
+  Dune::FieldVector<double, gridDim> topRightPos{L, h};
+  const auto globalIndices = utils::globalIndexFromGlobalPosition(basis.flat(), topRightPos);
+  auto pointLoad           = [&](const auto&, const auto& par, auto, auto, Eigen::VectorXd& vec) -> void {
+    auto loadFactor = par.parameter();
+    vec[globalIndices[1]] -= -loadFactor * 1.0;
+  };
+  sparseAssemblerAM->bind(pointLoad);
+
+  auto linSolver = LinearSolver(SolverTypeTag::sd_UmfPackLU);
+
+  AffordanceCollection elastoStaticsNoScalar(VectorAffordance::forces, MatrixAffordance::stiffness);
+
+  auto nonOp =
+      DifferentiableFunctionFactory::op(sparseAssemblerAM, elastoStaticsNoScalar, sparseAssemblerAM->dBCOption());
+
+  constexpr double tol = 1e-10;
+
+  auto nrConfig =
+      Ikarus::NewtonRaphsonConfig<decltype(linSolver)>{.parameters = {.tol = tol}, .linearSolver = linSolver};
+  NonlinearSolverFactory nrFactory(nrConfig);
+  auto nr = nrFactory.create(sparseAssemblerAM);
+
+  // Only when creating the control routine via the Factory, the elements get registered for correction update
+  // automatically.
+  auto lc = ControlRoutineFactory::create(LoadControlConfig{20, 0.0, 1.0}, nr, sparseFlatAssembler);
+
+  auto nonLinearSolverObserver = NonLinearSolverLogger();
+  auto controlLogger           = ControlLogger();
+  auto vtkWriter = ControlSubsamplingVertexVTKWriter<std::remove_cvref_t<decltype(basis.flat())>>(basis.flat(), d);
+  vtkWriter.setFileNamePrefix("CantileverNonlinearEAS");
+  vtkWriter.setFieldInfo("Displacement", Dune::VTK::FieldInfo::Type::vector, 2);
+
+  nonLinearSolverObserver.subscribeTo(lc.nonLinearSolver());
+  vtkWriter.subscribeTo(lc);
+  controlLogger.subscribeTo(lc);
+
+  const auto controlInfo = lc.run(req);
+  d                      = req.globalSolution();
+  lambda                 = req.parameter();
+
+  double expectedLambda                            = 1.0;
+  const auto [expectedIterations, expectedMaxDisp] = cantileverBeamResults<ES, MAT>();
+
+  t.check(controlInfo.success);
+  checkScalars(t, controlInfo.totalIterations, expectedIterations, " Number of iterations");
+  const auto maxDisp = std::ranges::max(d.cwiseAbs());
+
+  checkScalars(t, maxDisp, expectedMaxDisp, " Max. displacement", tol);
+  checkScalars(t, lambda, expectedLambda, " Load factor", tol);
+
+  return t;
+}
+
+int main(int argc, char** argv) {
+  Ikarus::init(argc, argv);
+  Dune::TestSuite t;
+  auto matParameter = toLamesFirstParameterAndShearModulus({.emodul = 100.0, .nu = 0.3});
+  Materials::StVenantKirchhoff matSVK(matParameter);
+  Materials::NeoHooke matNH(matParameter);
+  auto matBK         = Materials::makeBlatzKo(40.0);
+  auto reducedMatSVK = planeStrain(matSVK);
+  auto reducedMatNH  = planeStrain(matNH);
+  auto reducedMatBK  = planeStrain(matBK);
+
+  const auto testNonlinearEASFunc = [&]<typename ES>() {
+    t.subTest(cantileverBeamTest<ES>(reducedMatSVK));
+    t.subTest(cantileverBeamTest<ES>(reducedMatNH));
+    t.subTest(cantileverBeamTest<ES>(reducedMatBK));
+  };
+
+  testNonlinearEASFunc.operator()<EAS::GreenLagrangeStrain>();
+
+  return t.exit();
+}
diff --git a/tests/src/testeas.hh b/tests/src/testeas.hh
index d0067fece..a2b30f529 100644
--- a/tests/src/testeas.hh
+++ b/tests/src/testeas.hh
@@ -6,9 +6,10 @@
 #include "testcommon.hh"
 
 #include <ikarus/finiteelements/mechanics/enhancedassumedstrains.hh>
+#include <ikarus/finiteelements/mechanics/materials.hh>
 
 template <typename FE>
-requires(FE::template hasSkill<Ikarus::EnhancedAssumedStrainsPre::Skill>())
+requires(FE::template hasEAS<Ikarus::EAS::LinearStrain>)
 struct ElementTest<FE>
 {
   [[nodiscard]] static auto test() {
@@ -30,7 +31,7 @@ struct ElementTest<FE>
 
       for (auto& numberOfEASParameter : easParameters) {
         fe.setEASType(numberOfEASParameter);
-        auto messageIfFailed = "The numbers of EAS parameters are " + std::to_string(numberOfEASParameter) + ".";
+        auto messageIfFailed = "The number of EAS parameters are " + std::to_string(numberOfEASParameter) + ".";
         if (numberOfEASParameter == 0) {
           t.subTest(checkGradientOfElement(f, req, messageIfFailed));
           t.subTest(checkHessianOfElement(f, req, messageIfFailed));
@@ -68,10 +69,10 @@ struct ElementTest<FE>
           auto easVariantCopy    = fe.easVariant(); // This only test if the variant has a copy assignment operator
           const auto& easVariant = fe.easVariant();
           auto testM             = [&]<typename EAS>(const EAS& easFunction) {
-            typename EAS::MType MIntegrated;
+            typename EAS::AnsatzType MIntegrated;
             MIntegrated.setZero();
             for (const auto& gp : rule) {
-              const auto M = easFunction.calcM(gp.position());
+              const auto M = easFunction(gp.position());
 
               const double detJ = element.geometry().integrationElement(gp.position());
               MIntegrated += M * detJ * gp.weight();
diff --git a/tests/src/testhelpers.hh b/tests/src/testhelpers.hh
index c219ca2a0..a7c79b61e 100644
--- a/tests/src/testhelpers.hh
+++ b/tests/src/testhelpers.hh
@@ -85,6 +85,15 @@ void checkScalars(TestSuiteType& t, const ScalarType val, const ScalarType expec
                                 << " Actual:\t" << val << messageIfFailed;
 }
 
+template <typename TestSuiteType, int rows, int cols>
+requires(rows == cols)
+void checkSymmetricMatrix(TestSuiteType& t, const Eigen::Matrix<double, rows, cols>& mat, double tol,
+                          const std::string& matrixName = "") {
+  t.check(isApproxSame(mat, mat.transpose(), tol), "Matrix " + matrixName + " is not symmetric.\n")
+      << "Matrix " + matrixName + " is \n"
+      << mat;
+}
+
 template <typename TestSuiteType, typename ScalarType>
 requires(not std::is_integral_v<ScalarType>)
 void checkScalars(TestSuiteType& t, const ScalarType val, const ScalarType expectedVal,
diff --git a/tests/src/testnonlineareas.cpp b/tests/src/testnonlineareas.cpp
new file mode 100644
index 000000000..fbbcbb5f8
--- /dev/null
+++ b/tests/src/testnonlineareas.cpp
@@ -0,0 +1,219 @@
+// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
+// SPDX-License-Identifier: LGPL-3.0-or-later
+
+#include <config.h>
+
+#include "checkfebyautodiff.hh"
+#include "testcommon.hh"
+#include "testhelpers.hh"
+
+#include <dune/common/test/testsuite.hh>
+#include <dune/functions/functionspacebases/basistags.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <Eigen/Core>
+
+#include <ikarus/assembler/simpleassemblers.hh>
+#include <ikarus/finiteelements/autodifffe.hh>
+#include <ikarus/finiteelements/fefactory.hh>
+#include <ikarus/finiteelements/mechanics/enhancedassumedstrains.hh>
+#include <ikarus/finiteelements/mechanics/materials.hh>
+#include <ikarus/finiteelements/mechanics/nonlinearelastic.hh>
+#include <ikarus/solver/linearsolver/linearsolver.hh>
+#include <ikarus/utils/basis.hh>
+#include <ikarus/utils/dirichletvalues.hh>
+#include <ikarus/utils/functionhelper.hh>
+#include <ikarus/utils/init.hh>
+
+using namespace Ikarus;
+using Dune::TestSuite;
+
+template <int order, int gridDim, typename ES>
+auto easParametersFunc() {
+  const std::string& errorMessage =
+      "easParametersFunc is not implemented for the provided EnhancedStrainTag: " + ES::name();
+  if constexpr (order == 1) {
+    if constexpr (std::same_as<ES, EAS::GreenLagrangeStrain>) {
+      std::array<int, (gridDim == 2) ? 4 : 3> easParameters;
+      if constexpr (gridDim == 2)
+        easParameters = {0, 4, 5, 7};
+      else
+        easParameters = {0, 9, 21};
+      return easParameters;
+    } else
+      DUNE_THROW(Dune::NotImplemented, errorMessage);
+  } else {
+    if constexpr (std::same_as<ES, EAS::GreenLagrangeStrain>) {
+      static_assert(gridDim == 2);
+      std::array<int, 2> easParameters = {0, 11};
+      return easParameters;
+    } else
+      DUNE_THROW(Dune::NotImplemented, errorMessage);
+  }
+}
+
+template <int order, int gridDim, typename ES, typename TestSuitType, typename MAT>
+void easAutoDiffTest(TestSuitType& t, const MAT& mat) {
+  if constexpr (order == 2 and not std::same_as<ES, EAS::GreenLagrangeStrain>)
+    return;
+  else {
+    using namespace Dune::Functions::BasisFactory;
+    const auto easParameters = easParametersFunc<order, gridDim, ES>();
+    auto grid                = createUGGridFromCorners<gridDim>(CornerDistortionFlag::randomlyDistorted);
+    auto gridView            = grid->leafGridView();
+
+    for (const int numberOfEASParameters : easParameters) {
+      auto sk = skills(nonLinearElastic(mat), eas<ES>(numberOfEASParameters));
+      t.subTest(checkFESByAutoDiff(gridView, power<gridDim>(lagrange<order>()), sk,
+                                   AffordanceCollections::elastoStatics,
+                                   " (numberOfEASParameters = " + std::to_string(numberOfEASParameters) +
+                                       ") and EnhancedStrainTag = " + ES::name()));
+      if (numberOfEASParameters != 0) {
+        t.template checkThrow<Dune::NotImplemented>(
+            [&]() {
+              checkFESByAutoDiff(gridView, power<gridDim>(lagrange<3>()), sk, AffordanceCollections::elastoStatics,
+                                 " (numberOfEASParameters = " + std::to_string(numberOfEASParameters) +
+                                     ") and EnhancedStrainTag = " + ES::name());
+            },
+            "EAS method should have failed for a third-order basis.");
+      }
+    }
+  }
+}
+
+template <int order, int gridDim, typename ES, typename MAT>
+auto checkOrthogonalityCondition(const MAT& mat) {
+  TestSuite t("Orthogonality Condition Test for EAS Method");
+  if constexpr (order == 2 and not std::same_as<ES, EAS::GreenLagrangeStrain>)
+    return t;
+  else {
+    using namespace Dune::Functions::BasisFactory;
+    const auto easParameters  = easParametersFunc<order, gridDim, ES>();
+    auto grid                 = createUGGridFromCorners<gridDim>(CornerDistortionFlag::randomlyDistorted);
+    auto gridView             = grid->leafGridView();
+    auto basis                = Ikarus::makeBasis(gridView, power<gridDim>(lagrange<order>()));
+    constexpr auto affordance = AffordanceCollections::elastoStatics;
+
+    const auto orthogonalityFunc = [&]<typename FE>(const FE& fe) {
+      const auto& localView  = fe.localView();
+      const auto& element    = localView.element();
+      auto rule              = Dune::QuadratureRules<double, gridDim>::rule(localView.element().type(), 2 * order);
+      const auto& easVariant = fe.easVariant();
+      const int numberOfEASParameters = fe.numberOfEASParameters();
+      auto integrateFunctor           = [&]<typename EASF>(const EASF& easFunction) {
+        typename EASF::AnsatzType integratedStrain;
+        integratedStrain.setZero();
+        for (const auto& gp : rule) {
+          const auto M      = easFunction(gp.position());
+          const double detJ = element.geometry().integrationElement(gp.position());
+          if constexpr (std::same_as<ES, EAS::GreenLagrangeStrain>)
+            integratedStrain += M * detJ * gp.weight();
+          else
+            for (const auto& H : M)
+              integratedStrain += H * detJ * gp.weight();
+        }
+        t.check(integratedStrain.isZero()) << "Ansatz function for the enhanced strains in the EAS method should be "
+                                                        "zero when integrated over the domain. (numberOfEASParameters = " +
+                                                  std::to_string(numberOfEASParameters) +
+                                                  ") and EnhancedStrainTag = " + ES::name();
+      };
+      easVariant(integrateFunctor);
+    };
+
+    for (const int numberOfEASParameters : easParameters) {
+      auto sk      = skills(nonLinearElastic(mat), eas<ES>(numberOfEASParameters));
+      auto fe      = makeFE(basis, sk);
+      auto element = gridView.template begin<0>();
+      fe.bind(*element);
+      orthogonalityFunc(fe);
+    }
+    return t;
+  }
+}
+
+template <int order, int gridDim, typename ES, typename MAT>
+auto recoverNonlinearElastic(const MAT& mat) {
+  TestSuite t(
+      "Recover NonLinearElastic Test for EAS Method with numberOfEASParameters = 0 with enhanced strain type as " +
+      ES::name() + " and material type as " + Dune::className<MAT>());
+  using namespace Dune::Functions::BasisFactory;
+  auto grid                 = createUGGridFromCorners<gridDim>(CornerDistortionFlag::randomlyDistorted);
+  auto gridView             = grid->leafGridView();
+  auto basis                = Ikarus::makeBasis(gridView, power<gridDim>(lagrange<order>()));
+  constexpr auto affordance = AffordanceCollections::elastoStatics;
+  auto element              = gridView.template begin<0>();
+  auto nDOF                 = basis.flat().size();
+  const double tol          = 1e-10;
+
+  auto fe1 = makeFE(basis, skills(nonLinearElastic(mat)));
+  auto fe2 = makeFE(basis, skills(nonLinearElastic(mat), eas<ES>(0)));
+  fe1.bind(*element);
+  fe2.bind(*element); // updateState need not be called as numberOfEASParameters = 0
+
+  Eigen::VectorXd d;
+  d.setRandom(nDOF);
+  d *= 0.01; // to avoid a negative determinant of deformation gradient
+  double lambda = 0.0;
+  auto req      = typename decltype(fe1)::Requirement(d, lambda);
+
+  Eigen::MatrixXd K1, K2;
+  K1.setZero(nDOF, nDOF);
+  K2.setZero(nDOF, nDOF);
+
+  Eigen::VectorXd R1, R2;
+  R1.setZero(nDOF);
+  R2.setZero(nDOF);
+
+  calculateMatrix(fe1, req, affordance.matrixAffordance(), K1);
+  calculateMatrix(fe2, req, affordance.matrixAffordance(), K2);
+
+  calculateVector(fe1, req, affordance.vectorAffordance(), R1);
+  calculateVector(fe2, req, affordance.vectorAffordance(), R2);
+
+  checkApproxMatrices(t, K1, K2, testLocation() + "\nIncorrect stiffness matrices.", tol);
+  checkSymmetricMatrix(t, K1, tol, "K1");
+  checkSymmetricMatrix(t, K2, tol, "K2");
+  checkApproxVectors(t, R1, R2, testLocation() + "\nIncorrect residual vectors.", tol);
+  return t;
+}
+
+int main(int argc, char** argv) {
+  Ikarus::init(argc, argv);
+  Dune::TestSuite t("Nonlinear EAS Test");
+  auto matParameter = toLamesFirstParameterAndShearModulus({.emodul = 100.0, .nu = 0.3});
+  Materials::StVenantKirchhoff matSVK(matParameter);
+  Materials::NeoHooke matNH(matParameter);
+  auto matBK         = Materials::makeBlatzKo(40.0);
+  auto reducedMatSVK = planeStrain(matSVK);
+  auto reducedMatNH  = planeStrain(matNH);
+  auto reducedMatBK  = planeStrain(matBK);
+
+  const auto testNonlinearEASFunc = [&]<typename ES>() {
+    easAutoDiffTest<1, 2, ES>(t, reducedMatSVK);
+    easAutoDiffTest<1, 3, ES>(t, matSVK);
+    easAutoDiffTest<1, 2, ES>(t, reducedMatNH);
+    easAutoDiffTest<1, 3, ES>(t, matNH);
+    easAutoDiffTest<1, 2, ES>(t, reducedMatBK);
+    easAutoDiffTest<1, 3, ES>(t, matBK);
+    easAutoDiffTest<2, 2, ES>(t, reducedMatSVK);
+    easAutoDiffTest<2, 2, ES>(t, reducedMatNH);
+    easAutoDiffTest<2, 2, ES>(t, reducedMatBK);
+
+    t.subTest(recoverNonlinearElastic<1, 2, ES>(reducedMatSVK));
+    t.subTest(recoverNonlinearElastic<1, 3, ES>(matSVK));
+    t.subTest(recoverNonlinearElastic<1, 2, ES>(reducedMatNH));
+    t.subTest(recoverNonlinearElastic<1, 3, ES>(matNH));
+    t.subTest(recoverNonlinearElastic<2, 2, ES>(reducedMatSVK));
+    t.subTest(recoverNonlinearElastic<2, 2, ES>(reducedMatNH));
+
+    t.subTest(checkOrthogonalityCondition<1, 2, ES>(reducedMatSVK));
+    t.subTest(checkOrthogonalityCondition<1, 3, ES>(matSVK));
+    t.subTest(checkOrthogonalityCondition<2, 2, ES>(reducedMatSVK));
+  };
+
+  testNonlinearEASFunc.operator()<EAS::GreenLagrangeStrain>();
+
+  return t.exit();
+}
\ No newline at end of file
diff --git a/tests/src/testtruss.cpp b/tests/src/testtruss.cpp
index a218a670a..7b829bf60 100644
--- a/tests/src/testtruss.cpp
+++ b/tests/src/testtruss.cpp
@@ -120,7 +120,7 @@ static auto vonMisesTrussTest() {
   auto nr = nrFactory.create(denseFlatAssembler);
 
   /// Create Observer to write information of the non-linear solver
-  auto f = Ikarus::DifferentiableFunctionFactory::op(denseFlatAssembler);
+  auto f = DifferentiableFunctionFactory::op(denseFlatAssembler);
 
   t.check(utils::checkGradient(f, req, {.draw = false, .writeSlopeStatementIfFailed = true}))
       << "Check gradient failed";