diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 4f9719f639e6181fafecb8eb02f1ae3919572a63..4639f58dd3182e796c11c0faefe062bb71d55ab0 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -53,7 +53,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                                const Array<const double>& value_list)
       : m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list}
     {
-      Assert(m_value_list.size() == m_face_list.size());
+      Assert(m_value_list.size() == m_node_list.size());
     }
 
     ~DirichletBoundaryCondition() = default;
@@ -90,7 +90,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                              const Array<const double>& value_list)
       : m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list}
     {
-      Assert(m_value_list.size() == m_face_list.size());
+      Assert(m_value_list.size() == m_node_list.size());
     }
 
     ~NeumannBoundaryCondition() = default;
@@ -101,13 +101,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
    private:
     const Array<const double> m_coef_list;
     const Array<const double> m_value_list;
-    const Array<const FaceId> m_face_list;
+    const Array<const NodeId> m_node_list;
 
    public:
-    const Array<const FaceId>&
-    faceList() const
+    const Array<const NodeId>&
+    nodeList() const
     {
-      return m_face_list;
+      return m_node_list;
     }
 
     const Array<const double>&
@@ -123,13 +123,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     }
 
    public:
-    FourierBoundaryCondition(const Array<const FaceId>& face_list,
+    FourierBoundaryCondition(const Array<const NodeId>& node_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}
+      : m_coef_list{coef_list}, m_value_list{value_list}, m_node_list{node_list}
     {
-      Assert(m_coef_list.size() == m_face_list.size());
-      Assert(m_value_list.size() == m_face_list.size());
+      Assert(m_coef_list.size() == m_node_list.size());
+      Assert(m_value_list.size() == m_node_list.size());
     }
 
     ~FourierBoundaryCondition() = default;
@@ -173,13 +173,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             getMeshNodeBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
 
           const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-          MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
+          // 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());
+            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
+                                                                                                      mesh_node_boundary
+                                                                                                        .nodeList());
           boundary_condition_list.push_back(
             DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
 
@@ -203,13 +202,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             getMeshNodeBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
 
           const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
-          MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
+          //          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());
+            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
+                                                                                                      mesh_node_boundary
+                                                                                                        .nodeList());
 
           boundary_condition_list.push_back(
             NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
@@ -246,13 +244,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             getMeshNodeBoundary(*mesh, dirichlet_bc_prev_descriptor.boundaryDescriptor());
 
           const FunctionSymbolId g_id = dirichlet_bc_prev_descriptor.rhsSymbolId();
-          MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
+          //          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());
+            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
+                                                                                                      mesh_node_boundary
+                                                                                                        .nodeList());
           boundary_prev_condition_list.push_back(
             DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
 
@@ -276,13 +273,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             getMeshNodeBoundary(*mesh, neumann_bc_prev_descriptor.boundaryDescriptor());
 
           const FunctionSymbolId g_id = neumann_bc_prev_descriptor.rhsSymbolId();
-          MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
+          //          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());
+            InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
+                                                                                                      mesh_node_boundary
+                                                                                                        .nodeList());
 
           boundary_prev_condition_list.push_back(
             NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
@@ -315,8 +311,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
     const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
 
-    const auto& node_to_face_matrix               = mesh->connectivity().nodeToFaceMatrix();
-    const auto& face_to_node_matrix               = mesh->connectivity().faceToNodeMatrix();
+    const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
+    // const auto& face_to_node_matrix               = mesh->connectivity().faceToNodeMatrix();
     const auto& cell_to_node_matrix               = mesh->connectivity().cellToNodeMatrix();
     const auto& node_local_numbers_in_their_cells = mesh->connectivity().nodeLocalNumbersInTheirCells();
     const CellValue<const double> Vj              = mesh_data.Vj();
@@ -331,16 +327,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
           [&](auto&& bc) {
             using T = std::decay_t<decltype(bc)>;
             if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
-              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];
-                const auto& face_nodes = face_to_node_matrix[face_id];
-
-                for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                  const NodeId node_id = face_nodes[i_node];
-
-                  compute_node_is_neumann[node_id] = true;
-                }
+              const auto& node_list = bc.nodeList();
+              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                const NodeId node_id             = node_list[i_node];
+                compute_node_is_neumann[node_id] = true;
               }
             }
           },
@@ -357,17 +347,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
           [&](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& node_list = bc.nodeList();
 
-              for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                const FaceId face_id   = face_list[i_face];
-                const auto& face_nodes = face_to_node_matrix[face_id];
+              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                const NodeId node_id = node_list[i_node];
 
-                for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                  const NodeId node_id = face_nodes[i_node];
-
-                  compute_node_is_dirichlet[node_id] = true;
-                }
+                compute_node_is_dirichlet[node_id] = true;
               }
             }
           },
@@ -486,57 +471,28 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
           [&](auto&& bc) {
             using T = std::decay_t<decltype(bc)>;
             if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
-              const auto& face_list  = bc.faceList();
+              const auto& node_list  = bc.nodeList();
               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];
-                const auto& face_nodes = face_to_node_matrix[face_id];
-
-                for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                  const NodeId node_id = face_nodes[i_node];
-                  if (not node_is_dirichlet[node_id]) {
-                    if (node_is_angle[node_id]) {
-                      compute_node_boundary_values[node_id] +=
-                        value_list[i_face] * std::abs(dot(nl[face_id], exterior_normal[node_id]));
-                    } else {
-                      compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
-                      sum_mes_l[node_id] += mes_l[face_id];
-                    }
-                  }
-                }
+              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                const NodeId node_id                  = node_list[i_node];
+                compute_node_boundary_values[node_id] = value_list[i_node];
               }
+
             } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
-              const auto& face_list  = bc.faceList();
+              const auto& node_list  = bc.nodeList();
               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];
-                const auto& face_nodes = face_to_node_matrix[face_id];
-
-                for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                  const NodeId node_id = face_nodes[i_node];
-                  if (not node_is_neumann[node_id]) {
-                    compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
-                    // compute_node_boundary_values[node_id] = xr[node_id][0];
-                    sum_mes_l[node_id] += mes_l[face_id];
-                  } else {
-                    compute_node_boundary_values[node_id] = value_list[i_face];
-                  }
-                }
+              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                const NodeId node_id = node_list[i_node];
+
+                compute_node_boundary_values[node_id] = value_list[i_node];
               }
             }
           },
           boundary_condition);
       }
-      for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-        if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
-          compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
-        } else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
-          compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
-        }
-        // std::cout << node_id << "  " << compute_node_boundary_values[node_id] << "\n";
-      }
+      // std::cout << node_id << "  " << compute_node_boundary_values[node_id] << "\n";
       return compute_node_boundary_values;
     }();
 
@@ -550,51 +506,27 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
           [&](auto&& bc) {
             using T = std::decay_t<decltype(bc)>;
             if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
-              const auto& face_list  = bc.faceList();
+              const auto& node_list  = bc.nodeList();
               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];
-                const auto& face_nodes = face_to_node_matrix[face_id];
+              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                const NodeId node_id = node_list[i_node];
 
-                for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                  const NodeId node_id = face_nodes[i_node];
-                  if (not node_is_dirichlet[node_id]) {
-                    compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
-                    sum_mes_l[node_id] += mes_l[face_id];
-                  }
-                }
+                compute_node_boundary_values[node_id] = value_list[i_node];
               }
-
             } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
-              const auto& face_list  = bc.faceList();
+              const auto& node_list  = bc.nodeList();
               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];
-                const auto& face_nodes = face_to_node_matrix[face_id];
-
-                for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                  const NodeId node_id = face_nodes[i_node];
-                  if (not node_is_neumann[node_id]) {
-                    compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
-                    sum_mes_l[node_id] += mes_l[face_id];
-                  } else {
-                    compute_node_boundary_values[node_id] = value_list[i_face];
-                  }
-                }
+              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                const NodeId node_id = node_list[i_node];
+
+                compute_node_boundary_values[node_id] = value_list[i_node];
               }
             }
           },
           boundary_condition);
       }
-      for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-        if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id])) {
-          compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
-        } else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
-          compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
-        }
-      }
       return compute_node_boundary_values;
     }();