diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 9c442d0329900c6bf49f52ca050db4dc6cef511d..fa374bdfcb84a7680d2944af88aa4cc161f69c55 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -370,14 +370,14 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
                 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) {
+                  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];
                   }
                 }
               }
               for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-                if (not node_is_dirichlet[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];
                 }
               }
@@ -391,13 +391,19 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
                 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];
+                  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] * mes_l[face_id];
+                    compute_node_boundary_values[node_id] = value_list[i_face];
                   }
                 }
               }
+              for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+                if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
+                  compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
+                }
+              }
             }
           },
           boundary_condition);