First step to the implementation of the nodal solver for diffusion

src/language/modules/SchemeModule.cpp

@@ -42,6 +42,7 @@
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
 #include <scheme/ScalarDiamondScheme.hpp>
+#include <scheme/ScalarNodalScheme.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VectorDiamondScheme.hpp>
 #include <utils/Socket.hpp>
@@ -601,6 +602,24 @@ SchemeModule::SchemeModule()
         return ScalarDiamondSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
       }
 
+      ));
+  this->_addBuiltinFunction(
+    "nodalheat",
+    std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+      const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+                               const std::shared_ptr<const IDiscreteFunction>&,
+                               const std::shared_ptr<const IDiscreteFunction>&,
+                               const std::shared_ptr<const IDiscreteFunction>&,
+                               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&)>>(
+
+      [](const std::shared_ptr<const IDiscreteFunction>& alpha,
+         const std::shared_ptr<const IDiscreteFunction>& mub_dual,
+         const std::shared_ptr<const IDiscreteFunction>& mu_dual, const std::shared_ptr<const IDiscreteFunction>& f,
+         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+        -> const std::shared_ptr<const IDiscreteFunction> {
+        return ScalarNodalSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
+      }
+
       ));
 
   this->_addBuiltinFunction("unsteadyelasticity",

src/scheme/ScalarNodalScheme.cpp

@@ -3,7 +3,6 @@
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 
-
 class ScalarNodalSchemeHandler::IScalarNodalScheme
 {
  public:
@@ -56,6 +55,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
    private:
     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>&
@@ -64,14 +64,22 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       return m_face_list;
     }
 
+    const Array<const NodeId>&
+    nodeList() const
+    {
+      return m_node_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}
+    NeumannBoundaryCondition(const Array<const FaceId>& face_list,
+                             const Array<const NodeId>& node_list,
+                             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());
     }
@@ -118,7 +126,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     ~FourierBoundaryCondition() = default;
   };
 
-
  public:
   std::shared_ptr<const IDiscreteFunction>
   getSolution() const final
@@ -133,9 +140,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                     const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
                     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
   {
-
-    using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition,
-                                           NeumannBoundaryCondition>;
+    using BoundaryCondition =
+      std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition>;
 
     using BoundaryConditionList = std::vector<BoundaryCondition>;
 
@@ -179,6 +185,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         if constexpr (Dimension > 1) {
           MeshFaceBoundary<Dimension> mesh_face_boundary =
             getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
+          MeshNodeBoundary<Dimension> mesh_node_boundary =
+            getMeshNodeBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
 
           const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
           MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
@@ -189,10 +197,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                                                                                                       mesh_face_boundary
                                                                                                         .faceList());
 
-          boundary_condition_list.push_back(NeumannBoundaryCondition{mesh_face_boundary.faceList(), value_list});
+          boundary_condition_list.push_back(
+            NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
 
         } else {
-          throw NotImplementedError("Dirichlet BC in 1d");
+          throw NotImplementedError("Neumann BC in 1d");
         }
         break;
       }
@@ -207,89 +216,119 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       }
     }
 
-    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();
+    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 NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
+
+    const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
+
+    // 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();
+    const auto& node_to_cell_matrix               = mesh->connectivity().nodeToCellMatrix();
 
-      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());
+    const NodeValue<const bool> node_is_neumann = [&] {
+      NodeValue<bool> compute_node_is_neumann{mesh->connectivity()};
+      compute_node_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, DirichletBoundaryCondition>)) {
+            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];
-                  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++;
+                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;
                 }
               }
             }
           },
           boundary_condition);
       }
-
-        return compute_face_dof_number;
+      return compute_node_is_neumann;
     }();
 
-      const FaceValue<const bool> face_is_neumann = [&] {
-        FaceValue<bool> face_is_neumann{mesh->connectivity()};
-        face_is_neumann.fill(false);
+    const NodeValue<const bool> node_is_dirichlet = [&] {
+      NodeValue<bool> compute_node_is_dirichlet{mesh->connectivity()};
+      compute_node_is_dirichlet.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>)) {
+            if constexpr (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];
-                  face_is_neumann[face_id] = true;
+                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_dirichlet[node_id] = true;
+                }
               }
             }
           },
           boundary_condition);
       }
-
-        return face_is_neumann;
+      return compute_node_is_dirichlet;
     }();
 
+    const NodeValue<const bool> node_is_corner = [&] {
+      NodeValue<bool> compute_node_is_corner{mesh->connectivity()};
+      compute_node_is_corner.fill(false);
+      parallel_for(
+        mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+          if (is_boundary_node[node_id]) {
+            const auto& node_to_cell = node_to_cell_matrix[node_id];
 
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-      const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
-      const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-      const NodeValue<const TinyVector<Dimension>>& xr = mesh_data.xr();
-
-      const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
-      const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-      const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
-
-      const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
-
-      {
+            for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+              const CellId cell_id      = node_to_cell[i_cell];
+              const auto& cell_nodes    = cell_to_node_matrix[cell_id];
+              const NodeId prev_node_id = cell_to_node_matrix[cell_id][(node_id - 1) % cell_nodes.size()];
+              const NodeId next_node_id = cell_to_node_matrix[cell_id][(node_id + 1) % cell_nodes.size()];
 
+              if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
+                compute_node_is_corner[node_id] = true;
+              }
+            }
+          }
+        });
 
-        CellValue<const TinyMatrix<Dimension>> cell_kappaj  = cell_k->cellValues();
-        CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
+      return compute_node_is_corner;
+    }();
 
-        const CellValue<const double> Vj = mesh_data.Vj();
-        const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+    const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] {
+      NodeValue<TinyVector<Dimension>> compute_exterior_normal;
+      parallel_for(
+        mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+          compute_exterior_normal[node_id]          = zero;
+          const auto& node_to_cell                  = node_to_cell_matrix[node_id];
+          const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+          for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+            const CellId cell_id          = node_to_cell[i_cell];
+            const unsigned int i_node_loc = node_local_number_in_its_cell[cell_id];
+            compute_exterior_normal[node_id] += Cjr(cell_id, i_node_loc);
+          }
+          compute_exterior_normal[node_id] =
+            1. / l2Norm(compute_exterior_normal[node_id]) * compute_exterior_normal[node_id];
+        });
+      return compute_exterior_normal;
+    }();
 
     const FaceValue<const double> mes_l = [&] {
       if constexpr (Dimension == 1) {
@@ -301,8 +340,58 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       }
     }();
 
+    const NodeValue<const double> node_boundary_values = [&] {
+      NodeValue<double> compute_node_boundary_values;
+      compute_node_boundary_values.fill(0);
+      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>) {
+              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];
+                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 (node_is_dirichlet[node_id] == false) {
+                    compute_node_boundary_values[node_id] += 0.5 * value_list[i_face] * mes_l[face_id];
+                  }
+                }
+              }
+            } else 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];
+                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 (node_is_neumann[node_id] == false) {
+                    compute_node_boundary_values[node_id] += 0.5 * value_list[i_face] * mes_l[face_id];
+                  } else {
+                    compute_node_boundary_values[node_id] = value_list[i_face] * mes_l[face_id];
+                  }
+                }
+              }
+            }
+          },
+          boundary_condition);
+      }
+      return compute_node_boundary_values;
+    }();
+
+    {
+      CellValue<const TinyMatrix<Dimension>> cell_kappaj  = cell_k->cellValues();
+      CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
+
       const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
-          NodeValue<const TinyMatrix<Dimension>> kappa;
+        NodeValue<TinyMatrix<Dimension>> kappa;
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             kappa[node_id]           = zero;
@@ -310,18 +399,17 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             double weight            = 0;
             for (size_t j = 0; j < node_to_cell.size(); ++j) {
               const CellId J      = node_to_cell[j];
-                  double local_weight = 1 / l2Norm(xr[node_id]-xj[J]);
-                  kappa[node_id] += cell_kappaj[J] * local_weight;
+              double local_weight = 1. / l2Norm(xr[node_id] - xj[J]);
+              kappa[node_id] += local_weight * cell_kappaj[J];
               weight += local_weight;
             }
-                kappa[node_id] /= weight;
-            }
-          );
+            kappa[node_id] = 1. / weight * kappa[node_id];
+          });
         return kappa;
       }();
 
       const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] {
-          NodeValue<const TinyMatrix<Dimension>> kappa;
+        NodeValue<TinyMatrix<Dimension>> kappa;
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             kappa[node_id]           = zero;
@@ -329,17 +417,17 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             double weight            = 0;
             for (size_t j = 0; j < node_to_cell.size(); ++j) {
               const CellId J      = node_to_cell[j];
-                  double local_weight = 1 / l2Norm(xr[node_id]-xj[J]);
-                  kappa[node_id] += cell_kappajb[J] * local_weight;
+              double local_weight = 1. / l2Norm(xr[node_id] - xj[J]);
+              kappa[node_id] += local_weight * cell_kappajb[J];
               weight += local_weight;
             }
-                kappa[node_id] /= weight;
+            kappa[node_id] = 1. / weight * kappa[node_id];
           });
         return kappa;
       }();
 
       const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
-          NodeValue<const TinyMatrix<Dimension>> beta;
+        NodeValue<TinyMatrix<Dimension>> beta;
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
@@ -354,111 +442,142 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         return beta;
       }();
 
-
-        FaceValue<const double> alpha_l = [&] {
-          FaceValue<double> alpha_j{mesh->connectivity()};
-
+      const NodeValue<const TinyMatrix<Dimension>> kappar_invBetar = [&] {
+        NodeValue<TinyMatrix<Dimension>> kappa_invBeta;
         parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-              alpha_j[face_id] = 1;
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            kappa_invBeta[node_id] = node_kappar[node_id] * inverse(node_betar[node_id]);
           });
-
-          return alpha_j;
+        return kappa_invBeta;
       }();
 
-        FaceValue<const double> alphab_l = [&] {
-          FaceValue<double> alpha_lb{mesh->connectivity()};
+      const Array<int> non_zeros{mesh->numberOfCells()};
+      non_zeros.fill(Dimension);   // Modif pour optimiser
+      CRSMatrixDescriptor<double> S(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
 
-          parallel_for(
-            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-              alpha_lb[face_id] = 1; //Refaire
-            });
+      for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+        const auto& node_to_cell = node_to_cell_matrix[node_id];
 
-          return alpha_lb;
-        }();
+        if (not is_boundary_node[node_id]) {
+          for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
+            const CellId J                          = node_to_cell[cell_id_j];
+            const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
+            const size_t node_id_j                  = node_local_number_in_j_cell[node_id];
 
+            for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
+              const CellId K                          = node_to_cell[cell_id_k];
+              const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
+              const size_t node_id_k                  = node_local_number_in_k_cell[node_id];
 
-        const Array<int> non_zeros{number_of_dof};
-        non_zeros.fill(Dimension); //Modif pour optimiser
-        CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zeros);
-
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          const auto& face_to_cell = face_to_cell_matrix[face_id];
-          const double beta_l             = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
-          const double betab_l            = l2Norm(dualClj[face_id]) * alphab_l[face_id] * mes_l[face_id];
-          for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
-            const CellId cell_id1 = face_to_cell[i_cell];
-            for (size_t j_cell = 0; j_cell < face_to_cell.size(); ++j_cell) {
-              const CellId cell_id2 = face_to_cell[j_cell];
-              if (i_cell == j_cell) {
-                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
-                if (face_is_neumann[face_id]) {
-                  S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= betab_l;
-                }
-              } else {
-                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
+              S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), Cjr(J, node_id_j));
             }
           }
+        } else if (not node_is_corner[node_id]) {
         }
       }
 
-        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          const size_t j = cell_dof_number[cell_id];
-          S(j, j) += (*alpha)[cell_id] * Vj[cell_id];
+      for (const auto& boundary_condition : boundary_condition_list) {
+        std::visit(
+          [&](auto&& bc) {
+            using T = std::decay_t<decltype(bc)>;
+            if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+              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];
+                const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+                for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
+                  const CellId J                          = node_to_cell[cell_id_j];
+                  const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
+                  const size_t node_id_j                  = node_local_number_in_j_cell[node_id];
+
+                  for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
+                    const CellId K                          = node_to_cell[cell_id_k];
+                    const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
+                    const size_t node_id_k                  = node_local_number_in_k_cell[node_id];
+
+                    const TinyMatrix<Dimension> I   = identity;
+                    const TinyVector<Dimension> n   = exterior_normal[node_id];
+                    const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
+                    const TinyMatrix<Dimension> P   = I - nxn;
+
+                    S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), P * Cjr(J, node_id_j));
+                  }
+                }
               }
+            } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
+              const auto& node_list = bc.faceList();
 
+              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                const NodeId node_id     = node_list[i_node];
+                const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+                for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
+                  const CellId J                          = node_to_cell[cell_id_j];
+                  const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
+                  const size_t node_id_j                  = node_local_number_in_j_cell[node_id];
+
+                  for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
+                    const CellId K                          = node_to_cell[cell_id_k];
+                    const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
+                    const size_t node_id_k                  = node_local_number_in_k_cell[node_id];
 
-        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
-          if (face_is_dirichlet[face_id]) {
-            S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
+                    S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), Cjr(J, node_id_j));
+                  }
+                }
               }
             }
+          },
+          boundary_condition);
+      }
+
+      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+        S(cell_id, cell_id) += Vj[cell_id];
+      }
 
       CellValue<const double> fj = f->cellValues();
 
       CRSMatrix A{S.getCRSMatrix()};
-        Vector<double> b{number_of_dof};
+      Vector<double> b{mesh->numberOfCells()};
       b = zero;
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          b[cell_dof_number[cell_id]] = fj[cell_id] * Vj[cell_id];
+        b[cell_id] = fj[cell_id] * 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>) {
+            if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {   // To do
               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];
+                b[face_id] += value_list[i_face];
               }
-                }
-              },
-              boundary_condition);
-          }
-        }
-        // EL b^L
-        for (const auto& boundary_condition : boundary_condition_list) {
-          std::visit(
-            [&](auto&& bc) {
-              using T = std::decay_t<decltype(bc)>;
-              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
-                            (std::is_same_v<T, FourierBoundaryCondition>)) {
-                const auto& face_list  = bc.faceList();
+            } else if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+              const auto& node_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];
+
+              for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                NodeId node_id                = node_list[i_node];
+                const auto& node_to_cell      = node_to_cell_matrix[node_id];
+                const TinyVector<Dimension> n = exterior_normal[node_id];
+
+                for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+                  CellId J                              = node_to_cell[i_cell];
+                  const auto& node_local_number_in_cell = node_local_numbers_in_their_cells.itemArray(J);
+                  const size_t node_id_j                = node_local_number_in_cell[node_id];
+                  b[J] -= dot(Cjr(J, node_id_j), n) * value_list[i_node];
+                }
               }
             }
           },
           boundary_condition);
       }
 
-        Vector<double> T{number_of_dof};
+      Vector<double> T{mesh->numberOfCells()};
       T = zero;
 
       LinearSolver solver;
@@ -467,8 +586,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       m_solution     = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
       auto& solution = *m_solution;
       parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_dof_number[cell_id]]; });
-      }
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; });
     }
   }
 };
@@ -481,7 +599,7 @@ ScalarNodalSchemeHandler::solution() const
 
 ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
   const std::shared_ptr<const IDiscreteFunction>& alpha,
-  const std::shared_ptr<const IDiscreteFunction>& k_bound,
+  const std::shared_ptr<const IDiscreteFunction>& cell_k_bound,
   const std::shared_ptr<const IDiscreteFunction>& cell_k,
   const std::shared_ptr<const IDiscreteFunction>& f,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
@@ -490,30 +608,12 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
   if (not i_mesh) {
     throw NormalError("primal discrete functions are not defined on the same mesh");
   }
-  const std::shared_ptr i_dual_mesh = getCommonMesh({cell_k_bound, cell_k});
-  if (not i_dual_mesh) {
-    throw NormalError("dual discrete functions are not defined on the same mesh");
-  }
+
   checkDiscretizationType({alpha, cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
 
   switch (i_mesh->dimension()) {
   case 1: {
-    using MeshType                   = Mesh<Connectivity<1>>;
-    using DiscreteTensorFunctionType = DiscreteFunctionP0<1, TinyMatrix<1>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<1, double>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
-    m_scheme =
-      std::make_unique<ScalarNodalScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
-                                               std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
-                                               std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
-                                               bc_descriptor_list);
+    throw NormalError("The scheme is not implemented in dimension 1");
     break;
   }
   case 2: {
@@ -523,10 +623,6 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
     m_scheme =
       std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
                                              std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
@@ -542,10 +638,6 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
-    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
-      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
-    }
-
     m_scheme =
       std::make_unique<ScalarNodalScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
                                              std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),