From 244111df1fcfac885800805ac328e8adf50b37d1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 23 Feb 2023 00:19:17 +0100
Subject: [PATCH] Remove algorithm legacy files

---
 src/language/algorithms/CMakeLists.txt        |   8 +-
 .../algorithms/ElasticityDiamondAlgorithm.cpp | 874 ------------------
 .../algorithms/ElasticityDiamondAlgorithm.hpp |  32 -
 .../algorithms/Heat5PointsAlgorithm.cpp       | 207 -----
 .../algorithms/Heat5PointsAlgorithm.hpp       |  21 -
 .../algorithms/HeatDiamondAlgorithm.cpp       | 773 ----------------
 .../algorithms/HeatDiamondAlgorithm.hpp       |  29 -
 .../algorithms/HeatDiamondAlgorithm2.cpp      | 869 -----------------
 .../algorithms/HeatDiamondAlgorithm2.hpp      |  29 -
 src/language/algorithms/ParabolicHeat.cpp     | 568 ------------
 src/language/algorithms/ParabolicHeat.hpp     |  33 -
 .../algorithms/UnsteadyElasticity.cpp         | 613 ------------
 .../algorithms/UnsteadyElasticity.hpp         |  36 -
 src/language/modules/SchemeModule.cpp         | 187 ----
 14 files changed, 1 insertion(+), 4278 deletions(-)
 delete mode 100644 src/language/algorithms/ElasticityDiamondAlgorithm.cpp
 delete mode 100644 src/language/algorithms/ElasticityDiamondAlgorithm.hpp
 delete mode 100644 src/language/algorithms/Heat5PointsAlgorithm.cpp
 delete mode 100644 src/language/algorithms/Heat5PointsAlgorithm.hpp
 delete mode 100644 src/language/algorithms/HeatDiamondAlgorithm.cpp
 delete mode 100644 src/language/algorithms/HeatDiamondAlgorithm.hpp
 delete mode 100644 src/language/algorithms/HeatDiamondAlgorithm2.cpp
 delete mode 100644 src/language/algorithms/HeatDiamondAlgorithm2.hpp
 delete mode 100644 src/language/algorithms/ParabolicHeat.cpp
 delete mode 100644 src/language/algorithms/ParabolicHeat.hpp
 delete mode 100644 src/language/algorithms/UnsteadyElasticity.cpp
 delete mode 100644 src/language/algorithms/UnsteadyElasticity.hpp

diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt
index 50783e8ea..63c237498 100644
--- a/src/language/algorithms/CMakeLists.txt
+++ b/src/language/algorithms/CMakeLists.txt
@@ -1,13 +1,7 @@
 # ------------------- Source files --------------------
 
 add_library(PugsLanguageAlgorithms
-  ElasticityDiamondAlgorithm.cpp
-  UnsteadyElasticity.cpp
-  Heat5PointsAlgorithm.cpp
-  HeatDiamondAlgorithm.cpp
-  HeatDiamondAlgorithm2.cpp
-  ParabolicHeat.cpp
-#  INTERFACE # THIS SHOULD BE REMOVED IF FILES ARE FINALY PROVIDED
+  INTERFACE # THIS SHOULD BE REMOVED IF FILES ARE FINALY PROVIDED
   )
 
 
diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
deleted file mode 100644
index b14557ebd..000000000
--- a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
+++ /dev/null
@@ -1,874 +0,0 @@
-#include <language/algorithms/ElasticityDiamondAlgorithm.hpp>
-
-#include <algebra/CRSMatrix.hpp>
-#include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LinearSolver.hpp>
-#include <algebra/SmallMatrix.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/PrimalToDiamondDualConnectivityDataMapper.hpp>
-#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
-#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-
-template <size_t Dimension>
-ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
-  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,
-  const FunctionSymbolId& U_id)
-{
-  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) {
-        MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(*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 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("Dirichlet conditions are not supported in 1d");
-        }
-      } else if (dirichlet_bc_descriptor.name() == "normal_strain") {
-        if constexpr (Dimension > 1) {
-          MeshFaceBoundary 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 NodeValue<const TinyVector<Dimension>>& xr = 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                  = mesh->connectivity().nodeToCellMatrix();
-    const auto& node_to_face_matrix                  = mesh->connectivity().nodeToFaceMatrix();
-    CellValuePerNode<double> w_rj{mesh->connectivity()};
-    FaceValuePerNode<double> w_rl{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 < 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 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 (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 (primal_face_is_on_boundary[face_id]) {
-              if (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 (primal_face_is_on_boundary[face_id]) {
-            w_rl(i_node, i_face) = x[cpt_face++];
-          }
-        }
-      }
-    }
-
-    {
-      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>> Uj = InterpolateItemValue<TinyVector<Dimension>(
-        TinyVector<Dimension>)>::template interpolate<ItemType::cell>(U_id, 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;
-      }();
-
-      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 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) += 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) -= 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) -= 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) -=
-                        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) -=
-                            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;
-          }
-        }
-      }
-
-      CRSMatrix A{S.getCRSMatrix()};
-      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] += 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);
-      }
-
-      Vector<double> U{number_of_dof * Dimension};
-      U = zero;
-      CellValue<TinyVector<Dimension>> Speed{mesh->connectivity()};
-      FaceValue<TinyVector<Dimension>> Ul = InterpolateItemValue<TinyVector<Dimension>(
-        TinyVector<Dimension>)>::template interpolate<ItemType::face>(U_id, mesh_data.xl());
-      FaceValue<TinyVector<Dimension>> Speed_face{mesh->connectivity()};
-
-      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) {
-          Speed[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]) {
-            Speed_face[face_id][i] = U[(face_dof_number[face_id] * Dimension) + i];
-          } else {
-            Speed_face[face_id][i] = Ul[face_id][i];
-          }
-        }
-      }
-      Vector<double> Uexacte{mesh->numberOfCells() * Dimension};
-      for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
-        for (size_t l = 0; l < Dimension; ++l) {
-          Uexacte[(cell_dof_number[j] * Dimension) + l] = Uj[j][l];
-        }
-      }
-
-      Vector<double> error{mesh->numberOfCells() * Dimension};
-      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-        for (size_t i = 0; i < Dimension; ++i) {
-          error[(cell_id * Dimension) + i] = (Speed[cell_id][i] - Uj[cell_id][i]) * sqrt(primal_Vj[cell_id]);
-        }
-      }
-      Vector<double> error_face{mesh->numberOfFaces() * Dimension};
-      parallel_for(
-        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-          if (primal_face_is_on_boundary[face_id]) {
-            for (size_t i = 0; i < Dimension; ++i) {
-              error_face[face_id * Dimension + i] = (Speed_face[face_id][i] - Ul[face_id][i]) * sqrt(mes_l[face_id]);
-            }
-          } else {
-            error_face[face_id] = 0;
-          }
-        });
-
-      std::cout << "||Error||_2 (cell)= " << std::sqrt(dot(error, error)) << "\n";
-      std::cout << "||Error||_2 (face)= " << std::sqrt(dot(error_face, error_face)) << "\n";
-      std::cout << "||Error||_2 (total)= " << std::sqrt(dot(error, error)) + std::sqrt(dot(error_face, error_face))
-                << "\n";
-
-      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) * Speed[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) * Speed_face[face_id];
-            }
-          }
-          for (size_t i = 0; i < Dimension; ++i) {
-            ur3d[node_id][i] = x[i];
-          }
-        });
-    }
-  } else {
-    throw NotImplementedError("not done in 1d");
-  }
-}
-
-template <size_t Dimension>
-class ElasticityDiamondScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class ElasticityDiamondScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class ElasticityDiamondScheme<Dimension>::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;
-};
-
-template ElasticityDiamondScheme<1>::ElasticityDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template ElasticityDiamondScheme<2>::ElasticityDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template ElasticityDiamondScheme<3>::ElasticityDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.hpp b/src/language/algorithms/ElasticityDiamondAlgorithm.hpp
deleted file mode 100644
index 24988916e..000000000
--- a/src/language/algorithms/ElasticityDiamondAlgorithm.hpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef ELASTICITY_DIAMOND_ALGORITHM_HPP
-#define ELASTICITY_DIAMOND_ALGORITHM_HPP
-
-#include <algebra/TinyVector.hpp>
-#include <language/utils/FunctionSymbolId.hpp>
-#include <memory>
-#include <mesh/IMesh.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <variant>
-#include <vector>
-
-template <size_t Dimension>
-class ElasticityDiamondScheme
-{
- private:
-  class DirichletBoundaryCondition;
-  class NormalStrainBoundaryCondition;
-  class SymmetryBoundaryCondition;
-
- public:
-  ElasticityDiamondScheme(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,
-                          const FunctionSymbolId& U_id);
-};
-
-#endif   // ELASTICITY_DIAMOND_ALGORITHM2_HPP
diff --git a/src/language/algorithms/Heat5PointsAlgorithm.cpp b/src/language/algorithms/Heat5PointsAlgorithm.cpp
deleted file mode 100644
index ffd1e8102..000000000
--- a/src/language/algorithms/Heat5PointsAlgorithm.cpp
+++ /dev/null
@@ -1,207 +0,0 @@
-#include <language/algorithms/Heat5PointsAlgorithm.hpp>
-
-#include <algebra/CRSMatrix.hpp>
-#include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LinearSolver.hpp>
-#include <algebra/LinearSolverOptions.hpp>
-#include <algebra/SmallMatrix.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/PrimalToDiamondDualConnectivityDataMapper.hpp>
-
-template <size_t Dimension>
-Heat5PointsAlgorithm<Dimension>::Heat5PointsAlgorithm(
-  std::shared_ptr<const IMesh> i_mesh,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-  const FunctionSymbolId& T_id,
-  const FunctionSymbolId& kappa_id,
-  const FunctionSymbolId& f_id)
-{
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshData<Dimension>;
-
-  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-  std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
-
-  if constexpr (Dimension == 2) {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-    CellValue<double> Tj =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
-
-    NodeValue<double> Tr(mesh->connectivity());
-    const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
-    const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-    const auto& node_to_cell_matrix                  = mesh->connectivity().nodeToCellMatrix();
-    CellValuePerNode<double> w_rj{mesh->connectivity()};
-
-    for (NodeId i_node = 0; i_node < 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];
-
-      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);
-
-      Tr[i_node] = 0;
-      for (size_t j = 0; j < node_to_cell.size(); j++) {
-        Tr[i_node] += x[j] * Tj[node_to_cell[j]];
-        w_rj(i_node, j) = x[j];
-      }
-    }
-
-    {
-      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());
-
-      NodeValue<double> Trd{diamond_mesh->connectivity()};
-
-      mapper->toDualNode(Tr, Tj, Trd);
-
-      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 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_mes_l_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;
-      }();
-
-      const Array<int> non_zeros{mesh->numberOfCells()};
-      non_zeros.fill(Dimension);
-      CRSMatrixDescriptor<double> S(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
-
-      const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
-      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 = 0.5 * 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_id1, cell_id2) -= beta_l;
-            } else {
-              S(cell_id1, cell_id2) += beta_l;
-            }
-          }
-        }
-      }
-      CellValue<double> fj =
-        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
-
-      const CellValue<const double> primal_Vj = mesh_data.Vj();
-      CRSMatrix A{S.getCRSMatrix()};
-      Vector<double> b{mesh->numberOfCells()};
-      parallel_for(
-        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; });
-
-      Vector<double> T{mesh->numberOfCells()};
-      T = zero;
-
-      LinearSolver solver;
-      solver.solveLocalSystem(A, T, b);
-
-      CellValue<double> Temperature{mesh->connectivity()};
-
-      parallel_for(
-        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_id]; });
-
-      Vector<double> error{mesh->numberOfCells()};
-      parallel_for(
-        mesh->numberOfCells(),
-        PUGS_LAMBDA(CellId cell_id) { error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * primal_Vj[cell_id]; });
-
-      std::cout << "||Error||_2 = " << std::sqrt(dot(error, error)) << "\n";
-    }
-  } else {
-    throw NotImplementedError("not done in this dimension");
-  }
-}
-
-template Heat5PointsAlgorithm<1>::Heat5PointsAlgorithm(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template Heat5PointsAlgorithm<2>::Heat5PointsAlgorithm(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template Heat5PointsAlgorithm<3>::Heat5PointsAlgorithm(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
diff --git a/src/language/algorithms/Heat5PointsAlgorithm.hpp b/src/language/algorithms/Heat5PointsAlgorithm.hpp
deleted file mode 100644
index a061972d6..000000000
--- a/src/language/algorithms/Heat5PointsAlgorithm.hpp
+++ /dev/null
@@ -1,21 +0,0 @@
-#ifndef HEAT_5POINTS_ALGORITHM_HPP
-#define HEAT_5POINTS_ALGORITHM_HPP
-
-#include <language/utils/FunctionSymbolId.hpp>
-#include <mesh/IMesh.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-
-#include <memory>
-#include <vector>
-
-template <size_t Dimension>
-struct Heat5PointsAlgorithm
-{
-  Heat5PointsAlgorithm(std::shared_ptr<const IMesh> i_mesh,
-                       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-                       const FunctionSymbolId& T_id,
-                       const FunctionSymbolId& kappa_id,
-                       const FunctionSymbolId& f_id);
-};
-
-#endif   // HEAT_5POINTS_ALGORITHM_HPP
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
deleted file mode 100644
index d21c25062..000000000
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ /dev/null
@@ -1,773 +0,0 @@
-#include <language/algorithms/HeatDiamondAlgorithm.hpp>
-
-#include <algebra/CRSMatrix.hpp>
-#include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LinearSolver.hpp>
-#include <algebra/SmallMatrix.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/FourierBoundaryConditionDescriptor.hpp>
-#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-#include <utils/Timer.hpp>
-
-template <size_t Dimension>
-HeatDiamondScheme<Dimension>::HeatDiamondScheme(
-  std::shared_ptr<const IMesh> i_mesh,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-  const FunctionSymbolId& T_id,
-  const FunctionSymbolId& kappa_id,
-  const FunctionSymbolId& f_id)
-{
-  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::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
-  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});
-      } else {
-        throw NotImplementedError("not implemented 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("not implemented 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);
-
-    CellValue<double> Tj =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
-    FaceValue<double> Tl =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
-    NodeValue<double> Tr(mesh->connectivity());
-    const NodeValue<const TinyVector<Dimension>>& xr = 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                  = mesh->connectivity().nodeToCellMatrix();
-    const auto& node_to_face_matrix                  = mesh->connectivity().nodeToFaceMatrix();
-    CellValuePerNode<double> w_rj{mesh->connectivity()};
-    FaceValuePerNode<double> w_rl{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 < 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 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);
-
-        Tr[i_node] = 0;
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          Tr[i_node] += x[j] * Tj[node_to_cell[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 (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 (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);
-
-        Tr[i_node] = 0;
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          Tr[i_node] += x[j] * Tj[node_to_cell[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 (primal_face_is_on_boundary[face_id]) {
-            w_rl(i_node, i_face) = x[cpt_face++];
-            Tr[i_node] += w_rl(i_node, i_face) * Tl[face_id];
-          }
-        }
-      }
-    }
-
-    {
-      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());
-
-      NodeValue<double> Trd{diamond_mesh->connectivity()};
-
-      mapper->toDualNode(Tr, Tj, Trd);
-
-      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;
-      }();
-      Timer my_timer;
-      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;
-      }();
-
-      const Array<int> non_zeros{number_of_dof};
-      non_zeros.fill(Dimension * 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             = 1. / Dimension * alpha_l[face_id] * mes_l[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]) += beta_l;
-              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]) -= 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]) -= 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]) -= 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());
-
-      const CellValue<const double> primal_Vj = mesh_data.Vj();
-      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];   // sign
-                }
-              }
-            },
-            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];   // sign
-              }
-            }
-          },
-          boundary_condition);
-      }
-
-      Vector<double> T{number_of_dof};
-      T = zero;
-
-      LinearSolver solver;
-      solver.solveLocalSystem(A, T, b);
-
-      CellValue<double> Temperature{mesh->connectivity()};
-      FaceValue<double> Temperature_face{mesh->connectivity()};
-      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]];
-          } else {
-            Temperature_face[face_id] = Tl[face_id];
-          }
-        });
-      Vector<double> error{mesh->numberOfCells()};
-      CellValue<double> cell_error{mesh->connectivity()};
-      Vector<double> face_error{mesh->numberOfFaces()};
-      double error_max = 0.;
-      size_t cell_max  = 0;
-      parallel_for(
-        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-          error[cell_id]      = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
-          cell_error[cell_id] = (Temperature[cell_id] - Tj[cell_id]);
-        });
-      parallel_for(
-        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-          if (primal_face_is_on_boundary[face_id]) {
-            face_error[face_id] = (Temperature_face[face_id] - Tl[face_id]) * sqrt(mes_l[face_id]);
-          } else {
-            face_error[face_id] = 0;
-          }
-        });
-      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); cell_id++) {
-        if (error_max < std::abs(cell_error[cell_id])) {
-          error_max = std::abs(cell_error[cell_id]);
-          cell_max  = cell_id;
-        }
-      }
-
-      std::cout << " ||Error||_max (cell)= " << error_max << " on cell " << cell_max << "\n";
-      std::cout << "||Error||_2 (cell)= " << std::sqrt(dot(error, error)) << "\n";
-      std::cout << "||Error||_2 (face)= " << std::sqrt(dot(face_error, face_error)) << "\n";
-      std::cout << "||Error||_2 (total)= " << std::sqrt(dot(error, error)) + std::sqrt(dot(face_error, face_error))
-                << "\n";
-    }
-  } else {
-    throw NotImplementedError("not implemented in 1d");
-  }
-}
-
-template <size_t Dimension>
-class HeatDiamondScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class HeatDiamondScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class HeatDiamondScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class HeatDiamondScheme<Dimension>::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;
-};
-
-template HeatDiamondScheme<1>::HeatDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template HeatDiamondScheme<2>::HeatDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template HeatDiamondScheme<3>::HeatDiamondScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.hpp b/src/language/algorithms/HeatDiamondAlgorithm.hpp
deleted file mode 100644
index 9e29d91bc..000000000
--- a/src/language/algorithms/HeatDiamondAlgorithm.hpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#ifndef HEAT_DIAMOND_ALGORITHM_HPP
-#define HEAT_DIAMOND_ALGORITHM_HPP
-
-#include <language/utils/FunctionSymbolId.hpp>
-#include <mesh/IMesh.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-
-#include <memory>
-#include <variant>
-#include <vector>
-
-template <size_t Dimension>
-class HeatDiamondScheme
-{
- private:
-  class DirichletBoundaryCondition;
-  class FourierBoundaryCondition;
-  class NeumannBoundaryCondition;
-  class SymmetryBoundaryCondition;
-
- public:
-  HeatDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
-                    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-                    const FunctionSymbolId& T_id,
-                    const FunctionSymbolId& kappa_id,
-                    const FunctionSymbolId& f_id);
-};
-
-#endif   // HEAT_DIAMOND_ALGORITHM_HPP
diff --git a/src/language/algorithms/HeatDiamondAlgorithm2.cpp b/src/language/algorithms/HeatDiamondAlgorithm2.cpp
deleted file mode 100644
index c36e79678..000000000
--- a/src/language/algorithms/HeatDiamondAlgorithm2.cpp
+++ /dev/null
@@ -1,869 +0,0 @@
-#include <language/algorithms/HeatDiamondAlgorithm2.hpp>
-
-#include <algebra/CRSMatrix.hpp>
-#include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LinearSolver.hpp>
-#include <algebra/SmallMatrix.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/FourierBoundaryConditionDescriptor.hpp>
-#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-#include <utils/Timer.hpp>
-
-template <size_t Dimension>
-HeatDiamondScheme2<Dimension>::HeatDiamondScheme2(
-  std::shared_ptr<const IMesh> i_mesh,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-  const FunctionSymbolId& T_id,
-  const FunctionSymbolId& kappa_id,
-  const FunctionSymbolId& f_id)
-{
-  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::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<double> dirichlet_value{mesh->connectivity()};
-  dirichlet_value.fill(std::numeric_limits<double>::signaling_NaN());
-
-  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) {
-        MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-        const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-
-        MeshNodeBoundary<Dimension> mesh_node_boundary =
-          getMeshNodeBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-        Array<const double> node_value_list =
-          InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
-                                                                                                    mesh_node_boundary
-                                                                                                      .nodeList());
-
-        MeshFaceBoundary<Dimension> mesh_face_boundary =
-          getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
-
-        Array<const double> face_value_list =
-          InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
-                                                                                                    mesh_data.xl(),
-                                                                                                    mesh_face_boundary
-                                                                                                      .faceList());
-
-        for (size_t i_node = 0; i_node < mesh_node_boundary.nodeList().size(); ++i_node) {
-          NodeId node_id = mesh_node_boundary.nodeList()[i_node];
-
-          if (not is_dirichlet[node_id]) {
-            is_dirichlet[node_id]    = true;
-            dirichlet_value[node_id] = node_value_list[i_node];
-          }
-        }
-
-        boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), face_value_list,
-                                                                     mesh_node_boundary.nodeList(), node_value_list});
-
-      } else {
-        throw NotImplementedError("not implemented 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) {
-        MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-        const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
-
-        MeshFaceBoundary<Dimension> mesh_face_boundary =
-          getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
-
-        Array<const double> face_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(), face_value_list});
-      } else {
-        throw NotImplementedError("not implemented 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>)) {
-              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);
-
-    CellValue<double> Tj =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
-    FaceValue<double> Tl =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
-    NodeValue<double> Tr(mesh->connectivity());
-    const NodeValue<const TinyVector<Dimension>>& xr = 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                  = mesh->connectivity().nodeToCellMatrix();
-    const auto& node_to_face_matrix                  = mesh->connectivity().nodeToFaceMatrix();
-    CellValuePerNode<double> w_rj{mesh->connectivity()};
-    FaceValuePerNode<double> w_rl{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 < 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 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);
-
-        Tr[i_node] = 0;
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          Tr[i_node] += x[j] * Tj[node_to_cell[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 (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 (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);
-
-        Tr[i_node] = 0;
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          Tr[i_node] += x[j] * Tj[node_to_cell[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 (primal_face_is_on_boundary[face_id]) {
-            w_rl(i_node, i_face) = x[cpt_face++];
-            Tr[i_node] += w_rl(i_node, i_face) * Tl[face_id];
-          }
-        }
-      }
-    }
-
-    {
-      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());
-
-      NodeValue<double> Trd{diamond_mesh->connectivity()};
-
-      mapper->toDualNode(Tr, Tj, Trd);
-
-      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;
-      }();
-      Timer my_timer;
-      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;
-      }();
-
-      const Array<int> non_zeros{number_of_dof};
-      non_zeros.fill(Dimension * 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             = 1. / Dimension * alpha_l[face_id] * mes_l[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]) += beta_l;
-              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]) -= 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]) -= 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];
-                    // ELIMINATION METTRE AU SECOND MEMBRE
-                    if (primal_face_is_neumann[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) * a_ir;
-                      }
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
-      }
-
-      CellValue<double> fj =
-        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
-
-      const CellValue<const double> primal_Vj = mesh_data.Vj();
-      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 a Dirichlet face F contribute to the second member J via the node R thanks to A^{JR} A^{RF}
-      // for cell and A^{NR} A^{RF} for Neumann faces
-      {
-        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& value_list = bc.valueList();
-                const auto& face_list  = bc.faceList();
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  FaceId face_id = face_list[i_face];
-
-                  // loop on Nodes of Dirichlet faces
-                  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];
-                    // loop on faces of node
-                    const double w_face_node = [&] {
-                      for (size_t i_node_face = 0; i_node_face < node_to_face_matrix[node_id].size(); ++i_node_face) {
-                        FaceId l_id = node_to_face_matrix[node_id][i_node_face];
-                        if (l_id == face_id) {
-                          return w_rl(node_id, i_node_face);
-                        }
-                      }
-                      throw UnexpectedError("unable to get face node weight");
-                    }();
-
-                    for (size_t i_node_face = 0; i_node_face < node_to_face_matrix[node_id].size(); ++i_node_face) {
-                      FaceId l_id         = node_to_face_matrix[node_id][i_node_face];
-                      CellId dual_cell_id = face_dual_cell_id[l_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];
-                        // find the dual node of interest
-                        if (dual_node_primal_node_id[dual_node_id] == node_id) {
-                          // get Cjr
-                          const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-                          const double alpha_face_id      = mes_l[l_id] * alpha_l[l_id];
-                          for (size_t cell_id = 0; cell_id < face_to_cell_matrix[l_id].size(); ++cell_id) {
-                            CellId i_id                            = face_to_cell_matrix[l_id][cell_id];
-                            const bool is_face_reversed_for_cell_i = (dot(dualClj[l_id], xl[l_id] - xj[i_id]) < 0);
-
-                            const TinyVector<Dimension> nil = [&] {
-                              if (is_face_reversed_for_cell_i) {
-                                return -nlj[l_id];
-                              } else {
-                                return nlj[l_id];
-                              }
-                            }();
-
-                            const double a_ir = alpha_face_id * dot(nil, Clr);
-                            // loop on faces of cells
-                            b[cell_dof_number[i_id]] += w_face_node * a_ir * value_list[i_face];
-
-                            // If l_id is Neumann do...
-                            if (primal_face_is_neumann[l_id]) {
-                              b[face_dof_number[l_id]] -= w_face_node * a_ir * 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];
-                Assert(face_to_cell_matrix[face_id].size() == 1);
-                b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face];   // sign
-              }
-            }
-          },
-          boundary_condition);
-      }
-
-      Vector<double> T{number_of_dof};
-      T = zero;
-
-      LinearSolver solver;
-      solver.solveLocalSystem(A, T, b);
-
-      CellValue<double> Temperature{mesh->connectivity()};
-      FaceValue<double> Temperature_face{mesh->connectivity()};
-      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]];
-          } else {
-            Temperature_face[face_id] = Tl[face_id];
-          }
-        });
-      Vector<double> error{mesh->numberOfCells()};
-      CellValue<double> cell_error{mesh->connectivity()};
-      Vector<double> face_error{mesh->numberOfFaces()};
-      double error_max = 0.;
-      size_t cell_max  = 0;
-      parallel_for(
-        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-          error[cell_id]      = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
-          cell_error[cell_id] = (Temperature[cell_id] - Tj[cell_id]);
-        });
-      parallel_for(
-        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-          // ELIMINATION ENLEVER DIRICHLET
-          if (primal_face_is_neumann[face_id]) {
-            face_error[face_id] = (Temperature_face[face_id] - Tl[face_id]) * sqrt(mes_l[face_id]);
-          } else {
-            face_error[face_id] = 0;
-          }
-        });
-      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); cell_id++) {
-        if (error_max < std::abs(cell_error[cell_id])) {
-          error_max = std::abs(cell_error[cell_id]);
-          cell_max  = cell_id;
-        }
-      }
-
-      std::cout << " ||Error||_max (cell)= " << error_max << " on cell " << cell_max << "\n";
-      std::cout << "||Error||_2 (cell)= " << std::sqrt(dot(error, error)) << "\n";
-      std::cout << "||Error||_2 (face)= " << std::sqrt(dot(face_error, face_error)) << "\n";
-      std::cout << "||Error||_2 (total)= " << std::sqrt(dot(error, error)) + std::sqrt(dot(face_error, face_error))
-                << "\n";
-    }
-  } else {
-    throw NotImplementedError("not done in 1d");
-  }
-}
-
-template <size_t Dimension>
-class HeatDiamondScheme2<Dimension>::DirichletBoundaryCondition
-{
- private:
-  const Array<const double> m_value_list;
-  const Array<const FaceId> m_face_list;
-
-  const Array<const double> m_node_value_list;
-  const Array<const NodeId> m_node_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 NodeId>&
-  nodeList() const
-  {
-    return m_node_list;
-  }
-
-  const Array<const double>&
-  nodeValueList() const
-  {
-    return m_node_value_list;
-  }
-
-  DirichletBoundaryCondition(const Array<const FaceId>& face_list,
-                             const Array<const double>& value_list,
-                             const Array<const NodeId>& node_list,
-                             const Array<const double>& node_value_list)
-    : m_value_list{value_list}, m_face_list{face_list}, m_node_value_list(node_value_list), m_node_list(node_list)
-  {
-    Assert(m_value_list.size() == m_face_list.size());
-    Assert(m_node_value_list.size() == m_node_list.size());
-  }
-
-  ~DirichletBoundaryCondition() = default;
-};
-
-template <size_t Dimension>
-class HeatDiamondScheme2<Dimension>::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;
-};
-
-template <size_t Dimension>
-class HeatDiamondScheme2<Dimension>::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;
-};
-
-template <size_t Dimension>
-class HeatDiamondScheme2<Dimension>::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;
-};
-
-template HeatDiamondScheme2<1>::HeatDiamondScheme2(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template HeatDiamondScheme2<2>::HeatDiamondScheme2(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template HeatDiamondScheme2<3>::HeatDiamondScheme2(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
diff --git a/src/language/algorithms/HeatDiamondAlgorithm2.hpp b/src/language/algorithms/HeatDiamondAlgorithm2.hpp
deleted file mode 100644
index 2cf388ed1..000000000
--- a/src/language/algorithms/HeatDiamondAlgorithm2.hpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#ifndef HEAT_DIAMOND_ALGORITHM2_HPP
-#define HEAT_DIAMOND_ALGORITHM2_HPP
-
-#include <language/utils/FunctionSymbolId.hpp>
-#include <mesh/IMesh.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-
-#include <memory>
-#include <variant>
-#include <vector>
-
-template <size_t Dimension>
-class HeatDiamondScheme2
-{
- private:
-  class DirichletBoundaryCondition;
-  class FourierBoundaryCondition;
-  class NeumannBoundaryCondition;
-  class SymmetryBoundaryCondition;
-
- public:
-  HeatDiamondScheme2(std::shared_ptr<const IMesh> i_mesh,
-                     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-                     const FunctionSymbolId& T_id,
-                     const FunctionSymbolId& kappa_id,
-                     const FunctionSymbolId& f_id);
-};
-
-#endif   // HEAT_DIAMOND_ALGORITHM2_HPP
diff --git a/src/language/algorithms/ParabolicHeat.cpp b/src/language/algorithms/ParabolicHeat.cpp
deleted file mode 100644
index d535751ae..000000000
--- a/src/language/algorithms/ParabolicHeat.cpp
+++ /dev/null
@@ -1,568 +0,0 @@
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LinearSolver.hpp>
-#include <algebra/SmallMatrix.hpp>
-#include <algebra/TinyVector.hpp>
-#include <algebra/Vector.hpp>
-#include <language/algorithms/ParabolicHeat.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/MeshNodeBoundary.hpp>
-#include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
-#include <mesh/SubItemValuePerItem.hpp>
-#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
-#include <scheme/FourierBoundaryConditionDescriptor.hpp>
-#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
-#include <scheme/ScalarDiamondScheme.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-#include <utils/Timer.hpp>
-
-template <size_t Dimension>
-ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
-  std::shared_ptr<const IMesh> i_mesh,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-  const FunctionSymbolId& T_id,
-  const FunctionSymbolId& T_init_id,
-  const FunctionSymbolId& kappa_id,
-  const FunctionSymbolId& f_id,
-  const double& Tf,
-  const double& dt)
-{
-  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::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
-  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});
-      } else {
-        throw NotImplementedError("not implemented 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("not implemented 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);
-
-    CellValue<double> Tj =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
-    CellValue<double> Temperature =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_init_id,
-                                                                                                mesh_data.xj());
-    //{mesh->connectivity()};
-    FaceValue<double> Tl =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
-    FaceValue<double> Temperature_face =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_init_id,
-                                                                                                mesh_data.xl());
-
-    double time = 0;
-    Assert(dt > 0, "The time-step have to be positive!");
-    const CellValue<const double> primal_Vj = mesh_data.Vj();
-
-    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();
-    do {
-      double deltat = std::min(dt, Tf - time);
-      std::cout << "Current time = " << time << " time-step = " << deltat << " final time = " << Tf << "\n";
-      LegacyScalarDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, kappa_id, f_id, Temperature, Temperature_face,
-                                           Tf, deltat, w_rj, w_rl);
-      time += deltat;
-    } while (time < Tf && std::abs(time - Tf) > 1e-15);
-    {
-      Vector<double> error{mesh->numberOfCells()};
-      CellValue<double> cell_error{mesh->connectivity()};
-      Vector<double> face_error{mesh->numberOfFaces()};
-      double error_max                    = 0.;
-      size_t cell_max                     = 0;
-      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();
-        }
-      }();
-
-      parallel_for(
-        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
-          error[cell_id]      = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
-          cell_error[cell_id] = (Temperature[cell_id] - Tj[cell_id]);
-        });
-      parallel_for(
-        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-          if (primal_face_is_on_boundary[face_id]) {
-            face_error[face_id] = (Temperature_face[face_id] - Tl[face_id]) * sqrt(mes_l[face_id]);
-          } else {
-            face_error[face_id] = 0;
-          }
-        });
-      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); cell_id++) {
-        if (error_max < std::abs(cell_error[cell_id])) {
-          error_max = std::abs(cell_error[cell_id]);
-          cell_max  = cell_id;
-        }
-      }
-
-      std::cout << " ||Error||_max (cell)= " << error_max << " on cell " << cell_max << "\n";
-      std::cout << "||Error||_2 (cell)= " << std::sqrt(dot(error, error)) << "\n";
-      std::cout << "||Error||_2 (face)= " << std::sqrt(dot(face_error, face_error)) << "\n";
-      std::cout << "||Error||_2 (total)= " << std::sqrt(dot(error, error)) + std::sqrt(dot(face_error, face_error))
-                << "\n";
-    }
-  } else {
-    throw NotImplementedError("not done in 1d");
-  }
-}
-
-template <size_t Dimension>
-class ParabolicHeatScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class ParabolicHeatScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class ParabolicHeatScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class ParabolicHeatScheme<Dimension>::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;
-};
-
-template <size_t Dimension>
-class ParabolicHeatScheme<Dimension>::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:
-  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;
-  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;
-  }
-};
-
-template ParabolicHeatScheme<1>::ParabolicHeatScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const double&,
-  const double&);
-
-template ParabolicHeatScheme<2>::ParabolicHeatScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const double&,
-  const double&);
-
-template ParabolicHeatScheme<3>::ParabolicHeatScheme(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const double&,
-  const double&);
diff --git a/src/language/algorithms/ParabolicHeat.hpp b/src/language/algorithms/ParabolicHeat.hpp
deleted file mode 100644
index 84ee713ee..000000000
--- a/src/language/algorithms/ParabolicHeat.hpp
+++ /dev/null
@@ -1,33 +0,0 @@
-#ifndef PARABOLIC_HEAT_HPP
-#define PARABOLIC_HEAT_HPP
-
-#include <language/utils/FunctionSymbolId.hpp>
-#include <mesh/IMesh.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-
-#include <memory>
-#include <variant>
-#include <vector>
-
-template <size_t Dimension>
-class ParabolicHeatScheme
-{
- private:
-  class DirichletBoundaryCondition;
-  class FourierBoundaryCondition;
-  class NeumannBoundaryCondition;
-  class SymmetryBoundaryCondition;
-  class InterpolationWeightsManager;
-
- public:
-  ParabolicHeatScheme(std::shared_ptr<const IMesh> i_mesh,
-                      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-                      const FunctionSymbolId& T_id,
-                      const FunctionSymbolId& T_init_id,
-                      const FunctionSymbolId& kappa_id,
-                      const FunctionSymbolId& f_id,
-                      const double& final_time,
-                      const double& dt);
-};
-
-#endif   // HEAT_DIAMOND_ALGORITHM_HPP
diff --git a/src/language/algorithms/UnsteadyElasticity.cpp b/src/language/algorithms/UnsteadyElasticity.cpp
deleted file mode 100644
index 6eb2017c1..000000000
--- a/src/language/algorithms/UnsteadyElasticity.cpp
+++ /dev/null
@@ -1,613 +0,0 @@
-#include <language/algorithms/UnsteadyElasticity.hpp>
-
-#include <algebra/CRSMatrix.hpp>
-#include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/LeastSquareSolver.hpp>
-#include <algebra/LinearSolver.hpp>
-#include <algebra/SmallMatrix.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/PrimalToDiamondDualConnectivityDataMapper.hpp>
-#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
-#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-#include <scheme/VectorDiamondScheme.hpp>
-
-template <size_t Dimension>
-UnsteadyElasticity<Dimension>::UnsteadyElasticity(
-  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,
-  const FunctionSymbolId& U_id,
-  const FunctionSymbolId& U_init_id,
-  const double& Tf,
-  const double& dt)
-{
-  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) {
-        MeshFaceBoundary mesh_face_boundary = getMeshFaceBoundary(*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 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("Dirichlet conditions are not supported in 1d");
-        }
-      } else if (dirichlet_bc_descriptor.name() == "normal_strain") {
-        if constexpr (Dimension > 1) {
-          MeshFaceBoundary 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;
-    }();
-
-    FaceValue<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 auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
-    const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
-
-    {
-      CellValue<TinyVector<Dimension>> velocity = InterpolateItemValue<TinyVector<Dimension>(
-        TinyVector<Dimension>)>::template interpolate<ItemType::cell>(U_id, mesh_data.xj());
-      CellValue<TinyVector<Dimension>> Uj       = InterpolateItemValue<TinyVector<Dimension>(
-        TinyVector<Dimension>)>::template interpolate<ItemType::cell>(U_init_id, mesh_data.xj());
-
-      CellValue<TinyVector<Dimension>> fj            = InterpolateItemValue<TinyVector<Dimension>(
-        TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
-      FaceValue<TinyVector<Dimension>> velocity_face = InterpolateItemValue<TinyVector<Dimension>(
-        TinyVector<Dimension>)>::template interpolate<ItemType::face>(U_id, mesh_data.xl());
-      FaceValue<TinyVector<Dimension>> Ul            = InterpolateItemValue<TinyVector<Dimension>(
-        TinyVector<Dimension>)>::template interpolate<ItemType::face>(U_init_id, mesh_data.xl());
-
-      double time = 0;
-      Assert(dt > 0, "The time-step have to be positive!");
-      const CellValue<const double> primal_Vj = mesh_data.Vj();
-
-      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();
-      do {
-        double deltat = std::min(dt, Tf - time);
-        std::cout << "Current time = " << time << " time-step = " << deltat << " final time = " << Tf << "\n";
-        LegacyVectorDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, Uj, Ul, Tf, deltat,
-                                             w_rj, w_rl);
-        time += deltat;
-      } while (time < Tf && std::abs(time - Tf) > 1e-15);
-
-      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();
-        }
-      }();
-
-      Vector<double> Uexacte{mesh->numberOfCells() * Dimension};
-      for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
-        for (size_t l = 0; l < Dimension; ++l) {
-          Uexacte[(cell_dof_number[j] * Dimension) + l] = velocity[j][l];
-        }
-      }
-
-      Vector<double> error{mesh->numberOfCells() * Dimension};
-      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-        for (size_t i = 0; i < Dimension; ++i) {
-          error[(cell_id * Dimension) + i] = (Uj[cell_id][i] - velocity[cell_id][i]) * sqrt(primal_Vj[cell_id]);
-        }
-      }
-      Vector<double> error_face{mesh->numberOfFaces() * Dimension};
-      parallel_for(
-        mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-          if (primal_face_is_on_boundary[face_id]) {
-            for (size_t i = 0; i < Dimension; ++i) {
-              error_face[face_id * Dimension + i] = (Ul[face_id][i] - velocity_face[face_id][i]) * sqrt(mes_l[face_id]);
-            }
-          } else {
-            error_face[face_id] = 0;
-          }
-        });
-
-      std::cout << "||Error||_2 (cell)= " << std::sqrt(dot(error, error)) << "\n";
-      std::cout << "||Error||_2 (face)= " << std::sqrt(dot(error_face, error_face)) << "\n";
-      std::cout << "||Error||_2 (total)= " << std::sqrt(dot(error, error)) + std::sqrt(dot(error_face, error_face))
-                << "\n";
-
-      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 <size_t Dimension>
-class UnsteadyElasticity<Dimension>::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;
-};
-
-template <size_t Dimension>
-class UnsteadyElasticity<Dimension>::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;
-};
-
-template <size_t Dimension>
-class UnsteadyElasticity<Dimension>::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;
-};
-
-template <size_t Dimension>
-class UnsteadyElasticity<Dimension>::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;
-  FaceValue<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<bool> primal_face_is_on_boundary,
-                              NodeValue<bool> primal_node_is_on_boundary,
-                              FaceValue<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;
-  }
-};
-
-template UnsteadyElasticity<1>::UnsteadyElasticity(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const double&,
-  const double&);
-
-template UnsteadyElasticity<2>::UnsteadyElasticity(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const double&,
-  const double&);
-
-template UnsteadyElasticity<3>::UnsteadyElasticity(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const double&,
-  const double&);
diff --git a/src/language/algorithms/UnsteadyElasticity.hpp b/src/language/algorithms/UnsteadyElasticity.hpp
deleted file mode 100644
index f55ef9b32..000000000
--- a/src/language/algorithms/UnsteadyElasticity.hpp
+++ /dev/null
@@ -1,36 +0,0 @@
-#ifndef UNSTEADY_ELASTICITY_HPP
-#define UNSTEADY_ELASTICITY_HPP
-
-#include <algebra/TinyVector.hpp>
-#include <language/utils/FunctionSymbolId.hpp>
-#include <memory>
-#include <mesh/IMesh.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <variant>
-#include <vector>
-
-template <size_t Dimension>
-class UnsteadyElasticity
-{
- private:
-  class DirichletBoundaryCondition;
-  class NormalStrainBoundaryCondition;
-  class SymmetryBoundaryCondition;
-  class InterpolationWeightsManager;
-
- public:
-  UnsteadyElasticity(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,
-                     const FunctionSymbolId& U_id,
-                     const FunctionSymbolId& U_init_id,
-                     const double& Tf,
-                     const double& dt);
-};
-
-#endif   // ELASTICITY_DIAMOND_ALGORITHM2_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index c4d28fe5d..c1c8eebef 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -3,12 +3,6 @@
 #include <analysis/GaussLegendreQuadratureDescriptor.hpp>
 #include <analysis/GaussLobattoQuadratureDescriptor.hpp>
 #include <analysis/GaussQuadratureDescriptor.hpp>
-#include <language/algorithms/ElasticityDiamondAlgorithm.hpp>
-#include <language/algorithms/Heat5PointsAlgorithm.hpp>
-#include <language/algorithms/HeatDiamondAlgorithm.hpp>
-#include <language/algorithms/HeatDiamondAlgorithm2.hpp>
-#include <language/algorithms/ParabolicHeat.hpp>
-#include <language/algorithms/UnsteadyElasticity.hpp>
 #include <language/modules/BinaryOperatorRegisterForVh.hpp>
 #include <language/modules/MathFunctionRegisterForVh.hpp>
 #include <language/modules/UnaryOperatorRegisterForVh.hpp>
@@ -440,67 +434,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("heat", std::function(
-
-                                      [](std::shared_ptr<const IMesh> p_mesh,
-                                         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                           bc_descriptor_list,
-                                         const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
-                                         const FunctionSymbolId& f_id) -> void {
-                                        switch (p_mesh->dimension()) {
-                                        case 1: {
-                                          HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                          break;
-                                        }
-                                        case 2: {
-                                          HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                          break;
-                                        }
-                                        case 3: {
-                                          HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                          break;
-                                        }
-                                        default: {
-                                          throw UnexpectedError("invalid mesh dimension");
-                                        }
-                                        }
-                                      }
-
-                                      ));
-
-  this->_addBuiltinFunction("parabolicheat",
-                            std::function(
-
-                              [](std::shared_ptr<const IMesh> p_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const FunctionSymbolId& T_id, const FunctionSymbolId& T_init_id,
-                                 const FunctionSymbolId& kappa_id, const FunctionSymbolId& f_id,
-                                 const double& final_time, const double& dt) -> void {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  ParabolicHeatScheme<1>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
-                                                         f_id,   final_time,         dt};
-                                  break;
-                                }
-                                case 2: {
-                                  ParabolicHeatScheme<2>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
-                                                         f_id,   final_time,         dt};
-                                  break;
-                                }
-                                case 3: {
-                                  ParabolicHeatScheme<3>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
-                                                         f_id,   final_time,         dt};
-                                  break;
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
-                              }
-
-                              ));
-
   this->_addBuiltinFunction(
     "parabolicheat",
     std::function(
@@ -515,40 +448,6 @@ SchemeModule::SchemeModule()
 
       ));
 
-  this->_addBuiltinFunction("unsteadyelasticity",
-                            std::function(
-
-                              [](std::shared_ptr<const IMesh> p_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id,
-                                 const FunctionSymbolId& f_id, const FunctionSymbolId& U_id,
-                                 const FunctionSymbolId& U_init_id, const double& final_time,
-                                 const double& dt) -> void {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  UnsteadyElasticity<1>{p_mesh, bc_descriptor_list, lambda_id,  mu_id, f_id,
-                                                        U_id,   U_init_id,          final_time, dt};
-                                  break;
-                                }
-                                case 2: {
-                                  UnsteadyElasticity<2>{p_mesh, bc_descriptor_list, lambda_id,  mu_id, f_id,
-                                                        U_id,   U_init_id,          final_time, dt};
-                                  break;
-                                }
-                                case 3: {
-                                  UnsteadyElasticity<3>{p_mesh, bc_descriptor_list, lambda_id,  mu_id, f_id,
-                                                        U_id,   U_init_id,          final_time, dt};
-                                  break;
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
-                              }
-
-                              ));
-
   this->_addBuiltinFunction("unsteadyelasticity",
                             std::function(
 
@@ -603,92 +502,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("heat2", std::function(
-
-                                       [](std::shared_ptr<const IMesh> p_mesh,
-                                          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                            bc_descriptor_list,
-                                          const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
-                                          const FunctionSymbolId& f_id) -> void {
-                                         switch (p_mesh->dimension()) {
-                                         case 1: {
-                                           HeatDiamondScheme2<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                           break;
-                                         }
-                                         case 2: {
-                                           HeatDiamondScheme2<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                           break;
-                                         }
-                                         case 3: {
-                                           HeatDiamondScheme2<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                           break;
-                                         }
-                                         default: {
-                                           throw UnexpectedError("invalid mesh dimension");
-                                         }
-                                         }
-                                       }
-
-                                       ));
-
-  this->_addBuiltinFunction("elasticity",
-                            std::function(
-
-                              [](std::shared_ptr<const IMesh> p_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id,
-                                 const FunctionSymbolId& f_id, const FunctionSymbolId& U_id) -> void {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  ElasticityDiamondScheme<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
-                                  break;
-                                }
-                                case 2: {
-                                  ElasticityDiamondScheme<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
-                                  break;
-                                }
-                                case 3: {
-                                  ElasticityDiamondScheme<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
-                                  break;
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("heat5points",
-                            std::function(
-
-                              [](std::shared_ptr<const IMesh> p_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
-                                 const FunctionSymbolId& f_id) -> void {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  Heat5PointsAlgorithm<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                  break;
-                                }
-                                case 2: {
-                                  Heat5PointsAlgorithm<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                  break;
-                                }
-                                case 3: {
-                                  Heat5PointsAlgorithm<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
-                                  break;
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
-                              }
-
-                              ));
-
   this->_addBuiltinFunction("lagrangian",
                             std::function(
 
-- 
GitLab