Change the construction of the boundary values at the angles

src/scheme/ScalarNodalScheme.cpp

@@ -489,16 +489,30 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
         if (node_is_angle[node_id]) {
           const auto& node_to_face = node_to_face_matrix[node_id];
-          TinyMatrix<Dimension> A  = zero;
+          TinyMatrix<Dimension> A;
+          A(0, 0) = 1000;
+          A(0, 1) = 1000;
+          A(1, 0) = 1000;
+          A(1, 1) = 1000;
           for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
             const FaceId face_id = node_to_face[i_face];
             if (is_boundary_face[face_id]) {
-              if (A(0, 0) == 0) {
+              if (A(0, 0) == 1000) {
+                if (dot(nl[face_id], exterior_normal[node_id]) >= 0) {
                   A(0, 0) = nl[face_id][0] * mes_l[face_id];
-                A(0, 1) = nl[face_id][1] * mes_l[face_id];
+                  A(1, 0) = nl[face_id][1] * mes_l[face_id];
                 } else {
-                A(1, 0) = nl[face_id][0] * mes_l[face_id];
+                  A(0, 0) = -nl[face_id][0] * mes_l[face_id];
+                  A(1, 0) = -nl[face_id][1] * mes_l[face_id];
+                }
+              } else {
+                if (dot(nl[face_id], exterior_normal[node_id]) >= 0) {
+                  A(0, 1) = nl[face_id][0] * mes_l[face_id];
                   A(1, 1) = nl[face_id][1] * mes_l[face_id];
+                } else {
+                  A(0, 1) = -nl[face_id][0] * mes_l[face_id];
+                  A(1, 1) = -nl[face_id][1] * mes_l[face_id];
+                }
               }
             }
           }
@@ -528,9 +542,14 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                   const NodeId node_id = face_nodes[i_node];
                   if (not node_is_dirichlet[node_id]) {
                     if (node_is_angle[node_id]) {
+                      const auto& node_to_face = node_to_face_matrix[node_id];
+                      for (size_t i_face_node = 0; i_face_node < node_to_face.size(); ++i_face_node) {
+                        const FaceId face_node_id = node_to_face[i_face_node];
+                        if (face_id == face_node_id) {
                           compute_node_boundary_values[node_id] +=
-                        value_list[i_face] * std::abs(dot(nl[face_id], exterior_normal[node_id]));
-
+                            value_list[i_face] * mes_l[face_id] * normal_base[node_id][i_face_node];
+                        }
+                      }
                     } else {
                       compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
                       sum_mes_l[node_id] += mes_l[face_id];