From 2f338e16cf1f51f5af5c96aa57e979e88e808c2d Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Fri, 20 May 2022 15:37:55 +0200
Subject: [PATCH] Correction k_r on the border + construction second member for
 Neumann BC

---
 src/scheme/ScalarNodalScheme.cpp | 58 ++++++++++++++++++++++++++------
 1 file changed, 47 insertions(+), 11 deletions(-)

diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index fa374bdfc..e6f63ccaf 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/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,18 +593,45 @@ 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];
-      }
+      };
 
-      Vector<double> T{mesh->numberOfCells()};
-      T = 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_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];
+            }
+          }
+        };
 
-      LinearSolver solver;
-      solver.solveLocalSystem(A, T, b);
+        Vector<double> T{mesh->numberOfCells()};
+        T = zero;
 
-      m_solution     = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
-      auto& solution = *m_solution;
-      parallel_for(
-        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; });
+        LinearSolver solver;
+        solver.solveLocalSystem(A, T, b);
+
+        m_solution     = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+        auto& solution = *m_solution;
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; });
+      }
     }
   }
 };
-- 
GitLab