Correction k_r on the border + construction second member for Neumann BC

src/scheme/ScalarNodalScheme.cpp

@@ -238,7 +238,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
     const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
 
-    // const auto& node_to_face_matrix               = mesh->connectivity().nodeToFaceMatrix();
+    const auto& node_to_face_matrix               = mesh->connectivity().nodeToFaceMatrix();
     const auto& face_to_node_matrix               = mesh->connectivity().faceToNodeMatrix();
     const auto& cell_to_node_matrix               = mesh->connectivity().cellToNodeMatrix();
     const auto& node_local_numbers_in_their_cells = mesh->connectivity().nodeLocalNumbersInTheirCells();
@@ -476,6 +476,15 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         return kappa_invBeta;
       }();
 
+      const NodeValue<const TinyMatrix<Dimension>> kapparb_invBetar = [&] {
+        NodeValue<TinyMatrix<Dimension>> kappa_invBeta;
+        parallel_for(
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(node_betar[node_id]);
+          });
+        return kappa_invBeta;
+      }();
+
       const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
         NodeValue<TinyVector<Dimension>> compute_sum_Cjr;
         parallel_for(
@@ -496,7 +505,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         NodeValue<TinyVector<Dimension>> compute_v;
         parallel_for(
           mesh->numberOfNodes(),
-          PUGS_LAMBDA(NodeId node_id) { compute_v[node_id] = node_kappar[node_id] * exterior_normal[node_id]; });
+          PUGS_LAMBDA(NodeId node_id) { compute_v[node_id] = node_kapparb[node_id] * exterior_normal[node_id]; });
         return compute_v;
       }();
 
@@ -584,7 +593,33 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       b = zero;
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
         b[cell_id] = fj[cell_id] * Vj[cell_id];
+      };
+
+      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_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+        const TinyMatrix<Dimension> I              = identity;
+        const TinyVector<Dimension> n              = exterior_normal[node_id];
+        const TinyMatrix<Dimension> nxn            = tensorProduct(n, n);
+        const TinyMatrix<Dimension> P              = I - nxn;
+        if ((node_is_neumann[node_id]) && (not node_is_dirichlet[node_id])) {
+          if ((is_boundary_node[node_id]) and (not node_is_corner[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_cells[i_cell];
+
+              b[cell_id] += node_boundary_values[node_id] *
+                            (1. / (sum_theta[node_id] * l2Norm(node_kapparb[node_id] * n)) *
+                               dot(kapparb_invBetar[node_id] * sum_Cjr[node_id], P * Cjr(cell_id, i_node)) -
+                             dot(Cjr(cell_id, i_node), n));
+            }
+          } else if (node_is_corner[node_id]) {
+            const auto& node_to_face = node_to_face_matrix[node_id];
+            for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+              FaceId face_id = node_to_face[i_face];
+            }
           }
+        };
 
         Vector<double> T{mesh->numberOfCells()};
         T = zero;
@@ -598,6 +633,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
           mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; });
       }
     }
+  }
 };
 
 std::shared_ptr<const IDiscreteFunction>