Change of the angles calculation for Neumann BC

src/scheme/ScalarNodalScheme.cpp

@@ -313,6 +313,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
 
     const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
+    const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
 
     const auto& node_to_face_matrix               = mesh->connectivity().nodeToFaceMatrix();
     const auto& face_to_node_matrix               = mesh->connectivity().faceToNodeMatrix();
@@ -320,6 +321,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     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 auto& nl                                = mesh_data.nl();
 
     const NodeValue<const bool> node_is_neumann = [&] {
       NodeValue<bool> compute_node_is_neumann{mesh->connectivity()};
@@ -409,6 +411,35 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       return compute_node_is_corner;
     }();
 
+    const NodeValue<const bool> node_is_angle = [&] {
+      NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
+      compute_node_is_angle.fill(false);
+      parallel_for(
+        mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+          if (is_boundary_node[node_id]) {
+            const auto& node_to_face = node_to_face_matrix[node_id];
+
+            TinyVector<Dimension> n1 = zero;
+            TinyVector<Dimension> n2 = zero;
+
+            for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+              FaceId face_id = node_to_face[i_face];
+              if (is_boundary_face[face_id]) {
+                if (l2Norm(n1) == 0) {
+                  n1 = nl[face_id];
+                } else {
+                  n2 = nl[face_id];
+                }
+              }
+              if (n1 != n2) {
+                compute_node_is_angle[node_id] = true;
+              }
+            }
+          }
+        });
+      return compute_node_is_angle;
+    }();
+
     const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] {
       NodeValue<TinyVector<Dimension>> compute_exterior_normal{mesh->connectivity()};
       compute_exterior_normal.fill(zero);
@@ -429,6 +460,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       return compute_exterior_normal;
     }();
 
+    // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+    //  std::cout << node_id << "  " << exterior_normal[node_id] << "\n";
+    //};
+
     const FaceValue<const double> mes_l = [&] {
       if constexpr (Dimension == 1) {
         FaceValue<double> compute_mes_l{mesh->connectivity()};
@@ -459,11 +494,16 @@ 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 (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]));
+                    } else {
                       compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
                       sum_mes_l[node_id] += 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();
@@ -487,7 +527,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
           boundary_condition);
       }
       for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-        if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id])) {
+        if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
           compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
         } else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
           compute_node_boundary_values[node_id] /= sum_mes_l[node_id];