diff --git a/cmake/FindParMETIS.cmake b/cmake/FindParMETIS.cmake
index bf6d06c5a0de6c3cf9993fc587a1ab3fd26dd28b..d9b91d33d92617f3282d4eac2184549b62f62418 100644
--- a/cmake/FindParMETIS.cmake
+++ b/cmake/FindParMETIS.cmake
@@ -6,11 +6,11 @@ find_path(PARMETIS_INCLUDE_DIR parmetis.h
 
 if(EXISTS "${PARMETIS_INCLUDE_DIR}/parmetis.h")
   message(STATUS "Found parmetis.h in ${PARMETIS_INCLUDE_DIR}")
-  find_library(LIB_PARMETIS parmetis ${PARMETIS_LIBDIR})
+  find_library(LIB_PARMETIS parmetis $ENV{PARMETIS_LIBDIR})
   if("${LIB_PARMETIS}" STREQUAL "LIB_PARMETIS-NOTFOUND")
     message(WARNING "** Could not find parmetis library.\n** Is PARMETIS_LIBDIR correctly set (Actual: \"$ENV{PARMETIS_LIBDIR}\")?")
   endif()
-  find_library(LIB_METIS metis ${METIS_LIBDIR})
+  find_library(LIB_METIS metis $ENV{METIS_LIBDIR})
   if("${LIB_PARMETIS}" STREQUAL "LIB_METIS-NOTFOUND")
     message(WARNING "** Could not find metis library.\n** Is METIS_LIBDIR correctly set (Actual: \"$ENV{METIS_LIBDIR}\")?")
   endif()
@@ -22,12 +22,6 @@ if(EXISTS "${PARMETIS_INCLUDE_DIR}/parmetis.h")
     else()
       message(WARNING "** Could not find metis.h.\n** Is METIS_INCDIR correctly set (Actual: \"$ENV{METIS_INCDIR}\")?")
   endif()
-  include_directories(SYSTEM ${PARMETIS_INCLUDE_DIR})
-  include_directories(SYSTEM ${METIS_INCDIR})
-  link_directories(${PARMETIS_LIBDIR})
-  link_directories(${METIS_LIBDIR})
-  set(PARMETIS_LIBRARIES ${LIB_PARMETIS} ${LIB_METIS})
-  message(STATUS "Found parmetis/metis libraries ${PARMETIS_LIBRARIES}")
 else()
   message(WARNING "** Could not find parmetis.h.\n** Is PARMETIS_INCDIR correctly set (Actual: \"$ENV{PARMETIS_INCDIR}\")?")
 endif()
diff --git a/src/algebra/CMakeLists.txt b/src/algebra/CMakeLists.txt
index 847a47dcfb8830904f2f2912acf46f470c6eb34b..74f40290a2f9250cc2296ec62358b44c366ffc85 100644
--- a/src/algebra/CMakeLists.txt
+++ b/src/algebra/CMakeLists.txt
@@ -2,7 +2,6 @@
 
 add_library(
   PugsAlgebra
-  LeastSquareSolver.cpp
   EigenvalueSolver.cpp
   LinearSolver.cpp
   LinearSolverOptions.cpp
diff --git a/src/algebra/LeastSquareSolver.cpp b/src/algebra/LeastSquareSolver.cpp
deleted file mode 100644
index d376881c629a66d6c2864d882956ee83cd4fc5de..0000000000000000000000000000000000000000
--- a/src/algebra/LeastSquareSolver.cpp
+++ /dev/null
@@ -1,27 +0,0 @@
-#include <algebra/LeastSquareSolver.hpp>
-
-#include <algebra/CG.hpp>
-
-template <typename T>
-void
-LeastSquareSolver::solveLocalSystem(const SmallMatrix<T>& A, SmallVector<T>& x, const SmallVector<T>& b)
-{
-  if (A.numberOfRows() >= A.numberOfColumns()) {
-    const SmallMatrix tA = transpose(A);
-    const SmallMatrix B  = tA * A;
-    const SmallVector y  = tA * b;
-    CG{B, x, y, 1e-12, 10, false};
-  } else {
-    const SmallMatrix tA = transpose(A);
-    const SmallMatrix B  = A * tA;
-    SmallVector<double> y{A.numberOfRows()};
-    y = zero;
-    CG{B, y, b, 1e-12, 10, false};
-
-    x = tA * y;
-  }
-}
-
-template void LeastSquareSolver::solveLocalSystem(const SmallMatrix<double>&,
-                                                  SmallVector<double>&,
-                                                  const SmallVector<double>&);
diff --git a/src/algebra/LeastSquareSolver.hpp b/src/algebra/LeastSquareSolver.hpp
deleted file mode 100644
index b23ec9b7d160c43ee7a4f9e1b38251c8489a9f20..0000000000000000000000000000000000000000
--- a/src/algebra/LeastSquareSolver.hpp
+++ /dev/null
@@ -1,17 +0,0 @@
-#ifndef LEAST_SQUARE_SOLVER_HPP
-#define LEAST_SQUARE_SOLVER_HPP
-
-#include <algebra/SmallMatrix.hpp>
-#include <algebra/SmallVector.hpp>
-
-class LeastSquareSolver
-{
- public:
-  template <typename T>
-  void solveLocalSystem(const SmallMatrix<T>& A, SmallVector<T>& x, const SmallVector<T>& b);
-
-  LeastSquareSolver()  = default;
-  ~LeastSquareSolver() = default;
-};
-
-#endif   // LEAST_SQUARE_SOLVER_HPP
diff --git a/src/language/modules/CouplageModule.cpp b/src/language/modules/CouplageModule.cpp
index fc0553f8820f435bb84ed4543350e51bd521dbdc..3f53aa49a7bb067488df319e02abc650f9fffe18 100644
--- a/src/language/modules/CouplageModule.cpp
+++ b/src/language/modules/CouplageModule.cpp
@@ -1,3 +1,4 @@
+#include "mesh/DualMeshType.hpp"
 #include <utils/pugs_config.hpp>
 
 #ifdef PUGS_HAS_COSTO
@@ -10,10 +11,13 @@
 #include <memory>
 #include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 #include <mesh/MeshNodeInterface.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <utils/CommunicatorManager.hpp>
@@ -28,230 +32,150 @@ CouplageModule::CouplageModule()
 {
   this->_addBuiltinFunction("serialize",
                             std::function(
-                              [](std::shared_ptr<const IMesh> mesh,
+                              [](const std::shared_ptr<const MeshVariant> mesh_v,
                                  const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
-                                 const int64_t dest, const int64_t tag) -> void {
-                                /* const int tag = 20; */
-                                switch (mesh->dimension()) {
-                                case 1: {
-                                  return Serializer<1>().interfaces(mesh, i_interface_list, dest, tag);
-                                }
-                                case 2: {
-                                  return Serializer<2>().interfaces(mesh, i_interface_list, dest, tag);
-                                }
-                                case 3: {
-                                  return Serializer<3>().interfaces(mesh, i_interface_list, dest, tag);
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
+                                 const int64_t dest,
+				 const int64_t tag) -> void {
+				return Serializer().interfaces(mesh_v, i_interface_list, dest, tag);
                               }));
 
-  this->_addBuiltinFunction("serialize", std::function([](std::shared_ptr<const IMesh> mesh, const int64_t dest,
+  this->_addBuiltinFunction("serialize", std::function([](const std::shared_ptr<const MeshVariant> mesh_v,
+							  const int64_t dest,
                                                           const int64_t tag) -> void {
-                              switch (mesh->dimension()) {
-                              case 1: {
-                                return Serializer<1>().mesh(mesh, dest, tag);
-                              }
-                              case 2: {
-                                return Serializer<2>().mesh(mesh, dest, tag);
-                              }
-                              case 3: {
-                                return Serializer<3>().mesh(mesh, dest, tag);
-                              }
-                              default: {
-                                throw UnexpectedError("invalid mesh dimension");
-                              }
-                              }
+                                return Serializer().mesh(mesh_v, dest, tag);
+
                             }));
 
   this->_addBuiltinFunction("serialize", std::function([](const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-                                                          std::shared_ptr<const IMesh> mesh, int64_t location,
+							  const std::shared_ptr<const MeshVariant> mesh_v,
+							  const int64_t location,
                                                           const int64_t dest, const int64_t tag) -> void {
-                              std::cout << "\033[01;32m"
-                                        << "dest: " << dest << "\ntag: " << tag << "\033[00;00m" << std::endl;
-                              switch (mesh->dimension()) {
-                              case 1: {
-                                return Serializer<1>().BC(boundary, mesh, location, dest, tag);
-                              }
-                              case 2: {
-                                return Serializer<2>().BC(boundary, mesh, location, dest, tag);
-                              }
-                              case 3: {
-                                return Serializer<3>().BC(boundary, mesh, location, dest, tag);
-                              }
-                              default: {
-                                throw UnexpectedError("invalid mesh dimension");
-                              }
-                              }
-                            }));
-
-  this->_addBuiltinFunction("serialize_field",
-                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& field,
-                                             const int64_t dest, const int64_t tag) -> void {
-                              /* const int tag                = 300; */
-                              const std::shared_ptr i_mesh = getCommonMesh({field});
-                              if (not i_mesh) {
-                                throw NormalError("primal discrete functions are not defined on the same mesh");
-                              }
-
-                              switch (i_mesh->dimension()) {
-                              case 1: {
-                                return Serializer<1>().field(field, dest, tag);
-                              }
-                              case 2: {
-                                return Serializer<2>().field(field, dest, tag);
-                              }
-                              case 3: {
-                                return Serializer<3>().field(field, dest, tag);
-                              }
-                              default: {
-                                throw UnexpectedError("invalid mesh dimension");
-                              }
-                              }
+                                return Serializer().BC(boundary, mesh_v, location, dest, tag);
                             }));
 
   this->_addBuiltinFunction("serialize_field",
-                            std::function([](std::shared_ptr<const IMesh> mesh,
+                            std::function([](const std::shared_ptr<const MeshVariant> mesh,
                                              const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-                                             const std::shared_ptr<const ItemValueVariant>& field, const int64_t dest,
+                                             const std::shared_ptr<const ItemValueVariant>& field,
+					     const int64_t dest,
                                              const int64_t tag) -> void {
-                              /* const int tag = 300; */
-                              /* const std::shared_ptr i_mesh = getCommonMesh({mesh, field}); */
-                              /* if (not i_mesh) { */
-                              /*   throw NormalError("primal discrete functions are not defined on the same mesh"); */
-                              /* } */
-
-                              switch (mesh->dimension()) {
-                              case 1: {
-                                return Serializer<1>().BC_field(mesh, boundary, field, dest, tag);
-                              }
-                              case 2: {
-                                return Serializer<2>().BC_field(mesh, boundary, field, dest, tag);
-                              }
-                              case 3: {
-                                return Serializer<3>().BC_field(mesh, boundary, field, dest, tag);
-                              }
-                              default: {
-                                throw UnexpectedError("invalid mesh dimension");
-                              }
-                              }
+			      return Serializer().BC_field(mesh, boundary, field, dest, tag);
                             }));
 
-  this->_addBuiltinFunction("cpl_change_mesh_position",
-                            std::function([](std::shared_ptr<const IMesh> mesh, const int64_t src,
-                                             const int64_t tag) -> std::shared_ptr<const IMesh> {
-                              /* tag == 400; */
-                              /* src = 1 */
-                              switch (mesh->dimension()) {
-                              case 1: {
-                                return Serializer<1>().change_mesh_position(mesh, src, tag);
-                              }
-                              case 2: {
-                                return Serializer<2>().change_mesh_position(mesh, src, tag);
-                              }
-                              case 3: {
-                                return Serializer<3>().change_mesh_position(mesh, src, tag);
-                              }
-                              default: {
-                                throw UnexpectedError("invalid mesh dimension");
-                              }
-                              }
-                            }));
-
-  this->_addBuiltinFunction("get_flux_and_velocity",
-                            std::function(
-
-                              [](const std::shared_ptr<const IMesh>& i_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const std::shared_ptr<const ItemValueVariant>& ur,
-                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
-                                 const std::shared_ptr<const ItemValueVariant>& ur_saved,
-                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_saved)
-                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
-                                              std::shared_ptr<const SubItemValuePerItemVariant>,
-                                              std::shared_ptr<const ItemValueVariant>,
-                                              std::shared_ptr<const SubItemValuePerItemVariant>> {
-                                switch (i_mesh->dimension()) {
-                                case 1: {
-                                  auto myCouplingData = CouplingData<1>::getInstanceOrBuildIt();
-                                  return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
-                                                                              ur_saved, Fjr_saved);
-                                }
-                                case 2: {
-                                  auto myCouplingData = CouplingData<2>::getInstanceOrBuildIt();
-                                  return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
-                                                                              ur_saved, Fjr_saved);
-                                }
-                                case 3: {
-                                  auto myCouplingData = CouplingData<3>::getInstanceOrBuildIt();
-                                  return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
-                                                                              ur_saved, Fjr_saved);
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
-                                /* } */
-                                /* } */
-                              }));
-
-  this->_addBuiltinFunction("update_mesh", std::function(
-
-                                             [](const std::shared_ptr<const IMesh>& i_mesh,
-                                                const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-                                                const int64_t src, const int64_t tag) -> std::shared_ptr<const IMesh> {
-                                               switch (i_mesh->dimension()) {
-                                               case 1: {
-                                                 auto myCouplingData = CouplingData<1>::getInstanceOrBuildIt();
-                                                 return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
-                                               }
-                                               case 2: {
-                                                 auto myCouplingData = CouplingData<2>::getInstanceOrBuildIt();
-                                                 return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
-                                               }
-                                               case 3: {
-                                                 auto myCouplingData = CouplingData<3>::getInstanceOrBuildIt();
-                                                 return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
-                                               }
-                                               default: {
-                                                 throw UnexpectedError("invalid mesh dimension");
-                                               }
-                                               }
-                                             }));
-
-  this->_addBuiltinFunction("test", std::function(
-
-                                      [](const std::shared_ptr<const IMesh>& i_mesh) -> void {
-                                        auto dim = i_mesh->dimension();
-                                        std::cout << "in test dim: " << dim << std::endl;
-                                        using MeshType = Mesh<Connectivity<2>>;
-
-                                        std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-                                        auto nNodes            = p_mesh->numberOfNodes();
-                                        std::cout << "n nodes: " << nNodes << std::endl;
-                                        auto nCells = p_mesh->numberOfCells();
-                                        std::cout << "n cells: " << nCells << std::endl;
-                                        auto nEdges = p_mesh->numberOfEdges();
-                                        std::cout << "n edges: " << nEdges << std::endl;
-                                        auto nFaces = p_mesh->numberOfFaces();
-                                        std::cout << "n faces: " << nFaces << std::endl;
-                                        const auto nodes = p_mesh->xr();
-
-                                        /* if (CommunicatorManager::hasSplitColor()) { */
-                                        /* } */
-                                      }
-
-                                      ));
+  // this->_addBuiltinFunction("serialize_field",
+  //                           std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& field,
+  //                                            const int64_t dest, const int64_t tag) -> void {
+  //                             /* const int tag                = 300; */
+  //                             // const std::shared_ptr i_mesh = getCommonMesh({field});
+  //                             // if (not i_mesh) {
+  //                             //   throw NormalError("primal discrete functions are not defined on the same mesh");
+  //                             // }
+
+  //                             // switch (i_mesh->dimension()) {
+  //                             // case 1: {
+  //                             //   return Serializer<1>().field(field, dest, tag);
+  //                             // }
+  //                             // case 2: {
+  //                             //   return Serializer<2>().field(field, dest, tag);
+  //                             // }
+  //                             // case 3: {
+  //                             //   return Serializer<3>().field(field, dest, tag);
+  //                             // }
+  //                             // default: {
+  //                             //   throw UnexpectedError("invalid mesh dimension");
+  //                             // }
+  //                             // }
+  //                           }));
+
+  // this->_addBuiltinFunction("cpl_change_mesh_position",
+  //                           std::function([](std::shared_ptr<const MeshVariant> mesh, const int64_t src,
+  //                                            const int64_t tag) -> std::shared_ptr<const IMesh> {
+  //                             /* tag == 400; */
+  //                             /* src = 1 */
+  //                             // switch (mesh->dimension()) {
+  //                             // case 1: {
+  //                             //   return Serializer<1>().change_mesh_position(mesh, src, tag);
+  //                             // }
+  //                             // case 2: {
+  //                             //   return Serializer<2>().change_mesh_position(mesh, src, tag);
+  //                             // }
+  //                             // case 3: {
+  //                             //   return Serializer<3>().change_mesh_position(mesh, src, tag);
+  //                             // }
+  //                             // default: {
+  //                             //   throw UnexpectedError("invalid mesh dimension");
+  //                             // }
+  //                             // }
+  //                           }));
+
+  // this->_addBuiltinFunction("get_flux_and_velocity",
+  //                           std::function(
+
+  //                             [](std::shared_ptr<const MeshVariant> mesh,
+  //                                const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+  //                                  bc_descriptor_list,
+  //                                const std::shared_ptr<const ItemValueVariant>& ur,
+  //                                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+  //                                const std::shared_ptr<const ItemValueVariant>& ur_saved,
+  //                                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr_saved)
+  //                               -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+  //                                             std::shared_ptr<const SubItemValuePerItemVariant>,
+  //                                             std::shared_ptr<const ItemValueVariant>,
+  //                                             std::shared_ptr<const SubItemValuePerItemVariant>> {
+  //                               // switch (i_mesh->dimension()) {
+  //                               // case 1: {
+  //                               //   auto myCouplingData = CouplingData<1>::getInstanceOrBuildIt();
+  //                               //   return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
+  //                               //                                               ur_saved, Fjr_saved);
+  //                               // }
+  //                               // case 2: {
+  //                               //   auto myCouplingData = CouplingData<2>::getInstanceOrBuildIt();
+  //                               //   return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
+  //                               //                                               ur_saved, Fjr_saved);
+  //                               // }
+  //                               // case 3: {
+  //                               //   auto myCouplingData = CouplingData<3>::getInstanceOrBuildIt();
+  //                               //   return myCouplingData.get_flux_and_velocity(i_mesh, bc_descriptor_list, ur, Fjr,
+  //                               //                                               ur_saved, Fjr_saved);
+  //                               // }
+  //                               // default: {
+  //                               //   throw UnexpectedError("invalid mesh dimension");
+  //                               // }
+  //                               // }
+  //                               /* } */
+  //                               /* } */
+  //                             }));
+
+  // this->_addBuiltinFunction("update_mesh", std::function(
+
+  //                                            [](std::shared_ptr<const MeshVariant> mesh,
+  //                                               const std::shared_ptr<const IBoundaryDescriptor>& boundary,
+  //                                               const int64_t src, const int64_t tag) -> std::shared_ptr<const IMesh> {
+  //                                              // switch (i_mesh->dimension()) {
+  //                                              // case 1: {
+  //                                              //   auto myCouplingData = CouplingData<1>::getInstanceOrBuildIt();
+  //                                              //   return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
+  //                                              // }
+  //                                              // case 2: {
+  //                                              //   auto myCouplingData = CouplingData<2>::getInstanceOrBuildIt();
+  //                                              //   return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
+  //                                              // }
+  //                                              // case 3: {
+  //                                              //   auto myCouplingData = CouplingData<3>::getInstanceOrBuildIt();
+  //                                              //   return myCouplingData.update_mesh(i_mesh, boundary, src, tag);
+  //                                              // }
+  //                                              // default: {
+  //                                              //   throw UnexpectedError("invalid mesh dimension");
+  //                                              // }
+  //                                              // }
+  //                                            }));
 
   this->_addBuiltinFunction("hello", std::function(
 
                                        []() -> void { std::cout << "in test" << std::endl; }
 
                                        ));
-  // MathFunctionRegisterForVh{*this};
+
 }
 
 void
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 1d74740b2ebc51db875d1636832fff0368d8af40..3526c4d1d81c35fd7e2dd92b6c86867a09c83670 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -43,10 +43,8 @@
 #include <scheme/InflowBoundaryConditionDescriptor.hpp>
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
 #include <scheme/OutflowBoundaryConditionDescriptor.hpp>
-#include <scheme/ScalarDiamondScheme.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VariableBCDescriptor.hpp>
-#include <scheme/VectorDiamondScheme.hpp>
 #include <utils/CouplingData.hpp>
 #include <utils/Socket.hpp>
 
@@ -594,75 +592,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction(
-    "parabolicheat", std::function(
-
-                       [](const std::shared_ptr<const DiscreteFunctionVariant>& alpha,
-                          const std::shared_ptr<const DiscreteFunctionVariant>& mub_dual,
-                          const std::shared_ptr<const DiscreteFunctionVariant>& mu_dual,
-                          const std::shared_ptr<const DiscreteFunctionVariant>& f,
-                          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-                         -> std::shared_ptr<const DiscreteFunctionVariant> {
-                         return ScalarDiamondSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
-                       }
-
-                       ));
-
-  this->_addBuiltinFunction("unsteadyelasticity",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant> alpha,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> lambdab,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> mub,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> lambda,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> mu,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> f,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const ItemValueVariant>> {
-                                return VectorDiamondSchemeHandler{alpha, lambdab,           mub, lambda, mu,
-                                                                  f,     bc_descriptor_list}
-                                  .solution();
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("moleculardiffusion",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant> alpha,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> lambdab,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> mub,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> lambda,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> mu,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> f,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return VectorDiamondSchemeHandler{alpha, lambdab,           mub, lambda, mu,
-                                                                  f,     bc_descriptor_list}
-                                  .apply();
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("energybalance",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant> lambdab,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> mub,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> U,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> dual_U,
-                                 const std::shared_ptr<const DiscreteFunctionVariant> source,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return EnergyComputerHandler{lambdab, mub, U, dual_U, source, bc_descriptor_list}
-                                  .computeEnergyUpdate();
-                              }
-
-                              ));
-
   this->_addBuiltinFunction("hyperelastic_eucclhyd_fluxes",
                             std::function(
 
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 63334ab8b14f91e140ac78fac5588bfcca5e80a6..e39533dd16dbcdc1cec120d2002755559ea76b7a 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -688,9 +688,9 @@ AcousticSolverHandler::AcousticSolver<MeshType>::_applyExternalVelocityBC(const
   }
 }
 
-template <size_t Dimension>
+template <MeshConcept MeshType>
 void
-AcousticSolverHandler::AcousticSolver<Dimension>::_applyCouplingBC(const BoundaryConditionList& bc_list,
+AcousticSolverHandler::AcousticSolver<MeshType>::_applyCouplingBC(const BoundaryConditionList& bc_list,
                                                                    NodeValue<Rdxd>& Ar,
                                                                    NodeValue<Rd>& br) const
 {
@@ -956,11 +956,11 @@ class AcousticSolverHandler::AcousticSolver<MeshType>::SymmetryBoundaryCondition
   ~SymmetryBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class AcousticSolverHandler::AcousticSolver<Dimension>::CouplingBoundaryCondition
+template <MeshConcept MeshType>
+class AcousticSolverHandler::AcousticSolver<MeshType>::CouplingBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
  public:
   const Array<const NodeId>&
@@ -969,11 +969,9 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::CouplingBoundaryConditio
     return m_mesh_node_boundary.nodeList();
   }
 
-  CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary)
+  CouplingBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary)
     : m_mesh_node_boundary{mesh_node_boundary}
-  {
-    ;
-  }
+  {}
 
   ~CouplingBoundaryCondition() = default;
 };
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index eb9936dc371efeac7fc5df5d7591460e5511fc41..73f9522ba12ebb5f2570486a1bcb3a108b6b4270 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -82,4 +82,5 @@ class AcousticSolverHandler
 
   AcousticSolverHandler(const std::shared_ptr<const MeshVariant>& mesh_v);
 };
+
 #endif   // ACOUSTIC_SOLVER_HPP
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 194e33be93280a9ae3de5f6e76583c2036035a77..dadbe73bb24ad84e491248b89016baf85f0dbc31 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -9,8 +9,6 @@ add_library(
   DiscreteFunctionUtils.cpp
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
-  ScalarDiamondScheme.cpp
-  VectorDiamondScheme.cpp
   FluxingAdvectionSolver.cpp
 )
 
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index df662922a9f5db8584783f2ecd9848598721ad05..ea1b0f525beb3cea6f29e29c460610cda2290258 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -308,7 +308,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
           costo->recvData(recv, src, tag);
 
           if constexpr (Dimension == 1) {
-            MeshNodeBoundary<Dimension> mesh_node_boundary =
+            MeshNodeBoundary mesh_node_boundary =
               getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
             Array<double> face_values{mesh_node_boundary.nodeList().size()};
             for (size_t i_node = 0; i_node < mesh_node_boundary.nodeList().size(); ++i_node) {
@@ -319,7 +319,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
           } else {
             const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
 
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
+            MeshFaceBoundary mesh_face_boundary =
               getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
             Array<double> face_values{mesh_face_boundary.faceList().size()};
@@ -355,7 +355,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
 #ifdef PUGS_HAS_COSTO
         else if (dirichlet_bc_descriptor.name() == "cpl_forces") {
           auto costo = parallel::Messenger::getInstance().myCoupling;
-          MeshNodeBoundary<Dimension> mesh_node_boundary =
+          MeshNodeBoundary mesh_node_boundary =
             getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
           const int src = costo->myGlobalSize() - 1;
           const int tag = 200;
@@ -769,9 +769,9 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyNormalStressBC(co
   }
 }
 #ifdef PUGS_HAS_COSTO
-template <size_t Dimension>
+template <MeshConcept MeshType>
 void
-HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCplDirichletBC(const BoundaryConditionList& bc_list,
+HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyCplDirichletBC(const BoundaryConditionList& bc_list,
                                                                                NodeValue<Rd>& br) const
 {
   for (const auto& boundary_condition : bc_list) {
@@ -861,9 +861,9 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyVelocityBC(const
   }
 }
 
-template <size_t Dimension>
+template <MeshConcept MeshType>
 void
-HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC(const BoundaryConditionList& bc_list,
+HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyCouplingBC(const BoundaryConditionList& bc_list,
                                                                            NodeValue<Rdxd>& Ar,
                                                                            NodeValue<Rd>& br) const
 {
@@ -945,12 +945,12 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC(const
   }
 }
 
-template <size_t Dimension>
+template <MeshConcept MeshType>
 void
-HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC2(const MeshType& mesh,
-                                                                            const BoundaryConditionList& bc_list,
-                                                                            NodeValue<Rd>& ur,
-                                                                            NodeValuePerCell<Rd>& Fjr) const
+HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyCouplingBC2(const MeshType& mesh,
+									   const BoundaryConditionList& bc_list,
+									   NodeValue<Rd>& ur,
+									   NodeValuePerCell<Rd>& Fjr) const
 {
   for (const auto& boundary_condition : bc_list) {
     std::visit(
@@ -969,7 +969,7 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC2(cons
 
           auto costo          = parallel::Messenger::getInstance().myCoupling;
           auto rank           = costo->myGlobalRank();
-          auto myCouplingData = CouplingData<Dimension>::getInstanceOrBuildIt();
+          auto myCouplingData = CouplingData<MeshType>::getInstanceOrBuildIt();
 
           if (costo->m == 1) {
             /* std::cout << "\033[01;33m" */
@@ -1170,11 +1170,11 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::NormalStressBound
 };
 
 #ifdef PUGS_HAS_COSTO
-template <size_t Dimension>
-class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::CplDirichletBoundaryCondition
+template <MeshConcept MeshType>
+class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::CplDirichletBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
   const Array<const Rd> m_value_list;
 
  public:
@@ -1190,7 +1190,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::CplDirichletBoun
     return m_value_list;
   }
 
-  CplDirichletBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,
+  CplDirichletBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
                                 const Array<const Rd>& value_list)
     : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
   {}
@@ -1292,11 +1292,11 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::SymmetryBoundaryC
   ~SymmetryBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
-class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::CouplingBoundaryCondition
+template <MeshConcept MeshType>
+class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::CouplingBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
  public:
   const Array<const NodeId>&
@@ -1305,7 +1305,7 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::CouplingBoundary
     return m_mesh_node_boundary.nodeList();
   }
 
-  CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary)
+  CouplingBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary)
     : m_mesh_node_boundary{mesh_node_boundary}
   {
     ;
diff --git a/src/scheme/ScalarDiamondScheme.cpp b/src/scheme/ScalarDiamondScheme.cpp
deleted file mode 100644
index d0a9f7e48f7f93b9224a197312dbb2a0b2829847..0000000000000000000000000000000000000000
--- a/src/scheme/ScalarDiamondScheme.cpp
+++ /dev/null
@@ -1,844 +0,0 @@
-#include <scheme/ScalarDiamondScheme.hpp>
-
-#include <scheme/DiscreteFunctionP0.hpp>
-#include <scheme/DiscreteFunctionUtils.hpp>
-
-template <size_t Dimension>
-class ScalarDiamondSchemeHandler::InterpolationWeightsManager
-{
- private:
-  std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh;
-  FaceValue<bool> m_primal_face_is_on_boundary;
-  NodeValue<bool> m_primal_node_is_on_boundary;
-  CellValuePerNode<double> m_w_rj;
-  FaceValuePerNode<double> m_w_rl;
-
- public:
-  CellValuePerNode<double>&
-  wrj()
-  {
-    return m_w_rj;
-  }
-
-  FaceValuePerNode<double>&
-  wrl()
-  {
-    return m_w_rl;
-  }
-
-  void
-  compute()
-  {
-    using MeshDataType      = MeshData<Dimension>;
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
-
-    const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
-    const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-    const auto& node_to_cell_matrix                  = m_mesh->connectivity().nodeToCellMatrix();
-    const auto& node_to_face_matrix                  = m_mesh->connectivity().nodeToFaceMatrix();
-    CellValuePerNode<double> w_rj{m_mesh->connectivity()};
-    FaceValuePerNode<double> w_rl{m_mesh->connectivity()};
-
-    for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
-      w_rl[i] = std::numeric_limits<double>::signaling_NaN();
-    }
-
-    for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
-      SmallVector<double> b{Dimension + 1};
-      b[0] = 1;
-      for (size_t i = 1; i < Dimension + 1; i++) {
-        b[i] = xr[i_node][i - 1];
-      }
-      const auto& node_to_cell = node_to_cell_matrix[i_node];
-
-      if (not m_primal_node_is_on_boundary[i_node]) {
-        SmallMatrix<double> A{Dimension + 1, node_to_cell.size()};
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          A(0, j) = 1;
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J = node_to_cell[j];
-            A(i, j)        = xj[J][i - 1];
-          }
-        }
-        SmallVector<double> x{node_to_cell.size()};
-        x = zero;
-
-        LeastSquareSolver ls_solver;
-        ls_solver.solveLocalSystem(A, x, b);
-
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          w_rj(i_node, j) = x[j];
-        }
-      } else {
-        int nb_face_used = 0;
-        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[i_node][i_face];
-          if (m_primal_face_is_on_boundary[face_id]) {
-            nb_face_used++;
-          }
-        }
-        SmallMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
-        for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
-          A(0, j) = 1;
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J = node_to_cell[j];
-            A(i, j)        = xj[J][i - 1];
-          }
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          int cpt_face = 0;
-          for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-            FaceId face_id = node_to_face_matrix[i_node][i_face];
-            if (m_primal_face_is_on_boundary[face_id]) {
-              A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
-              cpt_face++;
-            }
-          }
-        }
-
-        SmallVector<double> x{node_to_cell.size() + nb_face_used};
-        x = zero;
-
-        LeastSquareSolver ls_solver;
-        ls_solver.solveLocalSystem(A, x, b);
-
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          w_rj(i_node, j) = x[j];
-        }
-        int cpt_face = node_to_cell.size();
-        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[i_node][i_face];
-          if (m_primal_face_is_on_boundary[face_id]) {
-            w_rl(i_node, i_face) = x[cpt_face++];
-          }
-        }
-      }
-    }
-    m_w_rj = w_rj;
-    m_w_rl = w_rl;
-  }
-
-  InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh,
-                              FaceValue<bool> primal_face_is_on_boundary,
-                              NodeValue<bool> primal_node_is_on_boundary)
-    : m_mesh(mesh),
-      m_primal_face_is_on_boundary(primal_face_is_on_boundary),
-      m_primal_node_is_on_boundary(primal_node_is_on_boundary)
-  {}
-
-  ~InterpolationWeightsManager() = default;
-};
-
-class ScalarDiamondSchemeHandler::IScalarDiamondScheme
-{
- public:
-  virtual std::shared_ptr<const DiscreteFunctionVariant> getSolution() const = 0;
-
-  IScalarDiamondScheme()          = default;
-  virtual ~IScalarDiamondScheme() = default;
-};
-
-template <size_t Dimension>
-class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSchemeHandler::IScalarDiamondScheme
-{
- private:
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshData<Dimension>;
-
-  DiscreteFunctionP0<Dimension, double> m_solution;
-
-  class DirichletBoundaryCondition
-  {
-   private:
-    const Array<const double> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const double>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~DirichletBoundaryCondition() = default;
-  };
-
-  class NeumannBoundaryCondition
-  {
-   private:
-    const Array<const double> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const double>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    NeumannBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~NeumannBoundaryCondition() = default;
-  };
-
-  class FourierBoundaryCondition
-  {
-   private:
-    const Array<const double> m_coef_list;
-    const Array<const double> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const double>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    const Array<const double>&
-    coefList() const
-    {
-      return m_coef_list;
-    }
-
-   public:
-    FourierBoundaryCondition(const Array<const FaceId>& face_list,
-                             const Array<const double>& coef_list,
-                             const Array<const double>& value_list)
-      : m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_coef_list.size() == m_face_list.size());
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~FourierBoundaryCondition() = default;
-  };
-
-  class SymmetryBoundaryCondition
-  {
-   private:
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-   public:
-    SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
-
-    ~SymmetryBoundaryCondition() = default;
-  };
-
- public:
-  std::shared_ptr<const DiscreteFunctionVariant>
-  getSolution() const final
-  {
-    return std::make_shared<DiscreteFunctionVariant>(m_solution);
-  }
-
-  ScalarDiamondScheme(const std::shared_ptr<const MeshType>& mesh,
-                      const DiscreteFunctionP0<Dimension, const double>& alpha,
-                      const DiscreteFunctionP0<Dimension, const double>& dual_mub,
-                      const DiscreteFunctionP0<Dimension, const double>& dual_mu,
-                      const DiscreteFunctionP0<Dimension, const double>& f,
-                      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-    : m_solution(mesh)
-  {
-    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mu.mesh(),
-           "diffusion coefficient is not defined on the dual mesh!");
-    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mub.mesh(),
-           "boundary diffusion coefficient is not defined on the dual mesh!");
-
-    using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition,
-                                           NeumannBoundaryCondition, SymmetryBoundaryCondition>;
-
-    using BoundaryConditionList = std::vector<BoundaryCondition>;
-
-    BoundaryConditionList boundary_condition_list;
-
-    for (const auto& bc_descriptor : bc_descriptor_list) {
-      bool is_valid_boundary_condition = true;
-
-      switch (bc_descriptor->type()) {
-      case IBoundaryConditionDescriptor::Type::symmetry: {
-        throw NotImplementedError("NIY");
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::dirichlet: {
-        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
-          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
-        if constexpr (Dimension > 1) {
-          MeshFaceBoundary<Dimension> mesh_face_boundary =
-            getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-          const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-          MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
-
-          Array<const double> value_list =
-            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
-                                                                                                      mesh_data.xl(),
-                                                                                                      mesh_face_boundary
-                                                                                                        .faceList());
-
-          boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-
-        } else {
-          throw NotImplementedError("Dirichlet BC in 1d");
-        }
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::fourier: {
-        throw NotImplementedError("NIY");
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::neumann: {
-        const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
-          dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor);
-
-        if constexpr (Dimension > 1) {
-          MeshFaceBoundary<Dimension> mesh_face_boundary =
-            getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
-
-          const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
-          MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
-
-          Array<const double> value_list =
-            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
-                                                                                                      mesh_data.xl(),
-                                                                                                      mesh_face_boundary
-                                                                                                        .faceList());
-
-          boundary_condition_list.push_back(NeumannBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-
-        } else {
-          throw NotImplementedError("Dirichlet BC in 1d");
-        }
-        break;
-      }
-      default: {
-        is_valid_boundary_condition = false;
-      }
-      }
-      if (not is_valid_boundary_condition) {
-        std::ostringstream error_msg;
-        error_msg << *bc_descriptor << " is an invalid boundary condition for heat equation";
-        throw NormalError(error_msg.str());
-      }
-    }
-
-    if constexpr (Dimension > 1) {
-      const CellValue<const size_t> cell_dof_number = [&] {
-        CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
-        return compute_cell_dof_number;
-      }();
-      size_t number_of_dof = mesh->numberOfCells();
-
-      const FaceValue<const size_t> face_dof_number = [&] {
-        FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
-        compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
-                            (std::is_same_v<T, FourierBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>) or
-                            (std::is_same_v<T, DirichletBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id = face_list[i_face];
-                  if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
-                    std::ostringstream os;
-                    os << "The face " << face_id << " is used at least twice for boundary conditions";
-                    throw NormalError(os.str());
-                  } else {
-                    compute_face_dof_number[face_id] = number_of_dof++;
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return compute_face_dof_number;
-      }();
-
-      const FaceValue<const bool> primal_face_is_neumann = [&] {
-        FaceValue<bool> face_is_neumann{mesh->connectivity()};
-        face_is_neumann.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
-                            (std::is_same_v<T, FourierBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id     = face_list[i_face];
-                  face_is_neumann[face_id] = true;
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return face_is_neumann;
-      }();
-
-      const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
-      const auto& face_to_cell_matrix        = mesh->connectivity().faceToCellMatrix();
-      NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
-      }
-
-      primal_node_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-            NodeId node_id                      = primal_face_to_node_matrix[face_id][i_node];
-            primal_node_is_on_boundary[node_id] = true;
-          }
-        }
-      }
-
-      FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
-      }
-
-      primal_face_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          primal_face_is_on_boundary[face_id] = true;
-        }
-      }
-
-      FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_neumann is incorrect");
-      }
-
-      primal_face_is_dirichlet.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
-      }
-
-      InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary);
-      iwm.compute();
-      CellValuePerNode<double> w_rj = iwm.wrj();
-      FaceValuePerNode<double> w_rl = iwm.wrl();
-
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-      const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
-      const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-      const auto& node_to_face_matrix                  = mesh->connectivity().nodeToFaceMatrix();
-
-      {
-        std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
-
-        MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
-
-        std::shared_ptr mapper =
-          DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
-
-        CellValue<const double> dual_kappaj  = dual_mu.cellValues();
-        CellValue<const double> dual_kappajb = dual_mub.cellValues();
-
-        const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
-
-        const FaceValue<const double> mes_l = [&] {
-          if constexpr (Dimension == 1) {
-            FaceValue<double> compute_mes_l{mesh->connectivity()};
-            compute_mes_l.fill(1);
-            return compute_mes_l;
-          } else {
-            return mesh_data.ll();
-          }
-        }();
-
-        const CellValue<const double> dual_mes_l_j = [=] {
-          CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
-          mapper->toDualCell(mes_l, compute_mes_j);
-
-          return compute_mes_j;
-        }();
-
-        FaceValue<const CellId> face_dual_cell_id = [=]() {
-          FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
-          CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
-
-          mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
-
-          return computed_face_dual_cell_id;
-        }();
-
-        NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
-          CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
-          cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> node_primal_id{mesh->connectivity()};
-
-          parallel_for(
-            mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
-
-          NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
-
-          mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
-
-          return computed_dual_node_primal_node_id;
-        }();
-
-        CellValue<NodeId> primal_cell_dual_node_id = [=]() {
-          CellValue<NodeId> cell_id{mesh->connectivity()};
-          NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
-          node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
-
-          CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
-
-          mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
-
-          return cell_id;
-        }();
-        const auto& dual_Cjr                     = diamond_mesh_data.Cjr();
-        FaceValue<TinyVector<Dimension>> dualClj = [&] {
-          FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
-          const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
-          const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
-          parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-              const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-              for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
-                CellId cell_id            = primal_face_to_cell[i];
-                const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
-                for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
-                     i_dual_cell++) {
-                  const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
-                  if (face_dual_cell_id[face_id] == dual_cell_id) {
-                    for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
-                         i_dual_node++) {
-                      const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                      if (final_dual_node_id == dual_node_id) {
-                        computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
-                      }
-                    }
-                  }
-                }
-              }
-            });
-          return computedClj;
-        }();
-
-        FaceValue<TinyVector<Dimension>> nlj = [&] {
-          FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
-          parallel_for(
-            mesh->numberOfFaces(),
-            PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
-          return computedNlj;
-        }();
-
-        FaceValue<const double> alpha_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        FaceValue<const double> alphab_l = [&] {
-          CellValue<double> alpha_jb{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_jb[diamond_cell_id] = dual_kappajb[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_lb{mesh->connectivity()};
-          mapper->fromDualCell(alpha_jb, computed_alpha_lb);
-          return computed_alpha_lb;
-        }();
-
-        const CellValue<const double> primal_Vj = mesh_data.Vj();
-
-        const Array<int> non_zeros{number_of_dof};
-        non_zeros.fill(Dimension);
-        CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zeros);
-
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-          const double beta_l             = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
-          const double betab_l            = l2Norm(dualClj[face_id]) * alphab_l[face_id] * mes_l[face_id];
-          for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
-            const CellId cell_id1 = primal_face_to_cell[i_cell];
-            for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
-              const CellId cell_id2 = primal_face_to_cell[j_cell];
-              if (i_cell == j_cell) {
-                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
-                if (primal_face_is_neumann[face_id]) {
-                  S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= betab_l;
-                }
-              } else {
-                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
-              }
-            }
-          }
-        }
-
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          const size_t j = cell_dof_number[cell_id];
-          S(j, j) += alpha[cell_id] * primal_Vj[cell_id];
-        }
-
-        const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
-
-        const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const double alpha_face_id  = mes_l[face_id] * alpha_l[face_id];
-          const double alphab_face_id = mes_l[face_id] * alphab_l[face_id];
-
-          for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
-            CellId i_id                            = face_to_cell_matrix[face_id][i_face_cell];
-            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
-
-            for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-              NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
-
-              const TinyVector<Dimension> nil = [&] {
-                if (is_face_reversed_for_cell_i) {
-                  return -nlj[face_id];
-                } else {
-                  return nlj[face_id];
-                }
-              }();
-
-              CellId dual_cell_id = face_dual_cell_id[face_id];
-
-              for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
-                const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                if (dual_node_primal_node_id[dual_node_id] == node_id) {
-                  const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-
-                  const double a_ir  = alpha_face_id * dot(nil, Clr);
-                  const double ab_ir = alphab_face_id * dot(nil, Clr);
-
-                  for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
-                    CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
-                    S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
-                    if (primal_face_is_neumann[face_id]) {
-                      S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * ab_ir;
-                    }
-                  }
-                  if (primal_node_is_on_boundary[node_id]) {
-                    for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
-                      FaceId l_id = node_to_face_matrix[node_id][l_face];
-                      if (primal_face_is_on_boundary[l_id]) {
-                        S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
-                        if (primal_face_is_neumann[face_id]) {
-                          S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * ab_ir;
-                        }
-                      }
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          if (primal_face_is_dirichlet[face_id]) {
-            S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
-          }
-        }
-
-        CellValue<const double> fj = f.cellValues();
-
-        CRSMatrix A{S.getCRSMatrix()};
-        Vector<double> b{number_of_dof};
-        b = zero;
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
-        }
-        // Dirichlet on b^L_D
-        {
-          for (const auto& boundary_condition : boundary_condition_list) {
-            std::visit(
-              [&](auto&& bc) {
-                using T = std::decay_t<decltype(bc)>;
-                if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
-                  const auto& face_list  = bc.faceList();
-                  const auto& value_list = bc.valueList();
-                  for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                    const FaceId face_id = face_list[i_face];
-                    b[face_dof_number[face_id]] += value_list[i_face];
-                  }
-                }
-              },
-              boundary_condition);
-          }
-        }
-        // EL b^L
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
-                            (std::is_same_v<T, FourierBoundaryCondition>)) {
-                const auto& face_list  = bc.faceList();
-                const auto& value_list = bc.valueList();
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  FaceId face_id = face_list[i_face];
-                  b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face];
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        Vector<double> T{number_of_dof};
-        T = zero;
-
-        LinearSolver solver;
-        solver.solveLocalSystem(A, T, b);
-
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { m_solution[cell_id] = T[cell_dof_number[cell_id]]; });
-      }
-    }
-  }
-};
-
-std::shared_ptr<const DiscreteFunctionVariant>
-ScalarDiamondSchemeHandler::solution() const
-{
-  return m_scheme->getSolution();
-}
-
-ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
-  const std::shared_ptr<const DiscreteFunctionVariant>& alpha,
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_mub,
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_mu,
-  const std::shared_ptr<const DiscreteFunctionVariant>& f,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-{
-  const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
-  if (not i_mesh) {
-    throw NormalError("primal discrete functions are not defined on the same mesh");
-  }
-  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_mub, dual_mu});
-  if (not i_dual_mesh) {
-    throw NormalError("dual discrete functions are not defined on the same mesh");
-  }
-  checkDiscretizationType({alpha, dual_mub, dual_mu, f}, DiscreteFunctionType::P0);
-
-  switch (i_mesh->dimension()) {
-  case 1: {
-    using MeshType                   = Mesh<Connectivity<1>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<1, const double>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_scheme = std::make_unique<ScalarDiamondScheme<1>>(mesh, alpha->get<DiscreteScalarFunctionType>(),
-                                                        dual_mub->get<DiscreteScalarFunctionType>(),
-                                                        dual_mu->get<DiscreteScalarFunctionType>(),
-                                                        f->get<DiscreteScalarFunctionType>(), bc_descriptor_list);
-    break;
-  }
-  case 2: {
-    using MeshType                   = Mesh<Connectivity<2>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<2, const double>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_scheme = std::make_unique<ScalarDiamondScheme<2>>(mesh, alpha->get<DiscreteScalarFunctionType>(),
-                                                        dual_mub->get<DiscreteScalarFunctionType>(),
-                                                        dual_mu->get<DiscreteScalarFunctionType>(),
-                                                        f->get<DiscreteScalarFunctionType>(), bc_descriptor_list);
-    break;
-  }
-  case 3: {
-    using MeshType                   = Mesh<Connectivity<3>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<3, const double>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_scheme = std::make_unique<ScalarDiamondScheme<3>>(mesh, alpha->get<DiscreteScalarFunctionType>(),
-                                                        dual_mub->get<DiscreteScalarFunctionType>(),
-                                                        dual_mu->get<DiscreteScalarFunctionType>(),
-                                                        f->get<DiscreteScalarFunctionType>(), bc_descriptor_list);
-    break;
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
-ScalarDiamondSchemeHandler::~ScalarDiamondSchemeHandler() = default;
diff --git a/src/scheme/ScalarDiamondScheme.hpp b/src/scheme/ScalarDiamondScheme.hpp
deleted file mode 100644
index 4040fe1ade92e3b599dcdf724e924057feacaf63..0000000000000000000000000000000000000000
--- a/src/scheme/ScalarDiamondScheme.hpp
+++ /dev/null
@@ -1,688 +0,0 @@
-#ifndef SCALAR_DIAMOND_SCHEME_HPP
-#define SCALAR_DIAMOND_SCHEME_HPP
-
-#include <algebra/CRSMatrix.hpp>
-#include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LinearSolver.hpp>
-#include <algebra/TinyVector.hpp>
-#include <algebra/Vector.hpp>
-#include <language/utils/InterpolateItemValue.hpp>
-#include <mesh/Connectivity.hpp>
-#include <mesh/DualConnectivityManager.hpp>
-#include <mesh/DualMeshManager.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
-#include <mesh/MeshFaceBoundary.hpp>
-#include <mesh/MeshNodeBoundary.hpp>
-#include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
-#include <mesh/SubItemValuePerItem.hpp>
-#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
-#include <scheme/DiscreteFunctionVariant.hpp>
-#include <scheme/FourierBoundaryConditionDescriptor.hpp>
-#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-
-class ScalarDiamondSchemeHandler
-{
- private:
-  class IScalarDiamondScheme;
-
-  template <size_t Dimension>
-  class ScalarDiamondScheme;
-
-  template <size_t Dimension>
-  class InterpolationWeightsManager;
-
- public:
-  std::unique_ptr<IScalarDiamondScheme> m_scheme;
-
-  std::shared_ptr<const DiscreteFunctionVariant> solution() const;
-
-  ScalarDiamondSchemeHandler(
-    const std::shared_ptr<const DiscreteFunctionVariant>& alpha,
-    const std::shared_ptr<const DiscreteFunctionVariant>& mu_dualb,
-    const std::shared_ptr<const DiscreteFunctionVariant>& mu_dual,
-    const std::shared_ptr<const DiscreteFunctionVariant>& f,
-    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
-
-  ~ScalarDiamondSchemeHandler();
-};
-
-template <size_t Dimension>
-class LegacyScalarDiamondScheme
-{
- private:
-  class DirichletBoundaryCondition
-  {
-   private:
-    const Array<const double> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const double>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~DirichletBoundaryCondition() = default;
-  };
-
-  class NeumannBoundaryCondition
-  {
-   private:
-    const Array<const double> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const double>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    NeumannBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~NeumannBoundaryCondition() = default;
-  };
-
-  class FourierBoundaryCondition
-  {
-   private:
-    const Array<const double> m_coef_list;
-    const Array<const double> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const double>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    const Array<const double>&
-    coefList() const
-    {
-      return m_coef_list;
-    }
-
-   public:
-    FourierBoundaryCondition(const Array<const FaceId>& face_list,
-                             const Array<const double>& coef_list,
-                             const Array<const double>& value_list)
-      : m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_coef_list.size() == m_face_list.size());
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~FourierBoundaryCondition() = default;
-  };
-
-  class SymmetryBoundaryCondition
-  {
-   private:
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-   public:
-    SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
-
-    ~SymmetryBoundaryCondition() = default;
-  };
-
- public:
-  LegacyScalarDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
-                            const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-                            const FunctionSymbolId& kappa_id,
-                            const FunctionSymbolId& f_id,
-                            CellValue<double>& Temperature,
-                            FaceValue<double>& Temperature_face,
-                            const double& Tf,
-                            const double& dt,
-                            const CellValuePerNode<double>& w_rj,
-                            const FaceValuePerNode<double>& w_rl)
-  {
-    using ConnectivityType = Connectivity<Dimension>;
-    using MeshType         = Mesh<ConnectivityType>;
-    using MeshDataType     = MeshData<Dimension>;
-
-    using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition,
-                                           NeumannBoundaryCondition, SymmetryBoundaryCondition>;
-
-    using BoundaryConditionList = std::vector<BoundaryCondition>;
-
-    BoundaryConditionList boundary_condition_list;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    for (const auto& bc_descriptor : bc_descriptor_list) {
-      bool is_valid_boundary_condition = true;
-
-      switch (bc_descriptor->type()) {
-      case IBoundaryConditionDescriptor::Type::symmetry: {
-        throw NotImplementedError("NIY");
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::dirichlet: {
-        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
-          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
-        if constexpr (Dimension > 1) {
-          MeshFaceBoundary<Dimension> mesh_face_boundary =
-            getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-          const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-          MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
-
-          Array<const double> value_list =
-            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
-                                                                                                      mesh_data.xl(),
-                                                                                                      mesh_face_boundary
-                                                                                                        .faceList());
-          boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-        }
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::fourier: {
-        throw NotImplementedError("NIY");
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::neumann: {
-        const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
-          dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor);
-
-        if constexpr (Dimension > 1) {
-          MeshFaceBoundary<Dimension> mesh_face_boundary =
-            getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
-
-          const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
-          MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
-
-          Array<const double> value_list =
-            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
-                                                                                                      mesh_data.xl(),
-                                                                                                      mesh_face_boundary
-                                                                                                        .faceList());
-
-          boundary_condition_list.push_back(NeumannBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-
-        } else {
-          throw NotImplementedError("Dirichlet BC in 1d");
-        }
-        break;
-      }
-      default: {
-        is_valid_boundary_condition = false;
-      }
-      }
-      if (not is_valid_boundary_condition) {
-        std::ostringstream error_msg;
-        error_msg << *bc_descriptor << " is an invalid boundary condition for heat equation";
-        throw NormalError(error_msg.str());
-      }
-    }
-
-    if constexpr (Dimension > 1) {
-      const CellValue<const size_t> cell_dof_number = [&] {
-        CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
-        return compute_cell_dof_number;
-      }();
-      size_t number_of_dof = mesh->numberOfCells();
-
-      const FaceValue<const size_t> face_dof_number = [&] {
-        FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
-        compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
-                            (std::is_same_v<T, FourierBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>) or
-                            (std::is_same_v<T, DirichletBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id = face_list[i_face];
-                  if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
-                    std::ostringstream os;
-                    os << "The face " << face_id << " is used at least twice for boundary conditions";
-                    throw NormalError(os.str());
-                  } else {
-                    compute_face_dof_number[face_id] = number_of_dof++;
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return compute_face_dof_number;
-      }();
-
-      const FaceValue<const bool> primal_face_is_neumann = [&] {
-        FaceValue<bool> face_is_neumann{mesh->connectivity()};
-        face_is_neumann.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
-                            (std::is_same_v<T, FourierBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id     = face_list[i_face];
-                  face_is_neumann[face_id] = true;
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return face_is_neumann;
-      }();
-
-      const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
-      const auto& face_to_cell_matrix        = mesh->connectivity().faceToCellMatrix();
-      NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
-      }
-
-      primal_node_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-            NodeId node_id                      = primal_face_to_node_matrix[face_id][i_node];
-            primal_node_is_on_boundary[node_id] = true;
-          }
-        }
-      }
-
-      FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
-      }
-
-      primal_face_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          primal_face_is_on_boundary[face_id] = true;
-        }
-      }
-
-      FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_neumann is incorrect");
-      }
-
-      primal_face_is_dirichlet.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
-      }
-
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-      const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
-      const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-      const auto& node_to_face_matrix                  = mesh->connectivity().nodeToFaceMatrix();
-
-      {
-        std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
-
-        MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
-
-        std::shared_ptr mapper =
-          DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
-
-        CellValue<double> kappaj =
-          InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
-                                                                                                    mesh_data.xj());
-
-        CellValue<double> dual_kappaj =
-          InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
-                                                                                                    diamond_mesh_data
-                                                                                                      .xj());
-
-        const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
-
-        const FaceValue<const double> mes_l = [&] {
-          if constexpr (Dimension == 1) {
-            FaceValue<double> compute_mes_l{mesh->connectivity()};
-            compute_mes_l.fill(1);
-            return compute_mes_l;
-          } else {
-            return mesh_data.ll();
-          }
-        }();
-
-        const CellValue<const double> dual_mes_l_j = [=] {
-          CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
-          mapper->toDualCell(mes_l, compute_mes_j);
-
-          return compute_mes_j;
-        }();
-
-        FaceValue<const CellId> face_dual_cell_id = [=]() {
-          FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
-          CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
-
-          mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
-
-          return computed_face_dual_cell_id;
-        }();
-
-        NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
-          CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
-          cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> node_primal_id{mesh->connectivity()};
-
-          parallel_for(
-            mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
-
-          NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
-
-          mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
-
-          return computed_dual_node_primal_node_id;
-        }();
-
-        CellValue<NodeId> primal_cell_dual_node_id = [=]() {
-          CellValue<NodeId> cell_id{mesh->connectivity()};
-          NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
-          node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
-
-          CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
-
-          mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
-
-          return cell_id;
-        }();
-        const auto& dual_Cjr                     = diamond_mesh_data.Cjr();
-        FaceValue<TinyVector<Dimension>> dualClj = [&] {
-          FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
-          const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
-          const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
-          parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-              const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-              for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
-                CellId cell_id            = primal_face_to_cell[i];
-                const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
-                for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
-                     i_dual_cell++) {
-                  const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
-                  if (face_dual_cell_id[face_id] == dual_cell_id) {
-                    for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
-                         i_dual_node++) {
-                      const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                      if (final_dual_node_id == dual_node_id) {
-                        computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
-                      }
-                    }
-                  }
-                }
-              }
-            });
-          return computedClj;
-        }();
-
-        FaceValue<TinyVector<Dimension>> nlj = [&] {
-          FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
-          parallel_for(
-            mesh->numberOfFaces(),
-            PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
-          return computedNlj;
-        }();
-
-        FaceValue<const double> alpha_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        double lambda      = (Tf == 0) ? 0 : 1;
-        double time_factor = (lambda == 0) ? 1 : dt;
-
-        const CellValue<const double> primal_Vj = mesh_data.Vj();
-        Array<int> non_zero{number_of_dof};
-        non_zero.fill(Dimension * Dimension);
-        CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zero);
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-          const double beta_l             = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
-          for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
-            const CellId cell_id1 = primal_face_to_cell[i_cell];
-            for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
-              const CellId cell_id2 = primal_face_to_cell[j_cell];
-              if (i_cell == j_cell) {
-                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) +=
-                  (time_factor * beta_l + lambda * primal_Vj[cell_id1]);
-                if (primal_face_is_neumann[face_id]) {
-                  S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= beta_l;
-                }
-              } else {
-                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= time_factor * beta_l;
-              }
-            }
-          }
-        }
-
-        const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
-
-        const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const double alpha_face_id = mes_l[face_id] * alpha_l[face_id];
-
-          for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
-            CellId i_id                            = face_to_cell_matrix[face_id][i_face_cell];
-            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
-
-            for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-              NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
-
-              const TinyVector<Dimension> nil = [&] {
-                if (is_face_reversed_for_cell_i) {
-                  return -nlj[face_id];
-                } else {
-                  return nlj[face_id];
-                }
-              }();
-
-              CellId dual_cell_id = face_dual_cell_id[face_id];
-
-              for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
-                const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                if (dual_node_primal_node_id[dual_node_id] == node_id) {
-                  const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-
-                  const double a_ir = alpha_face_id * dot(nil, Clr);
-
-                  for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
-                    CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
-                    S(cell_dof_number[i_id], cell_dof_number[j_id]) -= time_factor * w_rj(node_id, j_cell) * a_ir;
-                    if (primal_face_is_neumann[face_id]) {
-                      S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * a_ir;
-                    }
-                  }
-                  if (primal_node_is_on_boundary[node_id]) {
-                    for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
-                      FaceId l_id = node_to_face_matrix[node_id][l_face];
-                      if (primal_face_is_on_boundary[l_id]) {
-                        S(cell_dof_number[i_id], face_dof_number[l_id]) -= time_factor * w_rl(node_id, l_face) * a_ir;
-                        if (primal_face_is_neumann[face_id]) {
-                          S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * a_ir;
-                        }
-                      }
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          if (primal_face_is_dirichlet[face_id]) {
-            S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
-          }
-        }
-
-        CellValue<double> fj =
-          InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id,
-                                                                                                    mesh_data.xj());
-
-        CRSMatrix A{S.getCRSMatrix()};
-        Vector<double> b{number_of_dof};
-        b = zero;
-
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          b[cell_dof_number[cell_id]] =
-            (time_factor * fj[cell_id] + lambda * Temperature[cell_id]) * primal_Vj[cell_id];
-        }
-        // Dirichlet on b^L_D
-        {
-          for (const auto& boundary_condition : boundary_condition_list) {
-            std::visit(
-              [&](auto&& bc) {
-                using T = std::decay_t<decltype(bc)>;
-                if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
-                  const auto& face_list  = bc.faceList();
-                  const auto& value_list = bc.valueList();
-                  for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                    const FaceId face_id = face_list[i_face];
-                    b[face_dof_number[face_id]] += value_list[i_face];
-                  }
-                }
-              },
-              boundary_condition);
-          }
-        }
-        // EL b^L
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
-                            (std::is_same_v<T, FourierBoundaryCondition>)) {
-                const auto& face_list  = bc.faceList();
-                const auto& value_list = bc.valueList();
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  FaceId face_id = face_list[i_face];
-                  b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face];
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        Vector<double> T{number_of_dof};
-        T = zero;
-
-        LinearSolver solver;
-        solver.solveLocalSystem(A, T, b);
-
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_dof_number[cell_id]]; });
-        parallel_for(
-          mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-            if (primal_face_is_neumann[face_id]) {
-              Temperature_face[face_id] = T[face_dof_number[face_id]];
-            }
-          });
-      }
-    }
-  }
-};
-
-template LegacyScalarDiamondScheme<1>::LegacyScalarDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  CellValue<double>&,
-  FaceValue<double>&,
-  const double&,
-  const double&,
-  const CellValuePerNode<double>& w_rj,
-  const FaceValuePerNode<double>& w_rl);
-
-template LegacyScalarDiamondScheme<2>::LegacyScalarDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  CellValue<double>&,
-  FaceValue<double>&,
-  const double&,
-  const double&,
-  const CellValuePerNode<double>& w_rj,
-  const FaceValuePerNode<double>& w_rl);
-
-template LegacyScalarDiamondScheme<3>::LegacyScalarDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  CellValue<double>&,
-  FaceValue<double>&,
-  const double&,
-  const double&,
-  const CellValuePerNode<double>& w_rj,
-  const FaceValuePerNode<double>& w_rl);
-
-#endif   // SCALAR_DIAMOND_SCHEME_HPP
diff --git a/src/scheme/SymmetryBoundaryConditionDescriptor.hpp b/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
index 68ef6a55853e8a6665e582277c76a1c6c7d57f63..9364090dde18490a4803662c7eb770d9f6354b1c 100644
--- a/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
+++ b/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
@@ -12,20 +12,13 @@ class SymmetryBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
   std::ostream&
   _write(std::ostream& os) const final
   {
-    os << "symmetry(" << m_name << ',' << *m_boundary_descriptor << ")";
+    os << "symmetry(" << *m_boundary_descriptor << ")";
     return os;
   }
 
-  const std::string_view m_name;
-
   std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
 
  public:
-  std::string_view
-  name() const
-  {
-    return m_name;
-  }
   const IBoundaryDescriptor&
   boundaryDescriptor() const final
   {
diff --git a/src/scheme/VectorDiamondScheme.cpp b/src/scheme/VectorDiamondScheme.cpp
deleted file mode 100644
index 751488fc3acb2f35bfa12a404972a1b9100b0485..0000000000000000000000000000000000000000
--- a/src/scheme/VectorDiamondScheme.cpp
+++ /dev/null
@@ -1,2049 +0,0 @@
-#include <scheme/VectorDiamondScheme.hpp>
-
-#include <mesh/ItemValueVariant.hpp>
-#include <scheme/DiscreteFunctionP0.hpp>
-#include <scheme/DiscreteFunctionUtils.hpp>
-#include <utils/Messenger.hpp>
-
-template <size_t Dimension>
-class VectorDiamondSchemeHandler::InterpolationWeightsManager
-{
- private:
-  std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh;
-  FaceValue<const bool> m_primal_face_is_on_boundary;
-  NodeValue<const bool> m_primal_node_is_on_boundary;
-  FaceValue<const bool> m_primal_face_is_symmetry;
-  CellValuePerNode<double> m_w_rj;
-  FaceValuePerNode<double> m_w_rl;
-
- public:
-  InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh,
-                              FaceValue<const bool> primal_face_is_on_boundary,
-                              NodeValue<const bool> primal_node_is_on_boundary,
-                              FaceValue<const bool> primal_face_is_symmetry)
-    : m_mesh(mesh),
-      m_primal_face_is_on_boundary(primal_face_is_on_boundary),
-      m_primal_node_is_on_boundary(primal_node_is_on_boundary),
-      m_primal_face_is_symmetry(primal_face_is_symmetry)
-  {}
-  ~InterpolationWeightsManager() = default;
-  CellValuePerNode<double>&
-  wrj()
-  {
-    return m_w_rj;
-  }
-  FaceValuePerNode<double>&
-  wrl()
-  {
-    return m_w_rl;
-  }
-  void
-  compute()
-  {
-    using MeshDataType      = MeshData<Dimension>;
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
-
-    const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
-    const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-    const auto& node_to_cell_matrix                  = m_mesh->connectivity().nodeToCellMatrix();
-    const auto& node_to_face_matrix                  = m_mesh->connectivity().nodeToFaceMatrix();
-    const auto& face_to_cell_matrix                  = m_mesh->connectivity().faceToCellMatrix();
-
-    CellValuePerNode<double> w_rj{m_mesh->connectivity()};
-    FaceValuePerNode<double> w_rl{m_mesh->connectivity()};
-
-    const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
-    auto project_to_face = [&](const TinyVector<Dimension>& x, const FaceId face_id) -> const TinyVector<Dimension> {
-      TinyVector<Dimension> proj;
-      const TinyVector<Dimension> nil = primal_nlr(face_id, 0);
-      proj                            = x - dot((x - xl[face_id]), nil) * nil;
-      return proj;
-    };
-
-    for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
-      w_rl[i] = std::numeric_limits<double>::signaling_NaN();
-    }
-
-    for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
-      SmallVector<double> b{Dimension + 1};
-      b[0] = 1;
-      for (size_t i = 1; i < Dimension + 1; i++) {
-        b[i] = xr[i_node][i - 1];
-      }
-      const auto& node_to_cell = node_to_cell_matrix[i_node];
-
-      if (not m_primal_node_is_on_boundary[i_node]) {
-        SmallMatrix<double> A{Dimension + 1, node_to_cell.size()};
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          A(0, j) = 1;
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J = node_to_cell[j];
-            A(i, j)        = xj[J][i - 1];
-          }
-        }
-
-        SmallVector<double> x{node_to_cell.size()};
-        x = zero;
-
-        LeastSquareSolver ls_solver;
-        ls_solver.solveLocalSystem(A, x, b);
-
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          w_rj(i_node, j) = x[j];
-        }
-      } else {
-        int nb_face_used = 0;
-        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[i_node][i_face];
-          if (m_primal_face_is_on_boundary[face_id]) {
-            nb_face_used++;
-          }
-        }
-        SmallMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
-        for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
-          A(0, j) = 1;
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J = node_to_cell[j];
-            A(i, j)        = xj[J][i - 1];
-          }
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          int cpt_face = 0;
-          for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-            FaceId face_id = node_to_face_matrix[i_node][i_face];
-            if (m_primal_face_is_on_boundary[face_id]) {
-              if (m_primal_face_is_symmetry[face_id]) {
-                for (size_t j = 0; j < face_to_cell_matrix[face_id].size(); ++j) {
-                  const CellId cell_id                 = face_to_cell_matrix[face_id][j];
-                  TinyVector<Dimension> xproj          = project_to_face(xj[cell_id], face_id);
-                  A(i, node_to_cell.size() + cpt_face) = xproj[i - 1];
-                }
-              } else {
-                A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
-              }
-              cpt_face++;
-            }
-          }
-        }
-
-        SmallVector<double> x{node_to_cell.size() + nb_face_used};
-        x = zero;
-
-        LeastSquareSolver ls_solver;
-        ls_solver.solveLocalSystem(A, x, b);
-
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          w_rj(i_node, j) = x[j];
-        }
-        int cpt_face = node_to_cell.size();
-        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[i_node][i_face];
-          if (m_primal_face_is_on_boundary[face_id]) {
-            w_rl(i_node, i_face) = x[cpt_face++];
-          }
-        }
-      }
-    }
-    m_w_rj = w_rj;
-    m_w_rl = w_rl;
-  }
-};
-
-class VectorDiamondSchemeHandler::IVectorDiamondScheme
-{
- public:
-  virtual std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const ItemValueVariant>>
-  getSolution() const                                                            = 0;
-  virtual std::shared_ptr<const DiscreteFunctionVariant> getDualSolution() const = 0;
-  virtual std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>>
-  apply() const = 0;
-
-  IVectorDiamondScheme()          = default;
-  virtual ~IVectorDiamondScheme() = default;
-};
-
-template <size_t Dimension>
-class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSchemeHandler::IVectorDiamondScheme
-{
- private:
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshData<Dimension>;
-
-  DiscreteFunctionP0<Dimension, TinyVector<Dimension>> m_solution;
-  NodeValue<TinyVector<Dimension>> m_nodal_solution;
-  DiscreteFunctionP0<Dimension, TinyVector<Dimension>> m_dual_solution;
-
-  class DirichletBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const TinyVector<Dimension>>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    DirichletBoundaryCondition(const Array<const FaceId>& face_list,
-                               const Array<const TinyVector<Dimension>>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~DirichletBoundaryCondition() = default;
-  };
-
-  class NormalStrainBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const TinyVector<Dimension>>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    NormalStrainBoundaryCondition(const Array<const FaceId>& face_list,
-                                  const Array<const TinyVector<Dimension>>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~NormalStrainBoundaryCondition() = default;
-  };
-
-  class SymmetryBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-   public:
-    SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
-
-    ~SymmetryBoundaryCondition() = default;
-  };
-
- public:
-  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const ItemValueVariant>>
-  getSolution() const final
-  {
-    return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(m_solution),
-                           std::make_shared<ItemValueVariant>(m_nodal_solution));
-  }
-
-  std::shared_ptr<const DiscreteFunctionVariant>
-  getDualSolution() const final
-  {
-    return std::make_shared<DiscreteFunctionVariant>(m_dual_solution);
-  }
-
-  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>>
-  apply() const final
-  {
-    return {std::make_shared<DiscreteFunctionVariant>(m_solution),
-            std::make_shared<DiscreteFunctionVariant>(m_dual_solution)};
-  }
-
-  VectorDiamondScheme(const std::shared_ptr<const MeshType>& mesh,
-                      const DiscreteFunctionP0<Dimension, const double>& alpha,
-                      const DiscreteFunctionP0<Dimension, const double>& dual_lambdab,
-                      const DiscreteFunctionP0<Dimension, const double>& dual_mub,
-                      const DiscreteFunctionP0<Dimension, const double>& dual_lambda,
-                      const DiscreteFunctionP0<Dimension, const double>& dual_mu,
-                      const DiscreteFunctionP0<Dimension, const TinyVector<Dimension>>& source,
-                      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-    : m_solution(mesh), m_dual_solution(DualMeshManager::instance().getDiamondDualMesh(*mesh))
-  {
-    Assert(mesh == alpha.mesh());
-    Assert(mesh == source.mesh());
-    Assert(dual_lambda.mesh() == dual_mu.mesh());
-    Assert(dual_lambdab.mesh() == dual_mu.mesh());
-    Assert(dual_mub.mesh() == dual_mu.mesh());
-    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mu.mesh(),
-           "diffusion coefficient is not defined on the dual mesh!");
-
-    using MeshDataType = MeshData<Dimension>;
-
-    using BoundaryCondition =
-      std::variant<DirichletBoundaryCondition, NormalStrainBoundaryCondition, SymmetryBoundaryCondition>;
-
-    using BoundaryConditionList = std::vector<BoundaryCondition>;
-
-    BoundaryConditionList boundary_condition_list;
-
-    NodeValue<bool> is_dirichlet{mesh->connectivity()};
-    is_dirichlet.fill(false);
-    NodeValue<TinyVector<Dimension>> dirichlet_value{mesh->connectivity()};
-    {
-      TinyVector<Dimension> nan_tiny_vector;
-      for (size_t i = 0; i < Dimension; ++i) {
-        nan_tiny_vector[i] = std::numeric_limits<double>::signaling_NaN();
-      }
-      dirichlet_value.fill(nan_tiny_vector);
-    }
-
-    for (const auto& bc_descriptor : bc_descriptor_list) {
-      bool is_valid_boundary_condition = true;
-
-      switch (bc_descriptor->type()) {
-      case IBoundaryConditionDescriptor::Type::symmetry: {
-        const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-          dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-
-        if constexpr (Dimension > 1) {
-          MeshFlatFaceBoundary<Dimension> mesh_face_boundary =
-            getMeshFlatFaceBoundary(*mesh, sym_bc_descriptor.boundaryDescriptor());
-          boundary_condition_list.push_back(SymmetryBoundaryCondition{mesh_face_boundary.faceList()});
-        } else {
-          throw NotImplementedError("Symmetry conditions are not supported in 1d");
-        }
-
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::dirichlet: {
-        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
-          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
-        if (dirichlet_bc_descriptor.name() == "dirichlet") {
-          if constexpr (Dimension > 1) {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-            const FunctionSymbolId g_id                   = dirichlet_bc_descriptor.rhsSymbolId();
-            Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-              TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(),
-                                                                            mesh_face_boundary.faceList());
-            boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-          } else {
-            throw NotImplementedError("Neumann conditions are not supported in 1d");
-          }
-        } else if (dirichlet_bc_descriptor.name() == "normal_strain") {
-          if constexpr (Dimension > 1) {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-            const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-
-            Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-              TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(),
-                                                                            mesh_face_boundary.faceList());
-            boundary_condition_list.push_back(NormalStrainBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-
-          } else {
-            throw NotImplementedError("Normal strain conditions are not supported in 1d");
-          }
-        }
-#ifdef PUGS_HAS_COSTO
-        else if (dirichlet_bc_descriptor.name() == "cpl_normal_strain") {
-          if constexpr (Dimension > 1) {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-            const int src = 1;
-            const int tag = 200;
-
-            std::vector<double> recv;
-            Array<TinyVector<Dimension>> value_list{mesh_face_boundary.faceList().size()};
-
-            // RECV normal strain at ZOI face center
-            parallel::Messenger::getInstance().myCoupling->recvData(recv, src, tag);
-
-            for (size_t i_face = 0; i_face < mesh_face_boundary.faceList().size(); ++i_face) {
-              for (size_t i_dim = 0; i_dim < Dimension; ++i_dim) {
-                value_list[i_face][i_dim] = recv[i_face * Dimension + i_dim];
-              }
-            }
-            /* std::cout << "\033[01;31m" */
-            /*           << "Face list" */
-            /*           << "\033[00;00m" << std::endl; */
-            /* for (size_t i = 0; i < mesh_face_boundary.faceList().size(); ++i) { */
-            /*   for (size_t idim = 0; idim < Dimension; ++idim) { */
-            /*     std::cout << value_list[i][idim] << " "; */
-            /*   } */
-            /*   std::cout << "\n"; */
-            /* } */
-
-            boundary_condition_list.push_back(NormalStrainBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-          } else {
-            throw NotImplementedError("Normal strain conditions are not supported in 1d");
-          }
-        }
-#endif   // PUGS_HAS_COSTO
-        else {
-          is_valid_boundary_condition = false;
-        }
-        break;
-      }
-      default: {
-        is_valid_boundary_condition = false;
-      }
-      }
-      if (not is_valid_boundary_condition) {
-        std::ostringstream error_msg;
-        error_msg << *bc_descriptor << " is an invalid boundary condition for elasticity equation";
-        throw NormalError(error_msg.str());
-      }
-    }
-
-    if constexpr (Dimension > 1) {
-      const CellValue<const size_t> cell_dof_number = [&] {
-        CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
-        return compute_cell_dof_number;
-      }();
-      size_t number_of_dof = mesh->numberOfCells();
-
-      const FaceValue<const size_t> face_dof_number = [&] {
-        FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
-        compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>) or
-                            (std::is_same_v<T, DirichletBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id = face_list[i_face];
-                  if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
-                    std::ostringstream os;
-                    os << "The face " << face_id << " is used at least twice for boundary conditions";
-                    throw NormalError(os.str());
-                  } else {
-                    compute_face_dof_number[face_id] = number_of_dof++;
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return compute_face_dof_number;
-      }();
-
-      const auto& primal_face_to_node_matrix             = mesh->connectivity().faceToNodeMatrix();
-      const auto& face_to_cell_matrix                    = mesh->connectivity().faceToCellMatrix();
-      const FaceValue<const bool> primal_face_is_neumann = [&] {
-        FaceValue<bool> face_is_neumann{mesh->connectivity()};
-        face_is_neumann.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id     = face_list[i_face];
-                  face_is_neumann[face_id] = true;
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return face_is_neumann;
-      }();
-
-      const FaceValue<const bool> primal_face_is_symmetry = [&] {
-        FaceValue<bool> face_is_symmetry{mesh->connectivity()};
-        face_is_symmetry.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, SymmetryBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id      = face_list[i_face];
-                  face_is_symmetry[face_id] = true;
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return face_is_symmetry;
-      }();
-
-      NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
-      }
-
-      primal_node_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-            NodeId node_id                      = primal_face_to_node_matrix[face_id][i_node];
-            primal_node_is_on_boundary[node_id] = true;
-          }
-        }
-      }
-
-      FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
-      }
-
-      primal_face_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          primal_face_is_on_boundary[face_id] = true;
-        }
-      }
-
-      FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_neumann is incorrect");
-      }
-
-      primal_face_is_dirichlet.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] &&
-                                             (!primal_face_is_neumann[face_id]) && (!primal_face_is_symmetry[face_id]));
-      }
-
-      InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary,
-                                      primal_face_is_symmetry);
-      iwm.compute();
-      CellValuePerNode<double> w_rj = iwm.wrj();
-      FaceValuePerNode<double> w_rl = iwm.wrl();
-
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-      const FaceValue<const TinyVector<Dimension>>& xl               = mesh_data.xl();
-      const CellValue<const TinyVector<Dimension>>& xj               = mesh_data.xj();
-      const auto& node_to_cell_matrix                                = mesh->connectivity().nodeToCellMatrix();
-      const auto& node_to_face_matrix                                = mesh->connectivity().nodeToFaceMatrix();
-      const auto& cell_to_node_matrix                                = mesh->connectivity().cellToNodeMatrix();
-      const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
-
-      {
-        std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
-
-        MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
-
-        std::shared_ptr mapper =
-          DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
-
-        CellValue<const double> dual_muj      = dual_mu.cellValues();
-        CellValue<const double> dual_lambdaj  = dual_lambda.cellValues();
-        CellValue<const double> dual_mubj     = dual_mub.cellValues();
-        CellValue<const double> dual_lambdabj = dual_lambdab.cellValues();
-
-        CellValue<const TinyVector<Dimension>> fj = source.cellValues();
-
-        const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
-
-        const FaceValue<const double> mes_l = [&] {
-          if constexpr (Dimension == 1) {
-            FaceValue<double> compute_mes_l{mesh->connectivity()};
-            compute_mes_l.fill(1);
-            return compute_mes_l;
-          } else {
-            return mesh_data.ll();
-          }
-        }();
-
-        const CellValue<const double> dual_mes_l_j = [=] {
-          CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
-          mapper->toDualCell(mes_l, compute_mes_j);
-
-          return compute_mes_j;
-        }();
-
-        const CellValue<const double> primal_Vj = mesh_data.Vj();
-
-        CellValue<const FaceId> dual_cell_face_id = [=]() {
-          CellValue<FaceId> computed_dual_cell_face_id{diamond_mesh->connectivity()};
-          FaceValue<FaceId> primal_face_id{mesh->connectivity()};
-          parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { primal_face_id[face_id] = face_id; });
-
-          mapper->toDualCell(primal_face_id, computed_dual_cell_face_id);
-
-          return computed_dual_cell_face_id;
-        }();
-
-        FaceValue<const CellId> face_dual_cell_id = [=]() {
-          FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
-          CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
-
-          mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
-
-          return computed_face_dual_cell_id;
-        }();
-
-        NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
-          CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
-          cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> node_primal_id{mesh->connectivity()};
-
-          parallel_for(
-            mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
-
-          NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
-
-          mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
-
-          return computed_dual_node_primal_node_id;
-        }();
-
-        CellValue<NodeId> primal_cell_dual_node_id = [=]() {
-          CellValue<NodeId> cell_id{mesh->connectivity()};
-          NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
-          node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
-
-          CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
-
-          mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
-
-          return cell_id;
-        }();
-        const auto& dual_Cjr                     = diamond_mesh_data.Cjr();
-        FaceValue<TinyVector<Dimension>> dualClj = [&] {
-          FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
-          const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
-          const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
-          parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-              const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-              for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
-                CellId cell_id            = primal_face_to_cell[i];
-                const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
-                for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
-                     i_dual_cell++) {
-                  const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
-                  if (face_dual_cell_id[face_id] == dual_cell_id) {
-                    for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
-                         i_dual_node++) {
-                      const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                      if (final_dual_node_id == dual_node_id) {
-                        computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
-                      }
-                    }
-                  }
-                }
-              }
-            });
-          return computedClj;
-        }();
-
-        FaceValue<TinyVector<Dimension>> nlj = [&] {
-          FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
-          parallel_for(
-            mesh->numberOfFaces(),
-            PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
-          return computedNlj;
-        }();
-
-        FaceValue<const double> alpha_lambda_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_lambdaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        FaceValue<const double> alpha_mu_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_muj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        FaceValue<const double> alpha_lambdab_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_lambdabj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        FaceValue<const double> alpha_mub_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_mubj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        const TinyMatrix<Dimension> I = identity;
-
-        const Array<int> non_zeros{number_of_dof * Dimension};
-        non_zeros.fill(Dimension * Dimension);
-        CRSMatrixDescriptor<double> S(number_of_dof * Dimension, number_of_dof * Dimension, non_zeros);
-
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const double beta_mu_l          = l2Norm(dualClj[face_id]) * alpha_mu_l[face_id] * mes_l[face_id];
-          const double beta_lambda_l      = l2Norm(dualClj[face_id]) * alpha_lambda_l[face_id] * mes_l[face_id];
-          const double beta_mub_l         = l2Norm(dualClj[face_id]) * alpha_mub_l[face_id] * mes_l[face_id];
-          const double beta_lambdab_l     = l2Norm(dualClj[face_id]) * alpha_lambdab_l[face_id] * mes_l[face_id];
-          const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-          for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
-            const CellId i_id                      = primal_face_to_cell[i_cell];
-            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
-
-            const TinyVector<Dimension> nil = [&] {
-              if (is_face_reversed_for_cell_i) {
-                return -nlj[face_id];
-              } else {
-                return nlj[face_id];
-              }
-            }();
-            TinyMatrix<Dimension> M =
-              beta_mu_l * I + beta_mu_l * tensorProduct(nil, nil) + beta_lambda_l * tensorProduct(nil, nil);
-            TinyMatrix<Dimension> Mb =
-              beta_mub_l * I + beta_mub_l * tensorProduct(nil, nil) + beta_lambdab_l * tensorProduct(nil, nil);
-            TinyMatrix<Dimension> N = 1.e0 * tensorProduct(nil, nil);
-            double coef_adim        = beta_mu_l + beta_lambdab_l;
-            for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
-              const CellId j_id = primal_face_to_cell[j_cell];
-              if (i_cell == j_cell) {
-                for (size_t i = 0; i < Dimension; ++i) {
-                  for (size_t j = 0; j < Dimension; ++j) {
-                    S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) += M(i, j);
-                    if (primal_face_is_neumann[face_id]) {
-                      S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) -=
-                        1.e0 * Mb(i, j);
-                      // S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + j) +=
-                      //   1.e-10 * Mb(i, j);
-                    }
-                    if (primal_face_is_symmetry[face_id]) {
-                      S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
-                        ((i == j) ? -coef_adim : 0) + coef_adim * N(i, j);
-                      S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + j) +=
-                        (i == j) ? coef_adim : 0;
-                    }
-                  }
-                }
-              } else {
-                for (size_t i = 0; i < Dimension; ++i) {
-                  for (size_t j = 0; j < Dimension; ++j) {
-                    S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -= M(i, j);
-                  }
-                }
-              }
-            }
-          }
-        }
-
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          for (size_t i = 0; i < Dimension; ++i) {
-            const size_t j = cell_dof_number[cell_id] * Dimension + i;
-            S(j, j) += alpha[cell_id] * primal_Vj[cell_id];
-          }
-        }
-
-        const auto& dual_cell_to_node_matrix   = diamond_mesh->connectivity().cellToNodeMatrix();
-        const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const double alpha_mu_face_id      = mes_l[face_id] * alpha_mu_l[face_id];
-          const double alpha_lambda_face_id  = mes_l[face_id] * alpha_lambda_l[face_id];
-          const double alpha_mub_face_id     = mes_l[face_id] * alpha_mub_l[face_id];
-          const double alpha_lambdab_face_id = mes_l[face_id] * alpha_lambdab_l[face_id];
-
-          for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
-            CellId i_id                            = face_to_cell_matrix[face_id][i_face_cell];
-            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
-
-            for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-              NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
-
-              const TinyVector<Dimension> nil = [&] {
-                if (is_face_reversed_for_cell_i) {
-                  return -nlj[face_id];
-                } else {
-                  return nlj[face_id];
-                }
-              }();
-
-              CellId dual_cell_id = face_dual_cell_id[face_id];
-
-              for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
-                const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                if (dual_node_primal_node_id[dual_node_id] == node_id) {
-                  const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-
-                  TinyMatrix<Dimension> M = alpha_mu_face_id * dot(Clr, nil) * I +
-                                            alpha_mu_face_id * tensorProduct(Clr, nil) +
-                                            alpha_lambda_face_id * tensorProduct(nil, Clr);
-                  TinyMatrix<Dimension> Mb = alpha_mub_face_id * dot(Clr, nil) * I +
-                                             alpha_mub_face_id * tensorProduct(Clr, nil) +
-                                             alpha_lambdab_face_id * tensorProduct(nil, Clr);
-
-                  for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
-                    CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
-                    for (size_t i = 0; i < Dimension; ++i) {
-                      for (size_t j = 0; j < Dimension; ++j) {
-                        S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -=
-                          w_rj(node_id, j_cell) * M(i, j);
-                        if (primal_face_is_neumann[face_id]) {
-                          S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
-                            1.e0 * w_rj(node_id, j_cell) * Mb(i, j);
-                        }
-                      }
-                    }
-                  }
-                  if (primal_node_is_on_boundary[node_id]) {
-                    for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
-                      FaceId l_id = node_to_face_matrix[node_id][l_face];
-                      if (primal_face_is_on_boundary[l_id]) {
-                        for (size_t i = 0; i < Dimension; ++i) {
-                          for (size_t j = 0; j < Dimension; ++j) {
-                            // Mb?
-                            S(cell_dof_number[i_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) -=
-                              w_rl(node_id, l_face) * M(i, j);
-                          }
-                        }
-                        if (primal_face_is_neumann[face_id]) {
-                          for (size_t i = 0; i < Dimension; ++i) {
-                            for (size_t j = 0; j < Dimension; ++j) {
-                              S(face_dof_number[face_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) +=
-                                1.e0 * w_rl(node_id, l_face) * Mb(i, j);
-                            }
-                          }
-                        }
-                      }
-                    }
-                  }
-                }
-              }
-              //            }
-            }
-          }
-        }
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          if (primal_face_is_dirichlet[face_id]) {
-            for (size_t i = 0; i < Dimension; ++i) {
-              S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + i) += 1.e0;
-            }
-          }
-        }
-
-        Vector<double> b{number_of_dof * Dimension};
-        b = zero;
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          for (size_t i = 0; i < Dimension; ++i) {
-            b[(cell_dof_number[cell_id] * Dimension) + i] = primal_Vj[cell_id] * fj[cell_id][i];
-          }
-        }
-
-        // Dirichlet
-        NodeValue<bool> node_tag{mesh->connectivity()};
-        node_tag.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
-                const auto& face_list  = bc.faceList();
-                const auto& value_list = bc.valueList();
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id = face_list[i_face];
-
-                  for (size_t i = 0; i < Dimension; ++i) {
-                    b[(face_dof_number[face_id] * Dimension) + i] += 1.e0 * value_list[i_face][i];
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
-                const auto& face_list  = bc.faceList();
-                const auto& value_list = bc.valueList();
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  FaceId face_id = face_list[i_face];
-                  for (size_t i = 0; i < Dimension; ++i) {
-                    b[face_dof_number[face_id] * Dimension + i] +=
-                      1.e0 * mes_l[face_id] * value_list[i_face][i];   // sign
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        CRSMatrix A{S.getCRSMatrix()};
-        Vector<double> U{number_of_dof * Dimension};
-        U        = zero;
-        Vector r = A * U - b;
-        std::cout << "initial (real) residu = " << std::sqrt(dot(r, r)) << '\n';
-
-        LinearSolver solver;
-        solver.solveLocalSystem(A, U, b);
-
-        r = A * U - b;
-
-        std::cout << "final (real) residu = " << std::sqrt(dot(r, r)) << '\n';
-
-        m_nodal_solution = NodeValue<TinyVector<Dimension>>(mesh->connectivity());
-        m_nodal_solution.fill(zero);
-
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-            for (size_t i = 0; i < Dimension; ++i) {
-              m_solution[cell_id][i] = U[(cell_dof_number[cell_id] * Dimension) + i];
-            }
-          });
-
-        for (NodeId i_node = 0; i_node < mesh->numberOfNodes(); ++i_node) {
-          const auto& node_to_cell = node_to_cell_matrix[i_node];
-
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J           = node_to_cell[j];
-            const auto& cell_to_node = cell_to_node_matrix[J];
-            for (size_t k = 0; k < cell_to_node.size(); k++) {
-              const NodeId K = cell_to_node[k];
-              if (K == i_node) {
-                m_nodal_solution[i_node] += w_rj(i_node, j) * m_solution[J];
-              }
-            }
-          }
-          /* if (not primal_node_is_on_boundary[i_node]) { */
-          const auto& node_to_face = node_to_face_matrix[i_node];
-          for (size_t j = 0; j < node_to_face.size(); j++) {
-            const FaceId J           = node_to_face[j];
-            const auto& face_to_node = primal_face_to_node_matrix[J];
-            if (primal_face_is_on_boundary[J]) {
-              for (size_t k = 0; k < face_to_node.size(); k++) {
-                const NodeId K = face_to_node[k];
-                if (K == i_node) {
-                  for (size_t i = 0; i < Dimension; ++i) {
-                    m_nodal_solution[i_node][i] += w_rl(i_node, j) * U[(face_dof_number[J] * Dimension) + i];
-                  }
-                  /* m_nodal_solution[i_node] += w_rl(i_node, j) * solution[J]; */
-                }
-              }
-            }
-            /* m_nodal_solution[i_node][j] = w_rj(i_node, j) * solution[J][k]; */
-          }
-          /* } */
-        }
-
-        m_dual_solution.fill(zero);
-        const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
-        for (CellId cell_id = 0; cell_id < diamond_mesh->numberOfCells(); ++cell_id) {
-          const FaceId face_id = dual_cell_face_id[cell_id];
-          if (primal_face_is_on_boundary[face_id]) {
-            for (size_t i = 0; i < Dimension; ++i) {
-              m_dual_solution[cell_id][i] = U[(face_dof_number[face_id] * Dimension) + i];
-            }
-          } else {
-            CellId cell_id1 = face_to_cell_matrix[face_id][0];
-            CellId cell_id2 = face_to_cell_matrix[face_id][1];
-            for (size_t i = 0; i < Dimension; ++i) {
-              m_dual_solution[cell_id][i] =
-                0.5 * (U[(cell_dof_number[cell_id1] * Dimension) + i] + U[(cell_dof_number[cell_id2] * Dimension) + i]);
-            }
-          }
-        }
-      }
-      // provide a source for E?
-      // computeEnergyUpdate(mesh, dual_lambdab, dual_mub, m_solution,
-      //                     m_dual_solution,   // f,
-      //                     bc_descriptor_list);
-      // computeEnergyUpdate(mesh, alpha, dual_lambdab, dual_mub, dual_lambda, dual_mu, m_solution,
-      //                     m_dual_solution,   // f,
-      //                     bc_descriptor_list);
-    } else {
-      throw NotImplementedError("not done in 1d");
-    }
-  }
-};
-// NEW CLASS
-template <size_t Dimension>
-class EnergyComputerHandler::InterpolationWeightsManager
-{
- private:
-  std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh;
-  FaceValue<const bool> m_primal_face_is_on_boundary;
-  NodeValue<const bool> m_primal_node_is_on_boundary;
-  FaceValue<const bool> m_primal_face_is_symmetry;
-  CellValuePerNode<double> m_w_rj;
-  FaceValuePerNode<double> m_w_rl;
-
- public:
-  InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh,
-                              FaceValue<const bool> primal_face_is_on_boundary,
-                              NodeValue<const bool> primal_node_is_on_boundary,
-                              FaceValue<const bool> primal_face_is_symmetry)
-    : m_mesh(mesh),
-      m_primal_face_is_on_boundary(primal_face_is_on_boundary),
-      m_primal_node_is_on_boundary(primal_node_is_on_boundary),
-      m_primal_face_is_symmetry(primal_face_is_symmetry)
-  {}
-  ~InterpolationWeightsManager() = default;
-  CellValuePerNode<double>&
-  wrj()
-  {
-    return m_w_rj;
-  }
-  FaceValuePerNode<double>&
-  wrl()
-  {
-    return m_w_rl;
-  }
-  void
-  compute()
-  {
-    using MeshDataType      = MeshData<Dimension>;
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
-
-    const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
-    const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-    const auto& node_to_cell_matrix                  = m_mesh->connectivity().nodeToCellMatrix();
-    const auto& node_to_face_matrix                  = m_mesh->connectivity().nodeToFaceMatrix();
-    const auto& face_to_cell_matrix                  = m_mesh->connectivity().faceToCellMatrix();
-
-    CellValuePerNode<double> w_rj{m_mesh->connectivity()};
-    FaceValuePerNode<double> w_rl{m_mesh->connectivity()};
-
-    const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
-    auto project_to_face = [&](const TinyVector<Dimension>& x, const FaceId face_id) -> const TinyVector<Dimension> {
-      TinyVector<Dimension> proj;
-      const TinyVector<Dimension> nil = primal_nlr(face_id, 0);
-      proj                            = x - dot((x - xl[face_id]), nil) * nil;
-      return proj;
-    };
-
-    for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
-      w_rl[i] = std::numeric_limits<double>::signaling_NaN();
-    }
-
-    for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
-      SmallVector<double> b{Dimension + 1};
-      b[0] = 1;
-      for (size_t i = 1; i < Dimension + 1; i++) {
-        b[i] = xr[i_node][i - 1];
-      }
-      const auto& node_to_cell = node_to_cell_matrix[i_node];
-
-      if (not m_primal_node_is_on_boundary[i_node]) {
-        SmallMatrix<double> A{Dimension + 1, node_to_cell.size()};
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          A(0, j) = 1;
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J = node_to_cell[j];
-            A(i, j)        = xj[J][i - 1];
-          }
-        }
-
-        SmallVector<double> x{node_to_cell.size()};
-        x = zero;
-
-        LeastSquareSolver ls_solver;
-        ls_solver.solveLocalSystem(A, x, b);
-
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          w_rj(i_node, j) = x[j];
-        }
-      } else {
-        int nb_face_used = 0;
-        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[i_node][i_face];
-          if (m_primal_face_is_on_boundary[face_id]) {
-            nb_face_used++;
-          }
-        }
-        SmallMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
-        for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
-          A(0, j) = 1;
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J = node_to_cell[j];
-            A(i, j)        = xj[J][i - 1];
-          }
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          int cpt_face = 0;
-          for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-            FaceId face_id = node_to_face_matrix[i_node][i_face];
-            if (m_primal_face_is_on_boundary[face_id]) {
-              if (m_primal_face_is_symmetry[face_id]) {
-                for (size_t j = 0; j < face_to_cell_matrix[face_id].size(); ++j) {
-                  const CellId cell_id                 = face_to_cell_matrix[face_id][j];
-                  TinyVector<Dimension> xproj          = project_to_face(xj[cell_id], face_id);
-                  A(i, node_to_cell.size() + cpt_face) = xproj[i - 1];
-                }
-              } else {
-                A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
-              }
-              cpt_face++;
-            }
-          }
-        }
-
-        SmallVector<double> x{node_to_cell.size() + nb_face_used};
-        x = zero;
-
-        LeastSquareSolver ls_solver;
-        ls_solver.solveLocalSystem(A, x, b);
-
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          w_rj(i_node, j) = x[j];
-        }
-        int cpt_face = node_to_cell.size();
-        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[i_node][i_face];
-          if (m_primal_face_is_on_boundary[face_id]) {
-            w_rl(i_node, i_face) = x[cpt_face++];
-          }
-        }
-      }
-    }
-    m_w_rj = w_rj;
-    m_w_rl = w_rl;
-  }
-};
-class EnergyComputerHandler::IEnergyComputer
-{
- public:
-  virtual std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>>
-  apply() const = 0;
-
-  IEnergyComputer()          = default;
-  virtual ~IEnergyComputer() = default;
-};
-
-template <size_t Dimension>
-class EnergyComputerHandler::EnergyComputer : public EnergyComputerHandler::IEnergyComputer
-{
- private:
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshData<Dimension>;
-
-  DiscreteFunctionP0<Dimension, TinyVector<Dimension>> m_solution;
-  DiscreteFunctionP0<Dimension, TinyVector<Dimension>> m_dual_solution;
-  DiscreteFunctionP0<Dimension, double> m_energy_delta;
-
-  class DirichletBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const TinyVector<Dimension>>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    DirichletBoundaryCondition(const Array<const FaceId>& face_list,
-                               const Array<const TinyVector<Dimension>>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~DirichletBoundaryCondition() = default;
-  };
-
-  class NormalStrainBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const TinyVector<Dimension>>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    NormalStrainBoundaryCondition(const Array<const FaceId>& face_list,
-                                  const Array<const TinyVector<Dimension>>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~NormalStrainBoundaryCondition() = default;
-  };
-
-  class SymmetryBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-   public:
-    SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
-
-    ~SymmetryBoundaryCondition() = default;
-  };
-
- public:
-  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>>
-  apply() const final
-  {
-    return {std::make_shared<DiscreteFunctionVariant>(m_energy_delta),
-            std::make_shared<DiscreteFunctionVariant>(m_dual_solution)};
-  }
-
-  // compute the fluxes
-  EnergyComputer(const std::shared_ptr<const MeshType>& mesh,
-                 const DiscreteFunctionP0<Dimension, const double>& dual_lambdab,
-                 const DiscreteFunctionP0<Dimension, const double>& dual_mub,
-                 const DiscreteFunctionP0<Dimension, const TinyVector<Dimension>>& U,
-                 const DiscreteFunctionP0<Dimension, const TinyVector<Dimension>>& dual_U,
-                 const DiscreteFunctionP0<Dimension, const TinyVector<Dimension>>& source,
-                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-    : m_solution(mesh), m_dual_solution(DualMeshManager::instance().getDiamondDualMesh(*mesh)), m_energy_delta(mesh)
-  {
-    Assert(mesh == U.mesh());
-    Assert(dual_lambdab.mesh() == dual_U.mesh());
-    Assert(U.mesh() == source.mesh());
-    Assert(dual_lambdab.mesh() == dual_mub.mesh());
-    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mub.mesh(),
-           "diffusion coefficient is not defined on the dual mesh!");
-
-    using MeshDataType = MeshData<Dimension>;
-
-    using BoundaryCondition =
-      std::variant<DirichletBoundaryCondition, NormalStrainBoundaryCondition, SymmetryBoundaryCondition>;
-
-    using BoundaryConditionList = std::vector<BoundaryCondition>;
-
-    BoundaryConditionList boundary_condition_list;
-
-    NodeValue<bool> is_dirichlet{mesh->connectivity()};
-    is_dirichlet.fill(false);
-    NodeValue<TinyVector<Dimension>> dirichlet_value{mesh->connectivity()};
-    {
-      TinyVector<Dimension> nan_tiny_vector;
-      for (size_t i = 0; i < Dimension; ++i) {
-        nan_tiny_vector[i] = std::numeric_limits<double>::signaling_NaN();
-      }
-      dirichlet_value.fill(nan_tiny_vector);
-    }
-
-    for (const auto& bc_descriptor : bc_descriptor_list) {
-      bool is_valid_boundary_condition = true;
-
-      switch (bc_descriptor->type()) {
-      case IBoundaryConditionDescriptor::Type::symmetry: {
-        const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-          dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-
-        if constexpr (Dimension > 1) {
-          MeshFlatFaceBoundary<Dimension> mesh_face_boundary =
-            getMeshFlatFaceBoundary(*mesh, sym_bc_descriptor.boundaryDescriptor());
-          boundary_condition_list.push_back(SymmetryBoundaryCondition{mesh_face_boundary.faceList()});
-        } else {
-          throw NotImplementedError("Symmetry conditions are not supported in 1d");
-        }
-
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::dirichlet: {
-        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
-          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
-        if (dirichlet_bc_descriptor.name() == "dirichlet") {
-          if constexpr (Dimension > 1) {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-            const FunctionSymbolId g_id                   = dirichlet_bc_descriptor.rhsSymbolId();
-            Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-              TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(),
-                                                                            mesh_face_boundary.faceList());
-            boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-          } else {
-            throw NotImplementedError("Neumann conditions are not supported in 1d");
-          }
-        } else if (dirichlet_bc_descriptor.name() == "normal_strain") {
-          if constexpr (Dimension > 1) {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-            const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-
-            Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-              TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(),
-                                                                            mesh_face_boundary.faceList());
-            boundary_condition_list.push_back(NormalStrainBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-
-          } else {
-            throw NotImplementedError("Normal strain conditions are not supported in 1d");
-          }
-        } else {
-          is_valid_boundary_condition = false;
-        }
-        break;
-      }
-      default: {
-        is_valid_boundary_condition = false;
-      }
-      }
-      if (not is_valid_boundary_condition) {
-        std::ostringstream error_msg;
-        error_msg << *bc_descriptor << " is an invalid boundary condition for elasticity equation";
-        throw NormalError(error_msg.str());
-      }
-    }
-
-    if constexpr (Dimension > 1) {
-      const CellValue<const size_t> cell_dof_number = [&] {
-        CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
-        return compute_cell_dof_number;
-      }();
-      size_t number_of_dof = mesh->numberOfCells();
-
-      const FaceValue<const size_t> face_dof_number = [&] {
-        FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
-        compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>) or
-                            (std::is_same_v<T, DirichletBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id = face_list[i_face];
-                  if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
-                    std::ostringstream os;
-                    os << "The face " << face_id << " is used at least twice for boundary conditions";
-                    throw NormalError(os.str());
-                  } else {
-                    compute_face_dof_number[face_id] = number_of_dof++;
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return compute_face_dof_number;
-      }();
-
-      const auto& primal_face_to_node_matrix             = mesh->connectivity().faceToNodeMatrix();
-      const auto& face_to_cell_matrix                    = mesh->connectivity().faceToCellMatrix();
-      const FaceValue<const bool> primal_face_is_neumann = [&] {
-        FaceValue<bool> face_is_neumann{mesh->connectivity()};
-        face_is_neumann.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id     = face_list[i_face];
-                  face_is_neumann[face_id] = true;
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return face_is_neumann;
-      }();
-
-      const FaceValue<const bool> primal_face_is_symmetry = [&] {
-        FaceValue<bool> face_is_symmetry{mesh->connectivity()};
-        face_is_symmetry.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, SymmetryBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id      = face_list[i_face];
-                  face_is_symmetry[face_id] = true;
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return face_is_symmetry;
-      }();
-
-      NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
-      }
-
-      primal_node_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-            NodeId node_id                      = primal_face_to_node_matrix[face_id][i_node];
-            primal_node_is_on_boundary[node_id] = true;
-          }
-        }
-      }
-
-      FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
-      }
-
-      primal_face_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          primal_face_is_on_boundary[face_id] = true;
-        }
-      }
-
-      FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_neumann is incorrect");
-      }
-
-      primal_face_is_dirichlet.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] &&
-                                             (!primal_face_is_neumann[face_id]) && (!primal_face_is_symmetry[face_id]));
-      }
-
-      InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary,
-                                      primal_face_is_symmetry);
-      iwm.compute();
-      CellValuePerNode<double> w_rj = iwm.wrj();
-      FaceValuePerNode<double> w_rl = iwm.wrl();
-
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-      const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
-      const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-      // const auto& node_to_cell_matrix                                = mesh->connectivity().nodeToCellMatrix();
-      const auto& node_to_face_matrix                                = mesh->connectivity().nodeToFaceMatrix();
-      const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
-
-      {
-        std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
-
-        MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
-
-        std::shared_ptr mapper =
-          DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
-
-        CellValue<const double> dual_mubj     = dual_mub.cellValues();
-        CellValue<const double> dual_lambdabj = dual_lambdab.cellValues();
-
-        CellValue<const TinyVector<Dimension>> velocity = U.cellValues();
-
-        const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
-
-        const FaceValue<const double> mes_l = [&] {
-          if constexpr (Dimension == 1) {
-            FaceValue<double> compute_mes_l{mesh->connectivity()};
-            compute_mes_l.fill(1);
-            return compute_mes_l;
-          } else {
-            return mesh_data.ll();
-          }
-        }();
-
-        const CellValue<const double> dual_mes_l_j = [=] {
-          CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
-          mapper->toDualCell(mes_l, compute_mes_j);
-
-          return compute_mes_j;
-        }();
-
-        const CellValue<const double> primal_Vj   = mesh_data.Vj();
-        FaceValue<const CellId> face_dual_cell_id = [=]() {
-          FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
-          CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
-
-          mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
-
-          return computed_face_dual_cell_id;
-        }();
-
-        NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
-          CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
-          cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> node_primal_id{mesh->connectivity()};
-
-          parallel_for(
-            mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
-
-          NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
-
-          mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
-
-          return computed_dual_node_primal_node_id;
-        }();
-
-        CellValue<NodeId> primal_cell_dual_node_id = [=]() {
-          CellValue<NodeId> cell_id{mesh->connectivity()};
-          NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
-          node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
-
-          CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
-
-          mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
-
-          return cell_id;
-        }();
-        const auto& dual_Cjr                     = diamond_mesh_data.Cjr();
-        FaceValue<TinyVector<Dimension>> dualClj = [&] {
-          FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
-          const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
-          const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
-          parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-              const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-              for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
-                CellId cell_id            = primal_face_to_cell[i];
-                const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
-                for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
-                     i_dual_cell++) {
-                  const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
-                  if (face_dual_cell_id[face_id] == dual_cell_id) {
-                    for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
-                         i_dual_node++) {
-                      const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                      if (final_dual_node_id == dual_node_id) {
-                        computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
-                      }
-                    }
-                  }
-                }
-              }
-            });
-          return computedClj;
-        }();
-
-        FaceValue<TinyVector<Dimension>> nlj = [&] {
-          FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
-          parallel_for(
-            mesh->numberOfFaces(),
-            PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
-          return computedNlj;
-        }();
-
-        // FaceValue<const double> alpha_lambda_l = [&] {
-        //   CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-        //   parallel_for(
-        //     diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-        //       alpha_j[diamond_cell_id] = dual_lambdaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-        //     });
-
-        //   FaceValue<double> computed_alpha_l{mesh->connectivity()};
-        //   mapper->fromDualCell(alpha_j, computed_alpha_l);
-        //   return computed_alpha_l;
-        // }();
-
-        // FaceValue<const double> alpha_mu_l = [&] {
-        //   CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-        //   parallel_for(
-        //     diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-        //       alpha_j[diamond_cell_id] = dual_muj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-        //     });
-
-        //   FaceValue<double> computed_alpha_l{mesh->connectivity()};
-        //   mapper->fromDualCell(alpha_j, computed_alpha_l);
-        //   return computed_alpha_l;
-        // }();
-
-        FaceValue<const double> alpha_lambdab_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_lambdabj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        FaceValue<const double> alpha_mub_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_mubj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        const TinyMatrix<Dimension> I             = identity;
-        CellValue<const FaceId> dual_cell_face_id = [=]() {
-          CellValue<FaceId> computed_dual_cell_face_id{diamond_mesh->connectivity()};
-          FaceValue<FaceId> primal_face_id{mesh->connectivity()};
-          parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { primal_face_id[face_id] = face_id; });
-
-          mapper->toDualCell(primal_face_id, computed_dual_cell_face_id);
-
-          return computed_dual_cell_face_id;
-        }();
-
-        auto& dual_solution = dual_U;
-
-        CellValuePerFace<double> flux{mesh->connectivity()};
-        parallel_for(
-          flux.numberOfValues(), PUGS_LAMBDA(size_t jl) { flux[jl] = 0; });
-
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const double beta_mub_l         = l2Norm(dualClj[face_id]) * alpha_mub_l[face_id] * mes_l[face_id];
-          const double beta_lambdab_l     = l2Norm(dualClj[face_id]) * alpha_lambdab_l[face_id] * mes_l[face_id];
-          const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-          for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
-            const CellId i_id                      = primal_face_to_cell[i_cell];
-            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
-
-            const TinyVector<Dimension> nil = [&] {
-              if (is_face_reversed_for_cell_i) {
-                return -nlj[face_id];
-              } else {
-                return nlj[face_id];
-              }
-            }();
-            TinyMatrix<Dimension> M =
-              beta_mub_l * I + beta_mub_l * tensorProduct(nil, nil) + beta_lambdab_l * tensorProduct(nil, nil);
-            //            TinyMatrix<Dimension, double> Mb =
-            //  beta_mub_l * I + beta_mub_l * tensorProduct(nil, nil) + beta_lambdab_l * tensorProduct(nil, nil);
-            // TinyMatrix<Dimension, double> N = 1.e0 * tensorProduct(nil, nil);
-
-            for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
-              const CellId j_id = primal_face_to_cell[j_cell];
-              if (i_cell == j_cell) {
-                flux(face_id, i_cell) += dot(M * velocity[i_id], dual_solution[face_dual_cell_id[face_id]]);
-              } else {
-                flux(face_id, i_cell) -= dot(M * velocity[j_id], dual_solution[face_dual_cell_id[face_id]]);
-              }
-            }
-          }
-        }
-
-        const auto& dual_cell_to_node_matrix   = diamond_mesh->connectivity().cellToNodeMatrix();
-        const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          // const double alpha_mu_face_id     = mes_l[face_id] * alpha_mu_l[face_id];
-          // const double alpha_lambda_face_id = mes_l[face_id] * alpha_lambda_l[face_id];
-          const double alpha_mub_face_id     = mes_l[face_id] * alpha_mub_l[face_id];
-          const double alpha_lambdab_face_id = mes_l[face_id] * alpha_lambdab_l[face_id];
-
-          for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
-            CellId i_id                            = face_to_cell_matrix[face_id][i_face_cell];
-            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
-
-            for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-              NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
-
-              const TinyVector<Dimension> nil = [&] {
-                if (is_face_reversed_for_cell_i) {
-                  return -nlj[face_id];
-                } else {
-                  return nlj[face_id];
-                }
-              }();
-
-              CellId dual_cell_id = face_dual_cell_id[face_id];
-
-              for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
-                const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                if (dual_node_primal_node_id[dual_node_id] == node_id) {
-                  const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-
-                  TinyMatrix<Dimension> M = alpha_mub_face_id * dot(Clr, nil) * I +
-                                            alpha_mub_face_id * tensorProduct(Clr, nil) +
-                                            alpha_lambdab_face_id * tensorProduct(nil, Clr);
-                  // TinyMatrix<Dimension, double> Mb = alpha_mub_face_id * dot(Clr, nil) * I +
-                  //                                    alpha_mub_face_id * tensorProduct(Clr, nil) +
-                  //                                    alpha_lambdab_face_id * tensorProduct(nil, Clr);
-
-                  for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
-                    CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
-                    flux(face_id, i_face_cell) -=
-                      w_rj(node_id, j_cell) * dot(M * velocity[j_id], dual_solution[dual_cell_id]);
-                  }
-                  if (primal_node_is_on_boundary[node_id]) {
-                    for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
-                      FaceId l_id = node_to_face_matrix[node_id][l_face];
-                      if (primal_face_is_on_boundary[l_id]) {
-                        flux(face_id, i_face_cell) -=
-                          w_rl(node_id, l_face) *
-                          dot(M * dual_solution[face_dual_cell_id[l_id]], dual_solution[dual_cell_id]);
-                      }
-                    }
-                  }
-                }
-              }
-              //            }
-            }
-          }
-        }
-        // for (const auto& boundary_condition : boundary_condition_list) {
-        //   std::visit(
-        //     [&](auto&& bc) {
-        //       using T = std::decay_t<decltype(bc)>;
-        //       if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
-        //         const auto& face_list  = bc.faceList();
-        //         const auto& value_list = bc.valueList();
-        //         for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-        //           FaceId face_id      = face_list[i_face];
-        //           CellId dual_cell_id = face_dual_cell_id[face_id];
-        //           flux(face_id, 0)    = mes_l[face_id] * dot(value_list[i_face], dual_solution[dual_cell_id]); //
-        //           sign
-        //         }
-        //       }
-        //     },
-        //     boundary_condition);
-        // }
-        // for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        //   if (face_to_cell_matrix[face_id].size() == 2) {
-        //     CellId i_id = face_to_cell_matrix[face_id][0];
-        //     CellId j_id = face_to_cell_matrix[face_id][1];
-        //     if (flux(face_id, 0) != -flux(face_id, 1)) {
-        //       std::cout << "flux(" << i_id << "," << face_id << ")=" << flux(face_id, 0) << " not equal to
-        //       -flux("
-        //                 << j_id << "," << face_id << ")=" << -flux(face_id, 1) << "\n";
-        //     }
-        //   }
-        //   // exit(0);
-        // }
-        // Assemble
-        // CellValue<const TinyVector<Dimension>> fj = source->cellValues();
-
-        double sum_deltae = 0.;
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          m_energy_delta[cell_id] = 0.;   // dot(fj[cell_id], velocity[cell_id]);
-          sum_deltae += m_energy_delta[cell_id];
-        }
-        // CellValue<double>& deltae = m_energy_delta->cellValues();
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          for (size_t j = 0; j < face_to_cell_matrix[face_id].size(); j++) {
-            CellId i_id = face_to_cell_matrix[face_id][j];
-            m_energy_delta[i_id] -= flux(face_id, j) / primal_Vj[i_id];
-            sum_deltae -= flux(face_id, j);
-          }
-          // exit(0);
-        }
-
-        std::cout << "sum deltaej " << sum_deltae << "\n";
-      }
-    } else {
-      throw NotImplementedError("not done in 1d");
-    }
-    //    return m_energy_delta;
-  }
-};
-
-std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>>
-VectorDiamondSchemeHandler::apply() const
-{
-  return m_scheme->apply();
-}
-
-std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>>
-EnergyComputerHandler::computeEnergyUpdate() const
-{
-  return m_energy_computer->apply();
-}
-
-std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const ItemValueVariant>>
-VectorDiamondSchemeHandler::solution() const
-{
-  return m_scheme->getSolution();
-}
-
-std::shared_ptr<const DiscreteFunctionVariant>
-VectorDiamondSchemeHandler::dual_solution() const
-{
-  return m_scheme->getDualSolution();
-}
-
-VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
-  const std::shared_ptr<const DiscreteFunctionVariant>& alpha,
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_lambdab,
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_mub,
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_lambda,
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_mu,
-  const std::shared_ptr<const DiscreteFunctionVariant>& f,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-{
-  const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
-  if (not i_mesh) {
-    throw NormalError("primal discrete functions are not defined on the same mesh");
-  }
-  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_lambda, dual_lambdab, dual_mu, dual_mub});
-  if (not i_dual_mesh) {
-    throw NormalError("dual discrete functions are not defined on the same mesh");
-  }
-  checkDiscretizationType({alpha, dual_lambdab, dual_mub, dual_lambda, dual_mu, f}, DiscreteFunctionType::P0);
-
-  switch (i_mesh->dimension()) {
-  case 1: {
-    using MeshType                   = Mesh<Connectivity<1>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<1, const double>;
-    using DiscreteVectorFunctionType = DiscreteFunctionP0<1, const TinyVector<1>>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    m_scheme = std::make_unique<VectorDiamondScheme<1>>(mesh, alpha->get<DiscreteScalarFunctionType>(),
-                                                        dual_lambdab->get<DiscreteScalarFunctionType>(),
-                                                        dual_mub->get<DiscreteScalarFunctionType>(),
-                                                        dual_lambda->get<DiscreteScalarFunctionType>(),
-                                                        dual_mu->get<DiscreteScalarFunctionType>(),
-                                                        f->get<DiscreteVectorFunctionType>(), bc_descriptor_list);
-    break;
-  }
-  case 2: {
-    using MeshType                   = Mesh<Connectivity<2>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<2, const double>;
-    using DiscreteVectorFunctionType = DiscreteFunctionP0<2, const TinyVector<2>>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_scheme = std::make_unique<VectorDiamondScheme<2>>(mesh, alpha->get<DiscreteScalarFunctionType>(),
-                                                        dual_lambdab->get<DiscreteScalarFunctionType>(),
-                                                        dual_mub->get<DiscreteScalarFunctionType>(),
-                                                        dual_lambda->get<DiscreteScalarFunctionType>(),
-                                                        dual_mu->get<DiscreteScalarFunctionType>(),
-                                                        f->get<DiscreteVectorFunctionType>(), bc_descriptor_list);
-    break;
-  }
-  case 3: {
-    using MeshType                   = Mesh<Connectivity<3>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<3, const double>;
-    using DiscreteVectorFunctionType = DiscreteFunctionP0<3, const TinyVector<3>>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_scheme = std::make_unique<VectorDiamondScheme<3>>(mesh, alpha->get<DiscreteScalarFunctionType>(),
-                                                        dual_lambdab->get<DiscreteScalarFunctionType>(),
-                                                        dual_mub->get<DiscreteScalarFunctionType>(),
-                                                        dual_lambda->get<DiscreteScalarFunctionType>(),
-                                                        dual_mu->get<DiscreteScalarFunctionType>(),
-                                                        f->get<DiscreteVectorFunctionType>(), bc_descriptor_list);
-    break;
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
-VectorDiamondSchemeHandler::~VectorDiamondSchemeHandler() = default;
-
-EnergyComputerHandler::EnergyComputerHandler(
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_lambdab,
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_mub,
-  const std::shared_ptr<const DiscreteFunctionVariant>& U,
-  const std::shared_ptr<const DiscreteFunctionVariant>& dual_U,
-  const std::shared_ptr<const DiscreteFunctionVariant>& source,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-{
-  const std::shared_ptr i_mesh = getCommonMesh({U, source});
-  if (not i_mesh) {
-    throw NormalError("primal discrete functions are not defined on the same mesh");
-  }
-  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_lambdab, dual_mub, dual_U});
-  if (not i_dual_mesh) {
-    throw NormalError("dual discrete functions are not defined on the same mesh");
-  }
-  checkDiscretizationType({dual_lambdab, dual_mub, dual_U, U, source}, DiscreteFunctionType::P0);
-
-  switch (i_mesh->dimension()) {
-  case 1: {
-    using MeshType                   = Mesh<Connectivity<1>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<1, const double>;
-    using DiscreteVectorFunctionType = DiscreteFunctionP0<1, const TinyVector<1>>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_energy_computer =
-      std::make_unique<EnergyComputer<1>>(mesh, dual_lambdab->get<DiscreteScalarFunctionType>(),
-                                          dual_mub->get<DiscreteScalarFunctionType>(),
-                                          U->get<DiscreteVectorFunctionType>(),
-                                          dual_U->get<DiscreteVectorFunctionType>(),
-                                          source->get<DiscreteVectorFunctionType>(), bc_descriptor_list);
-    break;
-  }
-  case 2: {
-    using MeshType                   = Mesh<Connectivity<2>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<2, const double>;
-    using DiscreteVectorFunctionType = DiscreteFunctionP0<2, const TinyVector<2>>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_energy_computer =
-      std::make_unique<EnergyComputer<2>>(mesh, dual_lambdab->get<DiscreteScalarFunctionType>(),
-                                          dual_mub->get<DiscreteScalarFunctionType>(),
-                                          U->get<DiscreteVectorFunctionType>(),
-                                          dual_U->get<DiscreteVectorFunctionType>(),
-                                          source->get<DiscreteVectorFunctionType>(), bc_descriptor_list);
-
-    break;
-  }
-  case 3: {
-    using MeshType                   = Mesh<Connectivity<3>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<3, const double>;
-    using DiscreteVectorFunctionType = DiscreteFunctionP0<3, const TinyVector<3>>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_energy_computer =
-      std::make_unique<EnergyComputer<3>>(mesh, dual_lambdab->get<DiscreteScalarFunctionType>(),
-                                          dual_mub->get<DiscreteScalarFunctionType>(),
-                                          U->get<DiscreteVectorFunctionType>(),
-                                          dual_U->get<DiscreteVectorFunctionType>(),
-                                          source->get<DiscreteVectorFunctionType>(), bc_descriptor_list);
-
-    break;
-  }
-  default: {
-    throw UnexpectedError("invalid mesh dimension");
-  }
-  }
-}
-
-EnergyComputerHandler::~EnergyComputerHandler() = default;
diff --git a/src/scheme/VectorDiamondScheme.hpp b/src/scheme/VectorDiamondScheme.hpp
deleted file mode 100644
index 43213c94c80e0526e2b50471b948801d393b3775..0000000000000000000000000000000000000000
--- a/src/scheme/VectorDiamondScheme.hpp
+++ /dev/null
@@ -1,830 +0,0 @@
-#ifndef VECTOR_DIAMOND_SCHEME_HPP
-#define VECTOR_DIAMOND_SCHEME_HPP
-#include <algebra/CRSMatrix.hpp>
-#include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LinearSolver.hpp>
-#include <algebra/TinyVector.hpp>
-#include <algebra/Vector.hpp>
-#include <language/utils/InterpolateItemValue.hpp>
-#include <mesh/Connectivity.hpp>
-#include <mesh/DualConnectivityManager.hpp>
-#include <mesh/DualMeshManager.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
-#include <mesh/MeshFaceBoundary.hpp>
-#include <mesh/MeshFlatFaceBoundary.hpp>
-#include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
-#include <output/VTKWriter.hpp>
-#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
-#include <scheme/DiscreteFunctionVariant.hpp>
-#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-class ItemValueVariant;
-
-class VectorDiamondSchemeHandler
-{
- private:
-  class IVectorDiamondScheme;
-  template <size_t Dimension>
-  class VectorDiamondScheme;
-
-  template <size_t Dimension>
-  class InterpolationWeightsManager;
-
- public:
-  std::unique_ptr<IVectorDiamondScheme> m_scheme;
-
-  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const ItemValueVariant>> solution() const;
-
-  std::shared_ptr<const DiscreteFunctionVariant> dual_solution() const;
-
-  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> apply()
-    const;
-
-  VectorDiamondSchemeHandler(
-    const std::shared_ptr<const DiscreteFunctionVariant>& alphab,
-    const std::shared_ptr<const DiscreteFunctionVariant>& lambdab,
-    const std::shared_ptr<const DiscreteFunctionVariant>& alpha,
-    const std::shared_ptr<const DiscreteFunctionVariant>& lambda,
-    const std::shared_ptr<const DiscreteFunctionVariant>& mu,
-    const std::shared_ptr<const DiscreteFunctionVariant>& f,
-    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
-
-  ~VectorDiamondSchemeHandler();
-};
-
-class EnergyComputerHandler
-{
- private:
-  class IEnergyComputer;
-  template <size_t Dimension>
-  class EnergyComputer;
-
-  template <size_t Dimension>
-  class InterpolationWeightsManager;
-
- public:
-  std::unique_ptr<IEnergyComputer> m_energy_computer;
-  std::shared_ptr<const DiscreteFunctionVariant> dual_solution() const;
-
-  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>>
-  computeEnergyUpdate() const;
-
-  EnergyComputerHandler(const std::shared_ptr<const DiscreteFunctionVariant>& lambdab,
-                        const std::shared_ptr<const DiscreteFunctionVariant>& mub,
-                        const std::shared_ptr<const DiscreteFunctionVariant>& U,
-                        const std::shared_ptr<const DiscreteFunctionVariant>& dual_U,
-                        const std::shared_ptr<const DiscreteFunctionVariant>& source,
-                        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
-
-  ~EnergyComputerHandler();
-};
-
-template <size_t Dimension>
-class LegacyVectorDiamondScheme
-{
- private:
-  class DirichletBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const TinyVector<Dimension>>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    DirichletBoundaryCondition(const Array<const FaceId>& face_list,
-                               const Array<const TinyVector<Dimension>>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~DirichletBoundaryCondition() = default;
-  };
-
-  class NormalStrainBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-    const Array<const TinyVector<Dimension>>&
-    valueList() const
-    {
-      return m_value_list;
-    }
-
-    NormalStrainBoundaryCondition(const Array<const FaceId>& face_list,
-                                  const Array<const TinyVector<Dimension>>& value_list)
-      : m_value_list{value_list}, m_face_list{face_list}
-    {
-      Assert(m_value_list.size() == m_face_list.size());
-    }
-
-    ~NormalStrainBoundaryCondition() = default;
-  };
-
-  class SymmetryBoundaryCondition
-  {
-   private:
-    const Array<const TinyVector<Dimension>> m_value_list;
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-   public:
-    SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
-
-    ~SymmetryBoundaryCondition() = default;
-  };
-
- public:
-  LegacyVectorDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
-                            const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-                            const FunctionSymbolId& lambda_id,
-                            const FunctionSymbolId& mu_id,
-                            const FunctionSymbolId& f_id,
-                            CellValue<TinyVector<Dimension>>& Uj,
-                            FaceValue<TinyVector<Dimension>>& Ul,
-                            const double& Tf,
-                            const double& dt,
-                            CellValuePerNode<double>& w_rj,
-                            FaceValuePerNode<double>& w_rl)
-  {
-    using ConnectivityType = Connectivity<Dimension>;
-    using MeshType         = Mesh<ConnectivityType>;
-    using MeshDataType     = MeshData<Dimension>;
-
-    using BoundaryCondition =
-      std::variant<DirichletBoundaryCondition, NormalStrainBoundaryCondition, SymmetryBoundaryCondition>;
-
-    using BoundaryConditionList = std::vector<BoundaryCondition>;
-
-    BoundaryConditionList boundary_condition_list;
-
-    std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    NodeValue<bool> is_dirichlet{mesh->connectivity()};
-    is_dirichlet.fill(false);
-    NodeValue<TinyVector<Dimension>> dirichlet_value{mesh->connectivity()};
-    {
-      TinyVector<Dimension> nan_tiny_vector;
-      for (size_t i = 0; i < Dimension; ++i) {
-        nan_tiny_vector[i] = std::numeric_limits<double>::signaling_NaN();
-      }
-      dirichlet_value.fill(nan_tiny_vector);
-    }
-
-    for (const auto& bc_descriptor : bc_descriptor_list) {
-      bool is_valid_boundary_condition = true;
-
-      switch (bc_descriptor->type()) {
-      case IBoundaryConditionDescriptor::Type::symmetry: {
-        const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-          dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-
-        if constexpr (Dimension > 1) {
-          MeshFlatFaceBoundary<Dimension> mesh_face_boundary =
-            getMeshFlatFaceBoundary(*mesh, sym_bc_descriptor.boundaryDescriptor());
-          boundary_condition_list.push_back(SymmetryBoundaryCondition{mesh_face_boundary.faceList()});
-        } else {
-          throw NotImplementedError("Symmetry conditions are not supported in 1d");
-        }
-
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::dirichlet: {
-        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
-          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
-        if (dirichlet_bc_descriptor.name() == "dirichlet") {
-          if constexpr (Dimension > 1) {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-            const FunctionSymbolId g_id                   = dirichlet_bc_descriptor.rhsSymbolId();
-            Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-              TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(),
-                                                                            mesh_face_boundary.faceList());
-            boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-          } else {
-            throw NotImplementedError("Neumann conditions are not supported in 1d");
-          }
-
-        } else if (dirichlet_bc_descriptor.name() == "normal_strain") {
-          if constexpr (Dimension > 1) {
-            MeshFaceBoundary<Dimension> mesh_face_boundary =
-              getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-            const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-
-            Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-              TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(),
-                                                                            mesh_face_boundary.faceList());
-            boundary_condition_list.push_back(NormalStrainBoundaryCondition{mesh_face_boundary.faceList(), value_list});
-
-          } else {
-            throw NotImplementedError("Normal strain conditions are not supported in 1d");
-          }
-
-        } else {
-          is_valid_boundary_condition = false;
-        }
-        break;
-      }
-      default: {
-        is_valid_boundary_condition = false;
-      }
-      }
-      if (not is_valid_boundary_condition) {
-        std::ostringstream error_msg;
-        error_msg << *bc_descriptor << " is an invalid boundary condition for elasticity equation";
-        throw NormalError(error_msg.str());
-      }
-    }
-
-    if constexpr (Dimension > 1) {
-      const CellValue<const size_t> cell_dof_number = [&] {
-        CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
-        return compute_cell_dof_number;
-      }();
-      size_t number_of_dof = mesh->numberOfCells();
-
-      const FaceValue<const size_t> face_dof_number = [&] {
-        FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
-        compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>) or
-                            (std::is_same_v<T, DirichletBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id = face_list[i_face];
-                  if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
-                    std::ostringstream os;
-                    os << "The face " << face_id << " is used at least twice for boundary conditions";
-                    throw NormalError(os.str());
-                  } else {
-                    compute_face_dof_number[face_id] = number_of_dof++;
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return compute_face_dof_number;
-      }();
-
-      const auto& primal_face_to_node_matrix             = mesh->connectivity().faceToNodeMatrix();
-      const auto& face_to_cell_matrix                    = mesh->connectivity().faceToCellMatrix();
-      const FaceValue<const bool> primal_face_is_neumann = [&] {
-        FaceValue<bool> face_is_neumann{mesh->connectivity()};
-        face_is_neumann.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id     = face_list[i_face];
-                  face_is_neumann[face_id] = true;
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return face_is_neumann;
-      }();
-
-      const FaceValue<const bool> primal_face_is_symmetry = [&] {
-        FaceValue<bool> face_is_symmetry{mesh->connectivity()};
-        face_is_symmetry.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, SymmetryBoundaryCondition>)) {
-                const auto& face_list = bc.faceList();
-
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id      = face_list[i_face];
-                  face_is_symmetry[face_id] = true;
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        return face_is_symmetry;
-      }();
-
-      NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
-      }
-
-      primal_node_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-            NodeId node_id                      = primal_face_to_node_matrix[face_id][i_node];
-            primal_node_is_on_boundary[node_id] = true;
-          }
-        }
-      }
-
-      FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
-      }
-
-      primal_face_is_on_boundary.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        if (face_to_cell_matrix[face_id].size() == 1) {
-          primal_face_is_on_boundary[face_id] = true;
-        }
-      }
-
-      FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
-      if (parallel::size() > 1) {
-        throw NotImplementedError("Calculation of face_is_neumann is incorrect");
-      }
-
-      primal_face_is_dirichlet.fill(false);
-      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-        primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] &&
-                                             (!primal_face_is_neumann[face_id]) && (!primal_face_is_symmetry[face_id]));
-      }
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-      const FaceValue<const TinyVector<Dimension>>& xl               = mesh_data.xl();
-      const CellValue<const TinyVector<Dimension>>& xj               = mesh_data.xj();
-      const auto& node_to_cell_matrix                                = mesh->connectivity().nodeToCellMatrix();
-      const auto& node_to_face_matrix                                = mesh->connectivity().nodeToFaceMatrix();
-      const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
-
-      {
-        std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
-
-        MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
-
-        std::shared_ptr mapper =
-          DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
-
-        CellValue<double> dual_muj =
-          InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(mu_id,
-                                                                                                    diamond_mesh_data
-                                                                                                      .xj());
-
-        CellValue<double> dual_lambdaj =
-          InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(lambda_id,
-                                                                                                    diamond_mesh_data
-                                                                                                      .xj());
-
-        CellValue<TinyVector<Dimension>> fj = InterpolateItemValue<TinyVector<Dimension>(
-          TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
-
-        const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
-
-        const FaceValue<const double> mes_l = [&] {
-          if constexpr (Dimension == 1) {
-            FaceValue<double> compute_mes_l{mesh->connectivity()};
-            compute_mes_l.fill(1);
-            return compute_mes_l;
-          } else {
-            return mesh_data.ll();
-          }
-        }();
-
-        const CellValue<const double> dual_mes_l_j = [=] {
-          CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
-          mapper->toDualCell(mes_l, compute_mes_j);
-
-          return compute_mes_j;
-        }();
-
-        const CellValue<const double> primal_Vj   = mesh_data.Vj();
-        FaceValue<const CellId> face_dual_cell_id = [=]() {
-          FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
-          CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
-
-          mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
-
-          return computed_face_dual_cell_id;
-        }();
-
-        NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
-          CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
-          cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> node_primal_id{mesh->connectivity()};
-
-          parallel_for(
-            mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
-
-          NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
-
-          mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
-
-          return computed_dual_node_primal_node_id;
-        }();
-
-        CellValue<NodeId> primal_cell_dual_node_id = [=]() {
-          CellValue<NodeId> cell_id{mesh->connectivity()};
-          NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
-          node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
-          NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
-
-          CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
-
-          mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
-
-          return cell_id;
-        }();
-        const auto& dual_Cjr                     = diamond_mesh_data.Cjr();
-        FaceValue<TinyVector<Dimension>> dualClj = [&] {
-          FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
-          const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
-          const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
-          parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-              const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-              for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
-                CellId cell_id            = primal_face_to_cell[i];
-                const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
-                for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
-                     i_dual_cell++) {
-                  const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
-                  if (face_dual_cell_id[face_id] == dual_cell_id) {
-                    for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
-                         i_dual_node++) {
-                      const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                      if (final_dual_node_id == dual_node_id) {
-                        computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
-                      }
-                    }
-                  }
-                }
-              }
-            });
-          return computedClj;
-        }();
-
-        FaceValue<TinyVector<Dimension>> nlj = [&] {
-          FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
-          parallel_for(
-            mesh->numberOfFaces(),
-            PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
-          return computedNlj;
-        }();
-
-        FaceValue<const double> alpha_lambda_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_lambdaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        FaceValue<const double> alpha_mu_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
-          parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_muj[diamond_cell_id] / dual_Vj[diamond_cell_id];
-            });
-
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
-        }();
-
-        double unsteady    = (Tf == 0) ? 0 : 1;
-        double time_factor = (unsteady == 0) ? 1 : dt;
-
-        TinyMatrix<Dimension> I = identity;
-
-        const Array<int> non_zeros{number_of_dof * Dimension};
-        non_zeros.fill(Dimension * Dimension);
-        CRSMatrixDescriptor<double> S(number_of_dof * Dimension, number_of_dof * Dimension, non_zeros);
-
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const double beta_mu_l          = l2Norm(dualClj[face_id]) * alpha_mu_l[face_id] * mes_l[face_id];
-          const double beta_lambda_l      = l2Norm(dualClj[face_id]) * alpha_lambda_l[face_id] * mes_l[face_id];
-          const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
-          for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
-            const CellId i_id                      = primal_face_to_cell[i_cell];
-            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
-
-            const TinyVector<Dimension> nil = [&] {
-              if (is_face_reversed_for_cell_i) {
-                return -nlj[face_id];
-              } else {
-                return nlj[face_id];
-              }
-            }();
-            for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
-              const CellId j_id = primal_face_to_cell[j_cell];
-              TinyMatrix<Dimension> M =
-                beta_mu_l * I + beta_mu_l * tensorProduct(nil, nil) + beta_lambda_l * tensorProduct(nil, nil);
-              TinyMatrix<Dimension> N = tensorProduct(nil, nil);
-
-              if (i_cell == j_cell) {
-                for (size_t i = 0; i < Dimension; ++i) {
-                  for (size_t j = 0; j < Dimension; ++j) {
-                    S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) +=
-                      (time_factor * M(i, j) + unsteady * primal_Vj[i_id]);
-                    if (primal_face_is_neumann[face_id]) {
-                      S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) -= M(i, j);
-                    }
-                    if (primal_face_is_symmetry[face_id]) {
-                      S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
-                        -((i == j) ? 1 : 0) + N(i, j);
-                      S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + j) +=
-                        (i == j) ? 1 : 0;
-                    }
-                  }
-                }
-              } else {
-                for (size_t i = 0; i < Dimension; ++i) {
-                  for (size_t j = 0; j < Dimension; ++j) {
-                    S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -=
-                      time_factor * M(i, j);
-                  }
-                }
-              }
-            }
-          }
-        }
-
-        const auto& dual_cell_to_node_matrix   = diamond_mesh->connectivity().cellToNodeMatrix();
-        const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const double alpha_mu_face_id     = mes_l[face_id] * alpha_mu_l[face_id];
-          const double alpha_lambda_face_id = mes_l[face_id] * alpha_lambda_l[face_id];
-
-          for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
-            CellId i_id                            = face_to_cell_matrix[face_id][i_face_cell];
-            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
-
-            for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
-              NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
-
-              const TinyVector<Dimension> nil = [&] {
-                if (is_face_reversed_for_cell_i) {
-                  return -nlj[face_id];
-                } else {
-                  return nlj[face_id];
-                }
-              }();
-
-              CellId dual_cell_id = face_dual_cell_id[face_id];
-
-              for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
-                const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                if (dual_node_primal_node_id[dual_node_id] == node_id) {
-                  const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-
-                  TinyMatrix<Dimension> M = alpha_mu_face_id * dot(Clr, nil) * I +
-                                            alpha_mu_face_id * tensorProduct(Clr, nil) +
-                                            alpha_lambda_face_id * tensorProduct(nil, Clr);
-
-                  for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
-                    CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
-                    for (size_t i = 0; i < Dimension; ++i) {
-                      for (size_t j = 0; j < Dimension; ++j) {
-                        S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -=
-                          time_factor * w_rj(node_id, j_cell) * M(i, j);
-                        if (primal_face_is_neumann[face_id]) {
-                          S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
-                            w_rj(node_id, j_cell) * M(i, j);
-                        }
-                      }
-                    }
-                  }
-                  if (primal_node_is_on_boundary[node_id]) {
-                    for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
-                      FaceId l_id = node_to_face_matrix[node_id][l_face];
-                      if (primal_face_is_on_boundary[l_id]) {
-                        for (size_t i = 0; i < Dimension; ++i) {
-                          for (size_t j = 0; j < Dimension; ++j) {
-                            S(cell_dof_number[i_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) -=
-                              time_factor * w_rl(node_id, l_face) * M(i, j);
-                          }
-                        }
-                        if (primal_face_is_neumann[face_id]) {
-                          for (size_t i = 0; i < Dimension; ++i) {
-                            for (size_t j = 0; j < Dimension; ++j) {
-                              S(face_dof_number[face_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) +=
-                                w_rl(node_id, l_face) * M(i, j);
-                            }
-                          }
-                        }
-                      }
-                    }
-                  }
-                }
-              }
-              //            }
-            }
-          }
-        }
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          if (primal_face_is_dirichlet[face_id]) {
-            for (size_t i = 0; i < Dimension; ++i) {
-              S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + i) += 1;
-            }
-          }
-        }
-
-        Vector<double> b{number_of_dof * Dimension};
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          for (size_t i = 0; i < Dimension; ++i) {
-            b[(cell_dof_number[cell_id] * Dimension) + i] =
-              primal_Vj[cell_id] * (time_factor * fj[cell_id][i] + unsteady * Uj[cell_id][i]);
-          }
-        }
-
-        // Dirichlet
-        NodeValue<bool> node_tag{mesh->connectivity()};
-        node_tag.fill(false);
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
-                const auto& face_list  = bc.faceList();
-                const auto& value_list = bc.valueList();
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  const FaceId face_id = face_list[i_face];
-
-                  for (size_t i = 0; i < Dimension; ++i) {
-                    b[(face_dof_number[face_id] * Dimension) + i] += value_list[i_face][i];
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
-                const auto& face_list  = bc.faceList();
-                const auto& value_list = bc.valueList();
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  FaceId face_id = face_list[i_face];
-                  for (size_t i = 0; i < Dimension; ++i) {
-                    b[face_dof_number[face_id] * Dimension + i] += mes_l[face_id] * value_list[i_face][i];   // sign
-                  }
-                }
-              }
-            },
-            boundary_condition);
-        }
-
-        CRSMatrix A{S.getCRSMatrix()};
-        Vector<double> U{number_of_dof * Dimension};
-        U        = zero;
-        Vector r = A * U - b;
-        std::cout << "initial (real) residu = " << std::sqrt(dot(r, r)) << '\n';
-
-        LinearSolver solver;
-        solver.solveLocalSystem(A, U, b);
-
-        r = A * U - b;
-
-        std::cout << "final (real) residu = " << std::sqrt(dot(r, r)) << '\n';
-
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          for (size_t i = 0; i < Dimension; ++i) {
-            Uj[cell_id][i] = U[(cell_dof_number[cell_id] * Dimension) + i];
-          }
-        }
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          for (size_t i = 0; i < Dimension; ++i) {
-            if (primal_face_is_on_boundary[face_id]) {
-              Ul[face_id][i] = U[(face_dof_number[face_id] * Dimension) + i];
-            }
-          }
-        }
-        NodeValue<TinyVector<3>> ur3d{mesh->connectivity()};
-        ur3d.fill(zero);
-
-        parallel_for(
-          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-            TinyVector<Dimension> x = zero;
-            const auto node_cells   = node_to_cell_matrix[node_id];
-            for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
-              CellId cell_id = node_cells[i_cell];
-              x += w_rj(node_id, i_cell) * Uj[cell_id];
-            }
-            const auto node_faces = node_to_face_matrix[node_id];
-            for (size_t i_face = 0; i_face < node_faces.size(); ++i_face) {
-              FaceId face_id = node_faces[i_face];
-              if (primal_face_is_on_boundary[face_id]) {
-                x += w_rl(node_id, i_face) * Ul[face_id];
-              }
-            }
-            for (size_t i = 0; i < Dimension; ++i) {
-              ur3d[node_id][i] = x[i];
-            }
-          });
-      }
-    } else {
-      throw NotImplementedError("not done in 1d");
-    }
-  }
-};
-template LegacyVectorDiamondScheme<1>::LegacyVectorDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  CellValue<TinyVector<1>>&,
-  FaceValue<TinyVector<1>>&,
-  const double&,
-  const double&,
-  CellValuePerNode<double>&,
-  FaceValuePerNode<double>&);
-
-template LegacyVectorDiamondScheme<2>::LegacyVectorDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  CellValue<TinyVector<2>>&,
-  FaceValue<TinyVector<2>>&,
-  const double&,
-  const double&,
-  CellValuePerNode<double>&,
-  FaceValuePerNode<double>&);
-
-template LegacyVectorDiamondScheme<3>::LegacyVectorDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  CellValue<TinyVector<3>>&,
-  FaceValue<TinyVector<3>>&,
-  const double&,
-  const double&,
-  CellValuePerNode<double>&,
-  FaceValuePerNode<double>&);
-
-#endif   // VECTOR_DIAMOND_SCHEME_HPP
diff --git a/src/utils/CouplingData.hpp b/src/utils/CouplingData.hpp
index 485fe17b6a8b291def84512cfc8e40d84cb6e7d1..651b1539bbdc55a69b6a2f6ddb25050bb55cd449 100644
--- a/src/utils/CouplingData.hpp
+++ b/src/utils/CouplingData.hpp
@@ -16,11 +16,12 @@
 
 class IMesh;
 
-template <size_t Dimension>
+template <MeshConcept MeshType>
 class CouplingData
 {
  private:
-  using MeshType = Mesh<Connectivity<Dimension>>;
+  // using MeshType = Mesh<Connectivity<Dimension>>;
+  static constexpr size_t Dimension = MeshType::Dimension;
   using Rd       = TinyVector<Dimension>;
   class CouplingBoundaryCondition;
 
@@ -36,7 +37,7 @@ class CouplingData
   static auto&
   getInstanceOrBuildIt()
   {
-    static CouplingData<Dimension> cdata;
+    static CouplingData<MeshType> cdata;
     return cdata;
   }
   std::tuple<const std::shared_ptr<const ItemValueVariant>,
@@ -81,11 +82,11 @@ class CouplingData
   CouplingData();
 };
 
-template <size_t Dimension>
-class CouplingData<Dimension>::CouplingBoundaryCondition
+template <MeshConcept MeshType>
+class CouplingData<MeshType>::CouplingBoundaryCondition
 {
  private:
-  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const MeshNodeBoundary m_mesh_node_boundary;
 
  public:
   const Array<const NodeId>&
@@ -94,7 +95,7 @@ class CouplingData<Dimension>::CouplingBoundaryCondition
     return m_mesh_node_boundary.nodeList();
   }
 
-  CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary)
+  CouplingBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary)
     : m_mesh_node_boundary{mesh_node_boundary}
   {
     ;
@@ -103,9 +104,9 @@ class CouplingData<Dimension>::CouplingBoundaryCondition
   ~CouplingBoundaryCondition() = default;
 };
 
-template <size_t Dimension>
+template <MeshConcept MeshType>
 std::shared_ptr<const IMesh>
-CouplingData<Dimension>::update_mesh(const std::shared_ptr<const IMesh>& i_mesh,
+CouplingData<MeshType>::update_mesh(const std::shared_ptr<const IMesh>& i_mesh,
                                      const std::shared_ptr<const IBoundaryDescriptor>& boundary,
                                      const int64_t src,
                                      const int64_t tag)
@@ -120,7 +121,7 @@ CouplingData<Dimension>::update_mesh(const std::shared_ptr<const IMesh>& i_mesh,
 
   /* const std::shared_ptr given_mesh               = std::dynamic_pointer_cast<const MeshType>(*i_mesh); */
   const MeshType& mesh                           = dynamic_cast<const MeshType&>(*i_mesh);
-  MeshNodeBoundary<Dimension> mesh_node_boundary = getMeshNodeBoundary(mesh, *boundary);
+  MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, *boundary);
   const auto node_list                           = mesh_node_boundary.nodeList();
 
   NodeValue<Rd> new_xr = copy(mesh.xr());
@@ -140,12 +141,12 @@ CouplingData<Dimension>::update_mesh(const std::shared_ptr<const IMesh>& i_mesh,
   return std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
 }
 
-template <size_t Dimension>
+template <MeshConcept MeshType>
 std::tuple<const std::shared_ptr<const ItemValueVariant>,
            const std::shared_ptr<const SubItemValuePerItemVariant>,
            const std::shared_ptr<const ItemValueVariant>,
            const std::shared_ptr<const SubItemValuePerItemVariant>>
-CouplingData<Dimension>::get_flux_and_velocity(
+CouplingData<MeshType>::get_flux_and_velocity(
   const std::shared_ptr<const IMesh>& i_mesh,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
   const std::shared_ptr<const ItemValueVariant>& ur,
@@ -221,13 +222,14 @@ CouplingData<Dimension>::get_flux_and_velocity(
         }
 
         if (costo->ctm == costo->m) {
-          /* Serializer<Dimension>().BC(std::make_shared<const IBoundaryDescriptor>( */
+	  std::cout<<"ici"<<std::endl;
+          /* Serializer<MeshType>().BC(std::make_shared<const IBoundaryDescriptor>( */
           /*                              bc_descriptor_list[0]->boundaryDescriptor()), */
           /*                            i_mesh, 0, 2, 800); */
           /* auto ptr = std::make_shared<const IBoundaryDescriptor>(bc_descriptor_list[0]->boundaryDescriptor()); */
           /* Serializer<Dimension>().BC(ptr, i_mesh, 0, 2, 800); */
 
-          auto costo = parallel::Messenger::getInstance().myCoupling;
+          // auto costo = parallel::Messenger::getInstance().myCoupling;
           /* std::vector<int> shape; */
           /* std::vector<double> pts; */
           /* const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); */
@@ -268,8 +270,8 @@ CouplingData<Dimension>::get_flux_and_velocity(
           /* costo->sendData(shape, pts, 2, 800); */
           /* auto ptr = std::make_shared<const IBoundaryDescriptor>(bc_descriptor_list[0]->boundaryDescriptor()); */
 
-          Serializer<Dimension>().BC(bc_descriptor_list[0]->boundaryDescriptor(), i_mesh, 0, 2, 800);
-          costo->ctm = 0;
+          // Serializer<MeshType>().BC(bc_descriptor_list[0]->boundaryDescriptor(), i_mesh, 0, 2, 800);
+          // costo->ctm = 0;
         }
         return std::make_tuple(std::make_shared<const ItemValueVariant>(ur_ok),
                                std::make_shared<const SubItemValuePerItemVariant>(Fjr_ok),
@@ -280,15 +282,15 @@ CouplingData<Dimension>::get_flux_and_velocity(
   }
 }
 
-template <size_t Dimension>
-CouplingData<Dimension>::CouplingData()
+template <MeshConcept MeshType>
+CouplingData<MeshType>::CouplingData()
 {
   std::cout << "\033[01;31m"
-            << "COUPLING DATA CREATION" << Dimension << "\033[00;00m" << std::endl;
+            << "COUPLING DATA CREATION" << MeshType::Dimension << "\033[00;00m" << std::endl;
 }
 
-using CouplingData3D = CouplingData<3>;
-using CouplingData2D = CouplingData<2>;
-using CouplingData1D = CouplingData<1>;
+using CouplingData3D = CouplingData<Mesh<1>>;
+using CouplingData2D = CouplingData<Mesh<2>>;
+using CouplingData1D = CouplingData<Mesh<3>>;
 
 #endif   // PUGS_HAS_COSTO
diff --git a/src/utils/Serializer.cpp b/src/utils/Serializer.cpp
index 1bf1e080ba8a5bfe5c8c31268ae05353e68c4f2e..9227e2ec7c3e629afce9e4afb9d56c0549fcfcff 100644
--- a/src/utils/Serializer.cpp
+++ b/src/utils/Serializer.cpp
@@ -2,4 +2,334 @@
 
 #ifdef PUGS_HAS_COSTO
 
+// MESH node position
+void
+Serializer::mesh(const std::shared_ptr<const MeshVariant> mesh_v,
+		 const int64_t dest,
+		 const int64_t tag)
+{
+  if (not mesh_v) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  std::visit(
+    [&](auto&& p_mesh) {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      static const auto Dimension = MeshType::Dimension;
+      using Rd           = TinyVector<Dimension>;
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+	std::cout<<"Serialize polygonal mesh, dim:"<<Dimension<<std::endl;
+
+
+	auto costo = parallel::Messenger::getInstance().myCoupling;
+	std::vector<int> shape;
+	std::vector<double> pts;
+
+	const std::shared_ptr mesh    = std::dynamic_pointer_cast<const MeshType>(p_mesh);
+	const NodeValue<const Rd>& xr = mesh->xr();
+
+	const auto& connectivity = mesh->connectivity();
+	const auto node_is_owned = connectivity.nodeIsOwned().arrayView();
+
+	Array<TinyVector<Dimension>> positions(connectivity.numberOfNodes());
+	size_t sz = 0;
+	for (NodeId r = 0; r < node_is_owned.size(); ++r) {
+	  if (node_is_owned[r]){
+	    for (unsigned short i = 0; i < Dimension; ++i) {
+	      positions[sz][i] = xr[r][i];
+	    }
+	    sz += 1;
+	  }
+	}
+
+	/*TODO: serializer directement position pour eviter une copie*/
+	pts.resize(sz * Dimension);
+	for (unsigned short r = 0; r < sz; ++r) {
+	    for (unsigned short j = 0; j < Dimension; ++j) {
+	      pts[Dimension * r + j] = positions[r][j];
+	  }
+	}
+	shape.resize(3);
+	shape[0] = sz;
+	shape[1] = Dimension;
+	shape[2] = 0;
+	costo->sendData(shape, pts, dest, tag);
+
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    mesh_v->variant());
+}
+
+void Serializer::interfaces(const std::shared_ptr<const MeshVariant> mesh_v,
+			    const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
+			    const int64_t dest,
+			    const int64_t tag)
+{
+  if (not mesh_v) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  std::visit(
+    [&](auto&& p_mesh) {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      static const auto Dimension = MeshType::Dimension;
+      using Rd           = TinyVector<Dimension>;
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+	std::cout<<"Serialize polygonal mesh, dim:"<<Dimension<<std::endl;
+
+
+	auto costo = parallel::Messenger::getInstance().myCoupling;
+	std::vector<int> shape;
+	std::vector<double> pts;
+
+	// const std::shared_ptr mesh    = std::dynamic_pointer_cast<const MeshType>(p_mesh);
+	// const NodeValue<const Rd>& xr = mesh->xr();
+
+	// const auto& connectivity = mesh->connectivity();
+	// const auto node_is_owned = connectivity.nodeIsOwned().arrayView();
+
+	// // const auto node_is_owned_generic = connectivity.isOwned<ItemType::node>();
+	// // const auto node_is_owned = connectivity.nodeIsOwned();
+
+
+	// Array<TinyVector<Dimension>> positions(connectivity.numberOfNodes());
+	// size_t sz = 0;
+	// for (NodeId r = 0; r < node_is_owned.size(); ++r) {
+	//   if (node_is_owned[r]){
+	//     for (unsigned short i = 0; i < Dimension; ++i) {
+	//       positions[sz][i] = xr[r][i];
+	//     }
+	//     sz += 1;
+	//   }
+	// }
+
+	// /*TODO: serializer directement position pour eviter une copie*/
+	// pts.resize(sz * Dimension);
+	// for (unsigned short r = 0; r < sz; ++r) {
+	//     for (unsigned short j = 0; j < Dimension; ++j) {
+	//       pts[Dimension * r + j] = positions[r][j];
+	//   }
+	// }
+	// shape.resize(3);
+	// shape[0] = sz;
+	// shape[1] = Dimension;
+	// shape[2] = 0;
+	// costo->sendData(shape, pts, dest, tag);
+
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    mesh_v->variant());
+}
+
+
+
+// BC nodes position
+void
+Serializer::BC(const IBoundaryDescriptor& boundary,
+	       const std::shared_ptr<const MeshVariant> mesh_v,
+	       const int location,
+	       const int64_t dest,
+	       const int64_t tag)
+{
+  if (not mesh_v) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  std::visit(
+    [&](auto&& p_mesh) {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      static const auto Dimension = MeshType::Dimension;
+      using Rd           = TinyVector<Dimension>;
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+	std::cout<<"Serialize interface of polygonal mesh, dim:"<<Dimension<<std::endl;
+
+	auto costo = parallel::Messenger::getInstance().myCoupling;
+	std::vector<int> shape;
+	std::vector<double> pts;
+	std::cout << "\033[01;31m"
+		  << "L" << location << " ff " << at_face << "\033[00;00m" << std::endl;
+	const std::shared_ptr mesh    = std::dynamic_pointer_cast<const MeshType>(p_mesh);
+
+	if (location == at_face) {
+	  std::cout << "\033[01;31m"
+		    << "ICI AT FACE"
+		    << "\033[00;00m" << std::endl;
+
+	  // const MeshType& mesh_data   = *mesh->get<MeshType>();
+
+// 	  // const auto xl = mesh->xl();
+// 	  const auto xl = MeshDataManager::instance().getMeshData(mesh_data).xl();
+// // ::	  // const NodeValue<const Rd>& xl = mesh_data->xl();
+// // 	  // const auto xl = mesh_data.xl();
+
+// 	  const auto& connectivity = mesh->connectivity();
+// 	  const auto face_is_owned = connectivity.faceIsOwned().arrayView();
+
+// 	  Array<TinyVector<Dimension>> fpositions(connectivity.numberOfFaces());
+
+// 	  size_t sz = 0;
+// 	  for (FaceId r = 0; r < face_is_owned.size(); ++r) {
+// 	    if (face_is_owned[r]){
+// 	      for (unsigned short i = 0; i < Dimension; ++i) {
+// 		fpositions[sz][i] = xl[r][i];
+// 	      }
+// 	      sz += 1;
+// 	    }
+// 	  }
+
+
+// 	  /*TODO: serializer directement fposition pour eviter une copie*/
+// 	  MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(*mesh, boundary);
+// 	  /* mesh_face_boundary.faceList() */
+// 	  const auto face_list = mesh_face_boundary.faceList();
+// 	  pts.resize(face_list.size() * Dimension);
+
+// 	  for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+// 	    const FaceId face_id = face_list[i_face];
+// 	    if (face_is_owned[face_id]){
+
+// 	      for (unsigned short j = 0; j < Dimension; ++j) {
+// 		pts[Dimension * i_face + j] = fpositions[face_id][j];
+// 	      }
+
+// 	    }
+// 	  }
+
+// 	  shape.resize(3);
+// 	  shape[0] = face_list.size();
+// 	  shape[1] = Dimension;
+// 	  shape[2] = 0;
+// 	  costo->sendData(shape, pts, dest, tag);
+
+
+
+	} else if (location == at_node) {
+	  std::cout << "\033[01;31m"
+		    << "ICI AT node"
+		    << "\033[00;00m" << std::endl;
+
+	  const auto xr = mesh->xr();
+
+	  const auto& connectivity = mesh->connectivity();
+	  const auto node_is_owned = connectivity.nodeIsOwned().arrayView();
+
+	  MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(*mesh, boundary);
+	  const auto node_list = mesh_node_boundary.nodeList();
+
+	  Array<TinyVector<Dimension>> positions(node_list.size());
+
+	  size_t sz = 0;
+	  for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+	    const NodeId r = node_list[i_node];
+	    if (node_is_owned[r]){
+	      for (unsigned short i = 0; i < Dimension; ++i) {
+		positions[sz][i] = xr[r][i];
+	      }
+	      sz += 1;
+	    }
+	  }
+
+
+	  /*TODO: serializer directement position pour eviter une copie*/
+	  pts.resize(sz * Dimension);
+	  for (unsigned short r = 0; r < sz; ++r) {
+	    for (unsigned short j = 0; j < Dimension; ++j) {
+	      pts[Dimension * r + j] = positions[r][j];
+	    }
+	  }
+	  shape.resize(3);
+	  shape[0] = sz;
+	  shape[1] = Dimension;
+	  shape[2] = 0;
+	  costo->sendData(shape, pts, dest, tag);
+
+	} else {
+	  throw UnexpectedError(" location must be 0 or 1");
+	}
+
+
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    mesh_v->variant());
+}
+
+
+// Field at node
+void Serializer::BC_field(const std::shared_ptr<const MeshVariant> mesh_v,
+			  const std::shared_ptr<const IBoundaryDescriptor>& boundary,
+			  const std::shared_ptr<const ItemValueVariant>& field,
+			  const int64_t dest,
+			  const int64_t tag)
+{
+  if (not mesh_v) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  std::visit(
+    [&](auto&& p_mesh) {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      static const auto Dimension = MeshType::Dimension;
+      using Rd           = TinyVector<Dimension>;
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+	std::cout<<"Serialize polygonal field, dim:"<<Dimension<<std::endl;
+
+	// const std::shared_ptr i_mesh = getCommonMesh({mesh, field});
+	// if (not i_mesh) {
+	//   throw NormalError("primal discrete functions are not defined on the same mesh");
+	// }
+
+	auto costo = parallel::Messenger::getInstance().myCoupling;
+	std::vector<int> shape;
+	std::vector<double> data;
+
+	const std::shared_ptr mesh    = std::dynamic_pointer_cast<const MeshType>(p_mesh);
+	const NodeValue<const Rd>& xr = mesh->xr();
+
+	const auto& connectivity = mesh->connectivity();
+	const auto node_is_owned = connectivity.nodeIsOwned().arrayView();
+
+	const NodeValue<const TinyVector<Dimension>> nField = field->get<NodeValue<const TinyVector<Dimension>>>();
+	const auto node_list = getMeshNodeBoundary(*mesh, *boundary).nodeList();
+
+	Array<TinyVector<Dimension>> positions(node_list.size());
+
+
+	size_t sz = 0;
+	for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+	  const NodeId r = node_list[i_node];
+	  if (node_is_owned[r]){
+	    for (unsigned short i = 0; i < Dimension; ++i) {
+	      positions[sz][i] = nField[r][i];
+	    }
+	    sz += 1;
+	  }
+	}
+
+	/*TODO: serializer directement position pour eviter une copie*/
+	data.resize(sz * Dimension);
+	for (unsigned short r = 0; r < sz; ++r) {
+	    for (unsigned short j = 0; j < Dimension; ++j) {
+	      data[Dimension * r + j] = positions[r][j];
+	    }
+	}
+
+	shape.resize(3);
+	shape[0] = sz;
+	shape[1] = Dimension;
+	shape[2] = 0;
+	costo->sendData(shape, data, dest, tag);
+
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    mesh_v->variant());
+}
+
 #endif   // PUGS_HAS_COSTO
diff --git a/src/utils/Serializer.hpp b/src/utils/Serializer.hpp
index 7982581d1087be53abdc84abd96a28d2306d71bf..0060afd7937be5862d2bb4bfe4561b39c3644615 100644
--- a/src/utils/Serializer.hpp
+++ b/src/utils/Serializer.hpp
@@ -1,5 +1,6 @@
 #pragma once
 #include <utils/pugs_config.hpp>
+#include <cstdint>
 
 #ifdef PUGS_HAS_COSTO
 
@@ -9,6 +10,8 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
 #include <mesh/MeshFaceBoundary.hpp>
 #include <mesh/MeshFlatFaceBoundary.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
@@ -29,370 +32,149 @@
 #ifdef PUGS_HAS_MPI
 #include <mpi.h>
 #endif   // PUGS_HAS_MPI
+class MeshVariant;
 
-template <size_t Dimension>
 class Serializer
 {
- private:
-  using MeshType     = Mesh<Connectivity<Dimension>>;
-  using Rd           = TinyVector<Dimension>;
-  using MeshDataType = MeshData<Dimension>;
+private:
+  // std::shared_ptr<const MeshVariant> m_mesh_v;
 
- public:
+  // std::shared_ptr<const MeshVariant> m_mesh;
+  // static constexpr size_t Dimension = MeshType::Dimension;
+  // using Rd           = TinyVector<Dimension>;
+  // using MeshDataType = MeshData<MeshType>;
+
+public:
   enum locator : int
-  {
-    at_node,
-    at_face,
-    at_cell_center,
-  };
-  Serializer()  = default;
-  ~Serializer() = default;
-  void mesh(std::shared_ptr<const IMesh> i_mesh, const int64_t dest, const int64_t tag);
+    {
+      at_node,
+      at_face,
+      at_cell_center,
+    };
+  Serializer() = default;
+  // Serializer(const std::shared_ptr<const MeshVariant>& mesh_v) : m_mesh_v{mesh_v} {}
+
+  // ~Serializer() = dafaults;
+  void mesh(const std::shared_ptr<const MeshVariant> mesh_v,
+	    const int64_t dest,
+	    const int64_t tag);
+
   void
   BC(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-     std::shared_ptr<const IMesh> i_mesh,
+     const std::shared_ptr<const MeshVariant> mesh_v,
      const int location,
      const int64_t dest,
      const int64_t tag)
   {
-    this->BC(*boundary, i_mesh, location, dest, tag);
+    this->BC(*boundary, mesh_v, location, dest, tag);
   }
 
   void BC(const IBoundaryDescriptor& boundary,
-          std::shared_ptr<const IMesh> i_mesh,
+	  const std::shared_ptr<const MeshVariant> mesh_v,
           const int location,
           const int64_t dest,
           const int64_t tag);
 
-  void BC_field(std::shared_ptr<const IMesh> i_mesh,
+
+  void BC_field(const std::shared_ptr<const MeshVariant> mesh_v,
                 const std::shared_ptr<const IBoundaryDescriptor>& boundary,
                 const std::shared_ptr<const ItemValueVariant>& field,
                 const int64_t dest,
                 const int64_t tag);
 
-  void interfaces(std::shared_ptr<const IMesh> i_mesh,
+  void interfaces(const std::shared_ptr<const MeshVariant> mesh_v,
                   const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
                   const int64_t dest,
                   const int64_t tag);
 
-  void field(const std::shared_ptr<const DiscreteFunctionVariant>& field, const int64_t dest, const int64_t tag);
+  // void field(const std::shared_ptr<const DiscreteFunctionVariant>& field, const int64_t dest, const int64_t tag);
 
-  std::shared_ptr<const IMesh> change_mesh_position(std::shared_ptr<const IMesh> i_mesh,
-                                                    const int64_t src,
-                                                    const int64_t tag);
+  // std::shared_ptr<const MeshVariant> change_mesh_position(std::shared_ptr<const MeshVariant> mesh_v,
+  // 							  const int64_t src,
+  // 							  const int64_t tag);
 };
 
-// MESH node position
-template <size_t Dimension>
-void
-Serializer<Dimension>::mesh(std::shared_ptr<const IMesh> i_mesh, const int64_t dest, const int64_t tag)
-{
-  auto costo = parallel::Messenger::getInstance().myCoupling;
-  std::vector<int> shape;
-  std::vector<double> pts;
-
-  const std::shared_ptr mesh    = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-  const NodeValue<const Rd>& xr = mesh->xr();
-  Array<TinyVector<Dimension>> positions(mesh->numberOfNodes());
-  parallel_for(
-    mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-      for (unsigned short i = 0; i < Dimension; ++i) {
-        positions[r][i] = xr[r][i];
-      }
-    });
-
-  /*TODO: serializer directement position pour eviter une copie*/
-  pts.resize(mesh->numberOfNodes() * Dimension);
-  for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
-    for (unsigned short j = 0; j < Dimension; ++j) {
-      pts[Dimension * r + j] = positions[r][j];
-    }
-  }
-
-  shape.resize(3);
-  shape[0] = mesh->numberOfNodes();
-  shape[1] = Dimension;
-  shape[2] = 0;
-  costo->sendData(shape, pts, dest, tag);
-  /* const int src  = 1; */
-  /* std::vector<int> recv; */
-  /* parallel::Messenger::getInstance().myCoupling->recvData(recv, src, tag + 1); */
-}
-
-// BC nodes position
-template <size_t Dimension>
-void
-Serializer<Dimension>::BC(const IBoundaryDescriptor& boundary,
-                          std::shared_ptr<const IMesh> i_mesh,
-                          const int location,
-                          const int64_t dest,
-                          const int64_t tag)
-{
-  auto costo = parallel::Messenger::getInstance().myCoupling;
-  std::vector<int> shape;
-  std::vector<double> pts;
-  std::cout << "\033[01;31m"
-            << "L" << location << " ff " << at_face << "\033[00;00m" << std::endl;
-
-  if (location == at_face) {
-    std::cout << "\033[01;31m"
-              << "ICI AT FACE"
-              << "\033[00;00m" << std::endl;
-    const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-    /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
-    const auto xl = mesh_data.xl();
-
-    Array<TinyVector<Dimension>> fpositions(mesh->numberOfFaces());
-    parallel_for(
-      mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
-        for (unsigned short i = 0; i < Dimension; ++i) {
-          fpositions[l][i] = xl[l][i];
-        }
-        /* for (unsigned short i = Dimension; i < 3; ++i) { */
-        /*   positions[r][i] = 0; */
-        /* } */
-      });
-
-    /*TODO: serializer directement fposition pour eviter une copie*/
-    MeshFaceBoundary<Dimension> mesh_face_boundary = getMeshFaceBoundary(*mesh, boundary);
-    /* mesh_face_boundary.faceList() */
-    const auto face_list = mesh_face_boundary.faceList();
-    pts.resize(face_list.size() * Dimension);
-
-    for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-      const FaceId face_id = face_list[i_face];
-      for (unsigned short j = 0; j < Dimension; ++j) {
-        pts[Dimension * i_face + j] = fpositions[face_id][j];
-      }
-      /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-    }
-    shape.resize(3);
-    shape[0] = face_list.size();
-    shape[1] = Dimension;
-    shape[2] = 0;
-    costo->sendData(shape, pts, dest, tag);
-  } else if (location == at_node) {
-    std::cout << "\033[01;31m"
-              << "ICI AT node"
-              << "\033[00;00m" << std::endl;
-    const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    /* const ItemValue<const Rd>& xl = mesh_data.xl(); */
-    const auto xr = mesh->xr();
-
-    Array<TinyVector<Dimension>> positions(mesh->numberOfNodes());
-    parallel_for(
-      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        for (unsigned short i = 0; i < Dimension; ++i) {
-          positions[r][i] = xr[r][i];
-        }
-        /* for (unsigned short i = Dimension; i < 3; ++i) { */
-        /*   positions[r][i] = 0; */
-        /* } */
-      });
-
-    /*TODO: serializer directement position pour eviter une copie*/
-    MeshNodeBoundary<Dimension> mesh_node_boundary = getMeshNodeBoundary(*mesh, boundary);
-    /* mesh_face_boundary.faceList() */
-    const auto node_list = mesh_node_boundary.nodeList();
-    pts.resize(node_list.size() * Dimension);
-
-    for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-      const NodeId node_id = node_list[i_node];
-      for (unsigned short j = 0; j < Dimension; ++j) {
-        pts[Dimension * i_node + j] = positions[node_id][j];
-      }
-      /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-    }
-    shape.resize(3);
-    shape[0] = node_list.size();
-    shape[1] = Dimension;
-    shape[2] = 0;
-    costo->sendData(shape, pts, dest, tag);
-  } else {
-    throw UnexpectedError(" location must be 0 or 1");
-  }
-
-  /* const int src  = 1; */
-  /* std::vector<int> recv; */
-  /* parallel::Messenger::getInstance().myCoupling->recvData(recv, src, tag + 1); */
-}
-
-template <size_t Dimension>
-void
-Serializer<Dimension>::interfaces(std::shared_ptr<const IMesh> i_mesh,
-                                  const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
-                                  const int64_t dest,
-                                  const int64_t tag)
-{
-  auto costo = parallel::Messenger::getInstance().myCoupling;
-  /* const int64_t dest = costo->myGlobalSize() - 1; */
-  std::vector<int> shape;
-  std::vector<double> data;
-
-  const std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-  for (const auto& i_interface : i_interface_list) {
-    auto node_interface = getMeshNodeInterface(*p_mesh, *i_interface);
-    auto node_list      = node_interface.nodeList();
-
-    data.resize(node_list.size() * Dimension);
-
-    const auto xr = p_mesh->xr();
-
-    Array<TinyVector<Dimension>> positions(p_mesh->numberOfNodes());
-    parallel_for(
-      p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        for (unsigned short i = 0; i < Dimension; ++i) {
-          positions[r][i] = xr[r][i];
-        }
-        /* for (unsigned short i = Dimension; i < 3; ++i) { */
-        /*   positions[r][i] = 0; */
-        /* } */
-      });
-
-    /*TODO: serializer directement positions pour eviter une copie */
-    for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-      const NodeId node_id = node_list[i_node];
-      for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
-        data[i_node * Dimension + i_dim] = positions[node_id][i_dim];
-      }
-    }
-
-    shape.resize(3);
-    shape[0] = node_list.size();
-    shape[1] = Dimension;
-    shape[2] = 0;
-    costo->sendData(shape, data, dest, tag);
-  }
-}
-
-// Field at node
-template <size_t Dimension>
-void
-Serializer<Dimension>::BC_field(std::shared_ptr<const IMesh> i_mesh,
-                                const std::shared_ptr<const IBoundaryDescriptor>& boundary,
-                                const std::shared_ptr<const ItemValueVariant>& field,
-                                const int64_t dest,
-                                const int64_t tag)
-{
-  auto costo = parallel::Messenger::getInstance().myCoupling;
-  /* const int64_t dest = costo->myGlobalSize() - 1; */
-  std::vector<int> shape;
-  std::vector<double> data;
-
-  /* std::cout << "\033[01;31m" */
-  /*           << "ival" << ival.type() << "\033[00;00m" << std::endl; */
-  /* std::cout << "\033[01;31m" << ival << "\033[00;00m" << std::endl; */
-  /* const auto test = field->mesh(); */
-  /* const std::shared_ptr i_mesh           = getCommonMesh({field, field}); */
-
-  const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-  if (not i_mesh) {
-    throw NormalError("primal discrete functions are not defined on the same mesh");
-  }
-  const NodeValue<const TinyVector<Dimension>> nField = field->get<NodeValue<const TinyVector<Dimension>>>();
-  /* assert(MeshType(i_mesh) == field->mesh()); */
-
-  MeshNodeBoundary<Dimension> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
-
-  const auto node_list = mesh_node_boundary.nodeList();
-  data.resize(node_list.size() * Dimension);
-
-  for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
-    const NodeId node_id = node_list[i_node];
-    for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
-      data[i_node * Dimension + i_dim] = nField[node_id][i_dim];
-    }
-  }
-
-  shape.resize(3);
-  shape[0] = node_list.size();
-  shape[1] = Dimension;
-  shape[2] = 0;
-  costo->sendData(shape, data, dest, tag);
-}
-
-// Field (faces)
-template <size_t Dimension>
-void
-Serializer<Dimension>::field(const std::shared_ptr<const DiscreteFunctionVariant>& field,
-                             const int64_t dest,
-                             const int64_t tag)
-{
-  const std::shared_ptr i_mesh = getCommonMesh({field});
-  if (not i_mesh) {
-    throw NormalError("primal discrete functions are not defined on the same mesh");
-  }
-  /* assert(MeshType(i_mesh) == field->mesh()); */
-  const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-  /* const IBoundaryDescriptor& bd = *boundary; */
-  /* const Mesh<Connectivity<2>> mt = <const MeshType>(i_mesh); */
-
-  /* const auto& n_list = MeshNodeBoundary<2>(mt, bd); */
-
-  const NodeValue<const Rd>& xr = mesh->xr();
-  Array<TinyVector<Dimension>> positions(mesh->numberOfNodes());
-  parallel_for(
-    mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-      for (unsigned short i = 0; i < Dimension; ++i) {
-        positions[r][i] = xr[r][i];
-      }
-      /* for (unsigned short i = Dimension; i < 3; ++i) { */
-      /*   positions[r][i] = 0; */
-      /* } */
-    });
-
-  std::vector<double> pts;
-  pts.resize(mesh->numberOfNodes() * Dimension);
-  for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
-    for (unsigned short j = 0; j < Dimension; ++j) {
-      pts[Dimension * r + j] = positions[r][j];
-    }
-    /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
-  }
-
-  auto costo = parallel::Messenger::getInstance().myCoupling;
-  /* const int64_t dest = costo->myGlobalSize() - 1; */
-  std::vector<int> shape;
-  shape.resize(3);
-  shape[0] = mesh->numberOfNodes();
-  shape[1] = Dimension;
-  shape[2] = 0;
-  costo->sendData(shape, pts, dest, tag);
-
-  std::vector<int> recv;
-  const int src = costo->myGlobalSize() - 1;
-
-  costo->recvData(recv, src, tag + 1);
-}
-
-template <size_t Dimension>
-std::shared_ptr<const IMesh>
-Serializer<Dimension>::change_mesh_position(std::shared_ptr<const IMesh> i_mesh, const int64_t src, const int64_t tag)
-{
-  std::vector<double> recv_xr;
-  parallel::Messenger::getInstance().myCoupling->recvData(recv_xr, src, tag);
 
-  const std::shared_ptr given_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-  NodeValue<Rd> new_xr{given_mesh->connectivity()};
-  parallel_for(
-    given_mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-      for (unsigned short j = 0; j < Dimension; ++j) {
-        /* new_xr[node_id][j] = xr[node_id][j]; */
-        new_xr[node_id][j] = recv_xr[Dimension * node_id + j];
-      }
-    });
 
-  return std::make_shared<MeshType>(given_mesh->shared_connectivity(), new_xr);
-}
 
-using Serializer3D = Serializer<3>;
-using Serializer2D = Serializer<2>;
-using Serializer1D = Serializer<1>;
+// // Field (faces)
+// template <MeshConcept MeshType>void
+// Serializer<MeshType>::field(const std::shared_ptr<const DiscreteFunctionVariant>& field,
+//                              const int64_t dest,
+//                              const int64_t tag)
+// {
+//   const std::shared_ptr i_mesh = getCommonMesh({field});
+//   if (not i_mesh) {
+//     throw NormalError("primal discrete functions are not defined on the same mesh");
+//   }
+//   /* assert(MeshType(i_mesh) == field->mesh()); */
+//   const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+//   /* const IBoundaryDescriptor& bd = *boundary; */
+//   /* const Mesh<Connectivity<2>> mt = <const MeshType>(i_mesh); */
+
+//   /* const auto& n_list = MeshNodeBoundary<2>(mt, bd); */
+
+//   const NodeValue<const Rd>& xr = mesh->xr();
+//   Array<TinyVector<Dimension>> positions(mesh->numberOfNodes());
+//   parallel_for(
+//     mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+//       for (unsigned short i = 0; i < Dimension; ++i) {
+//         positions[r][i] = xr[r][i];
+//       }
+//       /* for (unsigned short i = Dimension; i < 3; ++i) { */
+//       /*   positions[r][i] = 0; */
+//       /* } */
+//     });
+
+//   std::vector<double> pts;
+//   pts.resize(mesh->numberOfNodes() * Dimension);
+//   for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
+//     for (unsigned short j = 0; j < Dimension; ++j) {
+//       pts[Dimension * r + j] = positions[r][j];
+//     }
+//     /* std::cout << "\033[01;31m" << pts[2 * r] << "; " << pts[2 * r + 1] << "\033[00;00m" << std::endl; */
+//   }
+
+//   auto costo = parallel::Messenger::getInstance().myCoupling;
+//   /* const int64_t dest = costo->myGlobalSize() - 1; */
+//   std::vector<int> shape;
+//   shape.resize(3);
+//   shape[0] = mesh->numberOfNodes();
+//   shape[1] = Dimension;
+//   shape[2] = 0;
+//   costo->sendData(shape, pts, dest, tag);
+
+//   std::vector<int> recv;
+//   const int src = costo->myGlobalSize() - 1;
+
+//   costo->recvData(recv, src, tag + 1);
+// }
+
+// template <MeshConcept MeshType>
+// std::shared_ptr<const MeshVariant>
+// Serializer<MeshType>::change_mesh_position(std::shared_ptr<const MeshVariant> mesh_v, const int64_t src, const int64_t tag)
+// {
+//   std::vector<double> recv_xr;
+//   parallel::Messenger::getInstance().myCoupling->recvData(recv_xr, src, tag);
+
+//   const std::shared_ptr given_mesh = std::dynamic_pointer_cast<const MeshType>(mesh_v);
+
+//   NodeValue<Rd> new_xr{given_mesh->connectivity()};
+//   parallel_for(
+//     given_mesh->numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+//       for (unsigned short j = 0; j < Dimension; ++j) {
+//         /* new_xr[node_id][j] = xr[node_id][j]; */
+//         new_xr[node_id][j] = recv_xr[Dimension * node_id + j];
+//       }
+//     });
+
+//   return std::make_shared<MeshVariant>(given_mesh->shared_connectivity(), new_xr);
+// }
+
+// using Serializer3D = Serializer<Mesh<1>>;
+// using Serializer2D = Serializer<Mesh<2>>;
+// using Serializer1D = Serializer<Mesh<3>>;
 
 #endif   // PUGS_HAS_COSTO
+// #endif   // SERIALIZER_HPP
diff --git a/tests/test_LeastSquareSolver.cpp b/tests/test_LeastSquareSolver.cpp
deleted file mode 100644
index 6022bbe151ff2c261dc21b015692b9045edcc9ca..0000000000000000000000000000000000000000
--- a/tests/test_LeastSquareSolver.cpp
+++ /dev/null
@@ -1,77 +0,0 @@
-#include <catch2/catch.hpp>
-
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LocalRectangularMatrix.hpp>
-#include <algebra/Vector.hpp>
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("LeastSquareSolver", "[algebra]")
-{
-  SECTION("Least Squares [under-determined]")
-  {
-    Vector<double> b{3};
-    b[0] = 1;
-    b[1] = 0;
-    b[2] = 0;
-
-    LocalRectangularMatrix<double> A{3, 4};
-    A(0, 0) = 1;
-    A(0, 1) = 1;
-    A(0, 2) = 1;
-    A(0, 3) = 1;
-    A(1, 0) = 1;
-    A(1, 1) = -1;
-    A(1, 2) = 0;
-    A(1, 3) = 0;
-    A(2, 0) = 0;
-    A(2, 1) = 0;
-    A(2, 2) = 1;
-    A(2, 3) = -1;
-
-    Vector<double> x_exact{4};
-
-    x_exact[0] = 0.25;
-    x_exact[1] = 0.25;
-    x_exact[2] = 0.25;
-    x_exact[3] = 0.25;
-
-    Vector<double> x{4};
-    x = 0;
-
-    LeastSquareSolver ls_solver;
-    ls_solver.solveLocalSystem(A, x, b);
-
-    Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
-  }
-
-  SECTION("Least Squares [over-determined]")
-  {
-    LocalRectangularMatrix<double> A{3, 2};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(1, 0) = 1;
-    A(1, 1) = 1;
-    A(2, 0) = 2;
-    A(2, 1) = 1;
-
-    Vector<double> x_exact{2};
-    x_exact[0] = -3;
-    x_exact[1] = 5;
-
-    Vector b{3};
-    b[0] = 6;
-    b[1] = 0;
-    b[2] = 0;
-
-    Vector<double> x{2};
-    x = 0;
-
-    LeastSquareSolver ls_solver;
-    ls_solver.solveLocalSystem(A, x, b);
-
-    Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
-  }
-}