diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index 1ecbe1db54c85a4282ed1a55a043115f09bc7525..486ccb2df2a920dd4488f47dfa39e46f198609bf 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -7,6 +7,7 @@ add_library(
   CoreModule.cpp
   DevUtilsModule.cpp
   LinearSolverModule.cpp
+  LocalFSIModule.cpp
   MathFunctionRegisterForVh.cpp
   MathModule.cpp
   MeshModule.cpp
diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..98d8d86cf98c26faf0975ea03c65ae3a6abe9fa5
--- /dev/null
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -0,0 +1,582 @@
+#include <language/modules/LocalFSIModule.hpp>
+
+#include <language/modules/SchemeModuleTypes.hpp>
+#include <language/modules/MeshModuleTypes.hpp>
+
+#include <language/modules/ModuleRepository.hpp>
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
+
+
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
+
+
+#include <scheme/LocalDtAcousticSolver.hpp>
+#include <scheme/LocalDtHyperelasticSolver.hpp>
+
+#include <scheme/Order2AcousticSolver.hpp>
+#include <scheme/Order2HyperelasticSolver.hpp>
+#include <scheme/Order2LocalDtHyperelasticSolver.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+
+#include <memory>
+
+LocalFSIModule::LocalFSIModule()
+{
+
+  this->_addBuiltinFunction("coupling",
+                            std::function(
+
+                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+                                 const size_t& label) -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                return std::make_shared<CouplingBoundaryConditionDescriptor>(boundary, label);
+                              }
+
+                              ));
+                              
+
+    this->_addBuiltinFunction("eucclhyd_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return Order2AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
+                                  .solver()
+                                  .apply(Order2AcousticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, c, p,
+                                         bc_descriptor_list);
+                              }
+
+                              ));
+
+    this->_addBuiltinFunction("glace_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return Order2AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
+                                  .solver()
+                                  .apply(Order2AcousticSolverHandler::SolverType::Glace, dt, rho, u, E, c, p,
+                                         bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_glace_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const size_t& q)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
+                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
+                                  .solver()
+                                  .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1, u2,
+                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_glace_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
+                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
+                                  .solver()
+                                  .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1, rho1, rho2, u1, u2, E1,
+                                         E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_eucclhyd_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const size_t& q)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
+                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
+                                  .solver()
+                                  .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1, u2,
+                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_eucclhyd_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
+                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
+                                  .solver()
+                                  .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1, u2,
+                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
+                                 const double& lambda,
+                                 const double& mu,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
+                                         sigma, bc_descriptor_list, chi_solid, lambda, mu, 0, 0, 0, 0);
+                              }
+
+                              ));
+
+    this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
+                                 const double& lambda,
+                                 const double& mu,
+                                 const double& gamma,
+                                 const double& p_inf,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
+                                         sigma, bc_descriptor_list, chi_solid, lambda, mu, gamma, p_inf, 0, 0);
+                              }
+
+                              ));
+
+        this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook2_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
+                                 const double& lambda,
+                                 const double& mu,
+                                 const double& gamma_f,
+                                 const double& p_inf_f,
+                                 const double& gamma_s,
+                                 const double& p_inf_s,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
+                                  .solver()
+                                  .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 2, dt, rho, u, E, CG, aL, aT,
+                                         sigma, bc_descriptor_list, chi_solid, lambda, mu, gamma_f, p_inf_f, gamma_s, p_inf_s);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& mu, const double& lambda, const double& dt1, const size_t& q)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         bc_descriptor_list2, mu, lambda);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& mu, const double& lambda, const double& dt1)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         bc_descriptor_list2, mu, lambda);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& lambda1, const double& mu1,
+                                 const double& lambda2, const double& mu2, 
+                                 const double& gamma1,  const double& gamma2 , const double& dt1)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, 0, gamma2, 0);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& lambda1, const double& mu1,
+                                 const double& lambda2, const double& mu2, 
+                                 const double& gamma1,  const double& p_inf1,
+                                 const double& gamma2,  const double& p_inf2, const double& dt1)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, p_inf1, gamma2, p_inf2);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& mu, const double& lambda, const double& dt1, const size_t& q)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         bc_descriptor_list2, mu, lambda);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& mu, const double& lambda, const double& dt1)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         bc_descriptor_list2, mu, lambda);
+                              }
+
+                              ));
+}
+
+void
+LocalFSIModule::registerOperators() const
+{
+
+}
+
+void
+LocalFSIModule::registerCheckpointResume() const
+{
+
+}
+
+ModuleRepository::Subscribe<LocalFSIModule> localFSI_module;
diff --git a/src/language/modules/LocalFSIModule.hpp b/src/language/modules/LocalFSIModule.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a3576c454618de1e391d67e5b7b702972cd7d047
--- /dev/null
+++ b/src/language/modules/LocalFSIModule.hpp
@@ -0,0 +1,25 @@
+#ifndef LOCALFSI_MODULE_HPP
+#define LOCALFSI_MODULE_HPP
+
+#include <language/modules/BuiltinModule.hpp>
+
+class LocalFSIModule : public BuiltinModule
+{
+  friend class MathFunctionRegisterForVh;
+
+ public:
+  std::string_view
+  name() const final
+  {
+    return "localFSI";
+  }
+
+  void registerOperators() const final;
+  void registerCheckpointResume() const final;
+
+  LocalFSIModule();
+
+  ~LocalFSIModule() = default;
+};
+
+#endif   // LOCALFSI_MODULE_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b9edcc21841ecf39cb3736bfea972facbd153d43..40cc703cc6dcc0029da81366fb95134fc1329599 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -26,7 +26,6 @@
 #include <mesh/MeshTraits.hpp>
 #include <scheme/AcousticSolver.hpp>
 #include <scheme/AxisBoundaryConditionDescriptor.hpp>
-#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
@@ -46,13 +45,8 @@
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/InflowBoundaryConditionDescriptor.hpp>
-#include <scheme/LocalDtAcousticSolver.hpp>
-#include <scheme/LocalDtHyperelasticSolver.hpp>
 #include <scheme/InflowListBoundaryConditionDescriptor.hpp>
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
-#include <scheme/Order2AcousticSolver.hpp>
-#include <scheme/Order2HyperelasticSolver.hpp>
-#include <scheme/Order2LocalDtHyperelasticSolver.hpp>
 #include <scheme/OutflowBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VariableBCDescriptor.hpp>
@@ -359,16 +353,6 @@ SchemeModule::SchemeModule()
 
                                           ));
 
-  this->_addBuiltinFunction("coupling",
-                            std::function(
-
-                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
-                                 const size_t& label) -> std::shared_ptr<const IBoundaryConditionDescriptor> {
-                                return std::make_shared<CouplingBoundaryConditionDescriptor>(boundary, label);
-                              }
-
-                              ));
-                              
   this->_addBuiltinFunction("wall", std::function(
 
                                       [](std::shared_ptr<const IBoundaryDescriptor> boundary)
@@ -539,184 +523,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-    this->_addBuiltinFunction("eucclhyd_solver_order2",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return Order2AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
-                                  .solver()
-                                  .apply(Order2AcousticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, c, p,
-                                         bc_descriptor_list);
-                              }
-
-                              ));
-
-    this->_addBuiltinFunction("glace_solver_order2",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return Order2AcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
-                                  .solver()
-                                  .apply(Order2AcousticSolverHandler::SolverType::Glace, dt, rho, u, E, c, p,
-                                         bc_descriptor_list);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_glace_solver",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2,
-                                 const size_t& q)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
-                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
-                                  .solver()
-                                  .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1, u2,
-                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_glace_solver",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
-                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
-                                  .solver()
-                                  .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1, rho1, rho2, u1, u2, E1,
-                                         E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_eucclhyd_solver",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2,
-                                 const size_t& q)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
-                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
-                                  .solver()
-                                  .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1, u2,
-                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_eucclhyd_solver",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const double& dt1, const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& c2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),
-                                                                    getCommonMesh({rho2, u2, E2, c2, p2})}
-                                  .solver()
-                                  .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1, u2,
-                                         E1, E2, c1, c2, p1, p2, bc_descriptor_list1, bc_descriptor_list2);
-                              }
-
-                              ));
-
   this->_addBuiltinFunction("apply_acoustic_fluxes",
                             std::function(
 
@@ -826,359 +632,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
-                                 const double& lambda,
-                                 const double& mu,
-                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
-                                  .solver()
-                                  .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
-                                         sigma, bc_descriptor_list, chi_solid, lambda, mu, 0, 0, 0, 0);
-                              }
-
-                              ));
-
-    this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook_solver_order2",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
-                                 const double& lambda,
-                                 const double& mu,
-                                 const double& gamma,
-                                 const double& p_inf,
-                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
-                                  .solver()
-                                  .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt, rho, u, E, CG, aL, aT,
-                                         sigma, bc_descriptor_list, chi_solid, lambda, mu, gamma, p_inf, 0, 0);
-                              }
-
-                              ));
-
-        this->_addBuiltinFunction("hyperelastic_eucclhyd_neohook2_solver_order2",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& chi_solid,
-                                 const double& lambda,
-                                 const double& mu,
-                                 const double& gamma_f,
-                                 const double& p_inf_f,
-                                 const double& gamma_s,
-                                 const double& p_inf_s,
-                                 const double& dt) -> std::tuple<std::shared_ptr<const MeshVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return Order2HyperelasticSolverHandler{getCommonMesh({rho, u, E, CG, aL, aT, sigma})}
-                                  .solver()
-                                  .apply(Order2HyperelasticSolverHandler::SolverType::Eucclhyd, 2, dt, rho, u, E, CG, aL, aT,
-                                         sigma, bc_descriptor_list, chi_solid, lambda, mu, gamma_f, p_inf_f, gamma_s, p_inf_s);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2,
-                                 const double& mu, const double& lambda, const double& dt1, const size_t& q)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
-                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
-                                                                        getCommonMesh(
-                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
-                                  .solver()
-                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1,
-                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
-                                         bc_descriptor_list2, mu, lambda);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2,
-                                 const double& mu, const double& lambda, const double& dt1)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
-                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
-                                                                        getCommonMesh(
-                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
-                                  .solver()
-                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1,
-                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
-                                         bc_descriptor_list2, mu, lambda);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2,
-                                 const double& lambda1, const double& mu1,
-                                 const double& lambda2, const double& mu2, 
-                                 const double& gamma1,  const double& gamma2 , const double& dt1)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
-                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
-                                                                        getCommonMesh(
-                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
-                                  .solver()
-                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
-                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
-                                         bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, 0, gamma2, 0);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_neohook_solver_order2",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2,
-                                 const double& lambda1, const double& mu1,
-                                 const double& lambda2, const double& mu2, 
-                                 const double& gamma1,  const double& p_inf1,
-                                 const double& gamma2,  const double& p_inf2, const double& dt1)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return Order2LocalDtHyperelasticSolverHandler{getCommonMesh(
-                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
-                                                                        getCommonMesh(
-                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
-                                  .solver()
-                                  .apply(Order2LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, 1, dt1, rho1, rho2, u1,
-                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
-                                         bc_descriptor_list2, lambda1, mu1, lambda2, mu2, gamma1, p_inf1, gamma2, p_inf2);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2,
-                                 const double& mu, const double& lambda, const double& dt1, const size_t& q)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
-                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
-                                                                        getCommonMesh(
-                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
-                                  .solver()
-                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1,
-                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
-                                         bc_descriptor_list2, mu, lambda);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list1,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list2,
-                                 const double& mu, const double& lambda, const double& dt1)
-                                -> std::tuple<
-                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>,
-                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
-                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
-                                                                        getCommonMesh(
-                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
-                                  .solver()
-                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, rho1, rho2, u1,
-                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
-                                         bc_descriptor_list2, mu, lambda);
-                              }
-
-                              ));
-
   this->_addBuiltinFunction("apply_hyperelastic_fluxes",
                             std::function(
 
diff --git a/src/scheme/Order2HyperelasticSolver.cpp b/src/scheme/Order2HyperelasticSolver.cpp
index e771f81b8abbdf628324cea480c2783c40fd6463..ec3ec68c8914b6b359bf484a278d102e1ecf726a 100644
--- a/src/scheme/Order2HyperelasticSolver.cpp
+++ b/src/scheme/Order2HyperelasticSolver.cpp
@@ -347,19 +347,30 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
     const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
     const NodeValue<const Rd>& xr        = mesh.xr();
 
-    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
-    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+    //const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+    //const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
     NodeValuePerCell<Rd> F{mesh.connectivity()};
+    //parallel_for(
+    //  mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+    //    const auto& node_to_cell                   = node_to_cell_matrix[r];
+    //    const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+    //
+    //    for (size_t j = 0; j < node_to_cell.size(); ++j) {
+    //      const CellId J       = node_to_cell[j];
+    //      const unsigned int R = node_local_number_in_its_cells[j];
+    //      F(J, R)              = -Ajr(J, R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma[J](xr[r]) * Cjr(J, R);
+    //    }
+    //  });
+
     parallel_for(
-      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        const auto& node_to_cell                   = node_to_cell_matrix[r];
-        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
 
-        for (size_t j = 0; j < node_to_cell.size(); ++j) {
-          const CellId J       = node_to_cell[j];
-          const unsigned int R = node_local_number_in_its_cells[j];
-          F(J, R)              = -Ajr(J, R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_sigma[J](xr[r]) * Cjr(J, R);
+        for(size_t r = 0; r < cell_nodes.size(); ++r) {
+          F(j,r) = -Ajr(j,r) * (DPk_uh[j](xr[cell_nodes[r]]) - ur[cell_nodes[r]]) + DPk_sigma[j](xr[cell_nodes[r]])*Cjr(j,r);
         }
       });
 
@@ -791,32 +802,19 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
     //
     //     }
     //);
-
-    std::cout << "Reconstruction"
-              << "\n";
     auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v, sigma_v);
     auto DPk_uh         = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
     auto DPk_sigmah     = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const Rdxd>>();
-    std::cout << "Reconstruction Ok ! "
-              << "\n";
 
-    std::cout << "Limitation"
-              << "\n";
     DiscreteFunctionDPk<Dimension, Rd> u_lim       = copy(DPk_uh);
     DiscreteFunctionDPk<Dimension, Rdxd> sigma_lim = copy(DPk_sigmah);
     this->_vector_limiter(mesh, u, u_lim);
     this->_tensor_limiter(mesh, sigma, sigma_lim);
-    std::cout << "Limitation OK !"
-              << "\n";
 
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
 
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
-    std::cout << "br"
-              << "\n";
     NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, sigma_lim);
-    std::cout << "br Ok !"
-              << "\n";
 
     const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
     this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
@@ -825,11 +823,7 @@ class Order2HyperelasticSolverHandler::Order2HyperelasticSolver final
     synchronize(br);
 
     NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
-    std::cout << "Fjr"
-              << "\n";
     NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim);
-    std::cout << "Fjr Ok ! "
-              << "\n";
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));