diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 9c8e47dacb4f36cc2f002d7155093d9d780bc6dd..c7fe9279c6fa0237f07cf803ae78dae518141595 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/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) { - A(0, 0) = nl[face_id][0] * mes_l[face_id]; - A(0, 1) = nl[face_id][1] * mes_l[face_id]; + 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(1, 0) = nl[face_id][1] * mes_l[face_id]; + } else { + A(0, 0) = -nl[face_id][0] * 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(1, 1) = nl[face_id][1] * mes_l[face_id]; + 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]) { - compute_node_boundary_values[node_id] += - value_list[i_face] * std::abs(dot(nl[face_id], exterior_normal[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] * 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];