Fix angle detection + calculus of ur

src/scheme/ScalarNodalScheme.cpp

@@ -414,6 +414,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     const NodeValue<const bool> node_is_angle = [&] {
       NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
       compute_node_is_angle.fill(false);
+      // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
       parallel_for(
         mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
           if (is_boundary_node[node_id]) {
@@ -431,10 +432,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                   n2 = nl[face_id];
                 }
               }
-              if (n1 != n2) {
-                compute_node_is_angle[node_id] = true;
             }
+            if (l2Norm(n1 - n2) > (1.E-15) and l2Norm(n1 + n2) > (1.E-15) and not node_is_corner[node_id]) {
+              compute_node_is_angle[node_id] = true;
             }
+            // std::cout << node_id << "  " << n1 << "  " << n2 << "  " << compute_node_is_angle[node_id] << "\n";
           }
         });
       return compute_node_is_angle;
@@ -461,7 +463,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     }();
 
     // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-    //  std::cout << node_id << "  " << exterior_normal[node_id] << "\n";
+    //  std::cout << node_id << "  " << exterior_normal[node_id] << "  " << node_is_angle[node_id] << "\n";
     //};
 
     const FaceValue<const double> mes_l = [&] {
@@ -920,8 +922,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               FaceId face_id = node_to_face[i_face];
               sum_mes_l += mes_l[face_id];
             }
-            b[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_values[node_id] * sum_mes_l;
-            b_prev[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_prev_values[node_id] * sum_mes_l;
+            b[cell_id] += 0.5 / Vj[cell_id] * node_boundary_values[node_id] * sum_mes_l;
+            b_prev[cell_id] += 0.5 / Vj[cell_id] * node_boundary_prev_values[node_id] * sum_mes_l;
           }
         } else if (node_is_dirichlet[node_id]) {
           if (not node_is_corner[node_id]) {
@@ -980,6 +982,34 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
       // std::cout << "g = " << node_boundary_values << "\n";
 
+      NodeValue<TinyVector<Dimension>> ur{mesh->connectivity()};
+      ur.fill(zero);
+      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_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+        if (not is_boundary_node[node_id]) {
+          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];
+            ur[node_id] += Pj[cell_id] * Cjr(cell_id, i_node);
+          }
+          ur[node_id] = -inverse(node_betar[node_id]) * ur[node_id];
+          // std::cout << node_id << "; ur = " << ur[node_id] << "\n";
+        } else if (is_boundary_node[node_id] and node_is_dirichlet[node_id]) {
+          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];
+            ur[node_id] += (node_boundary_values[node_id] - Pj[cell_id]) * Cjr(cell_id, i_node);
+          }
+          if (not node_is_corner[node_id]) {
+            ur[node_id] = inverse(node_betar[node_id]) * ur[node_id];
+          } else {
+            ur[node_id] = inverse(corner_betar[node_id]) * ur[node_id];
+          }
+        }
+        // std::cout << node_id << "; ur = " << ur[node_id] << "\n";
+      }
+
       LinearSolver solver;
       solver.solveLocalSystem(A1, T, B);