diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index edbce1e3f20264e1a5da949314b1f2a4123c133e..8800b37b4ba8dd1d69387f3046e8e5cced712e77 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -603,24 +603,23 @@ SchemeModule::SchemeModule()
       }
 
       ));
-  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>>&)>>(
+  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::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();
-      }
+                              [](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{mub_dual, mu_dual, f, bc_descriptor_list}.solution();
+                              }
 
-      ));
+                              ));
 
   this->_addBuiltinFunction("unsteadyelasticity",
                             std::make_shared<BuiltinFunctionEmbedder<
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 215cf8510ec40f3cf19268d70c561f1034fc6ac6..d1f22b1442754f35a06f8e7bae2e1a67d149bdf6 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -27,8 +27,15 @@ 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 NodeId>&
+    nodeList() const
+    {
+      return m_node_list;
+    }
+
     const Array<const FaceId>&
     faceList() const
     {
@@ -41,8 +48,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       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}
+    DirichletBoundaryCondition(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());
     }
@@ -134,7 +143,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
   }
 
   ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
-                    const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
                     const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_bound,
                     const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
                     const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
@@ -157,6 +165,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         if constexpr (Dimension > 1) {
           MeshFaceBoundary<Dimension> mesh_face_boundary =
             getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+          MeshNodeBoundary<Dimension> mesh_node_boundary =
+            getMeshNodeBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
 
           const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
           MeshDataType& mesh_data     = MeshDataManager::instance().getMeshData(*mesh);
@@ -167,7 +177,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                                                                                                       mesh_face_boundary
                                                                                                         .faceList());
 
-          boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list});
+          boundary_condition_list.push_back(
+            DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
 
         } else {
           throw NotImplementedError("Dirichlet BC in 1d");
@@ -451,85 +462,100 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         return kappa_invBeta;
       }();
 
+      const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
+        NodeValue<TinyVector<Dimension>> compute_sum_Cjr;
+        parallel_for(
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            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);
+            compute_sum_Cjr[node_id]                  = zero;
+            for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+              const CellId cell_id = node_to_cell[i_cell];
+              const size_t i_node  = node_local_number_in_its_cell[i_cell];
+              compute_sum_Cjr[node_id] += Cjr(cell_id, i_node);
+            }
+          });
+        return compute_sum_Cjr;
+      }();
+
+      const NodeValue<const TinyVector<Dimension>> v = [&] {
+        NodeValue<TinyVector<Dimension>> compute_v;
+        parallel_for(
+          mesh->numberOfNodes(),
+          PUGS_LAMBDA(NodeId node_id) { compute_v[node_id] = node_kappar[node_id] * exterior_normal[node_id]; });
+        return compute_v;
+      }();
+
+      const NodeValuePerCell<const double> theta = [&] {
+        NodeValuePerCell<double> compute_theta;
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+            const auto& cell_nodes = cell_to_node_matrix[cell_id];
+            for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+              const NodeId node_id           = cell_nodes[i_node];
+              compute_theta(cell_id, i_node) = dot(inverse(node_betar[node_id]) * Cjr(cell_id, i_node), v[node_id]);
+            }
+          });
+        return compute_theta;
+      }();
+
+      const NodeValue<const double> sum_theta = [&] {
+        NodeValue<double> compute_sum_theta;
+        parallel_for(
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            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);
+            compute_sum_theta[node_id]                = 0;
+            for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+              const CellId cell_id = node_to_cell[i_cell];
+              const size_t i_node  = node_local_number_in_its_cell[i_cell];
+              compute_sum_theta[node_id] += theta(cell_id, i_node);
+            }
+          });
+        return compute_sum_theta;
+      }();
+
       const Array<int> non_zeros{mesh->numberOfCells()};
       non_zeros.fill(Dimension);   // Modif pour optimiser
       CRSMatrixDescriptor<double> S(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
 
       for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-        const auto& node_to_cell = node_to_cell_matrix[node_id];
-
+        const auto& node_to_cell                   = node_to_cell_matrix[node_id];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
         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 i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+            const CellId cell_id_j = node_to_cell[i_cell_j];
+            const size_t i_node_j  = node_local_number_in_its_cells[i_cell_j];
 
-            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 (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+              const CellId cell_id_k = node_to_cell[i_cell_k];
+              const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
-              S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), Cjr(J, node_id_j));
+              S(cell_id_j, cell_id_k) +=
+                dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
             }
           }
         } else if (not node_is_corner[node_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];
-
-                    S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), Cjr(J, node_id_j));
-                  }
-                }
-              }
+          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;
+
+          for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+            const CellId cell_id_j = node_to_cell[i_cell_j];
+            const size_t i_node_j  = node_local_number_in_its_cells[i_cell_j];
+
+            for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+              const CellId cell_id_k = node_to_cell[i_cell_k];
+              const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
+
+              S(cell_id_j, cell_id_k) +=
+                dot(kappar_invBetar[node_id] *
+                      (Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
+                    P * Cjr(cell_id_k, i_node_j));
             }
-          },
-          boundary_condition);
+          }
+        }
       }
 
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
@@ -545,38 +571,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         b[cell_id] = fj[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, 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_id] += value_list[i_face];
-              }
-            } 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_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{mesh->numberOfCells()};
       T = zero;
 
@@ -598,18 +592,17 @@ ScalarNodalSchemeHandler::solution() const
 }
 
 ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
-  const std::shared_ptr<const IDiscreteFunction>& alpha,
   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)
 {
-  const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
+  const std::shared_ptr i_mesh = getCommonMesh({cell_k, f});
   if (not i_mesh) {
     throw NormalError("primal discrete functions are not defined on the same mesh");
   }
 
-  checkDiscretizationType({alpha, cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
+  checkDiscretizationType({cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
 
   switch (i_mesh->dimension()) {
   case 1: {
@@ -624,7 +617,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
     m_scheme =
-      std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+      std::make_unique<ScalarNodalScheme<2>>(mesh,
                                              std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
                                              std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
                                              std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
@@ -639,7 +632,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
     m_scheme =
-      std::make_unique<ScalarNodalScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+      std::make_unique<ScalarNodalScheme<3>>(mesh,
                                              std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
                                              std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
                                              std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
diff --git a/src/scheme/ScalarNodalScheme.hpp b/src/scheme/ScalarNodalScheme.hpp
index 253beaf716e0fa5f54a51b1c615275c2aa954284..27ec1818efd0c17dc861aaf7b1c00413b5ec542d 100644
--- a/src/scheme/ScalarNodalScheme.hpp
+++ b/src/scheme/ScalarNodalScheme.hpp
@@ -40,15 +40,12 @@ class ScalarNodalSchemeHandler
 
   std::shared_ptr<const IDiscreteFunction> solution() const;
 
-  ScalarNodalSchemeHandler(
-    const std::shared_ptr<const IDiscreteFunction>& alpha,
-    const std::shared_ptr<const IDiscreteFunction>& nod_k_bound,
-    const std::shared_ptr<const IDiscreteFunction>& nod_k,
-    const std::shared_ptr<const IDiscreteFunction>& f,
-    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+  ScalarNodalSchemeHandler(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);
 
   ~ScalarNodalSchemeHandler();
 };
 
-
 #endif   // SCALAR_NODAL_SCHEME_HPP