From 1a66c65f19538cb34d729b5dca2ab2b79e6790ed Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Fri, 20 May 2022 14:44:49 +0200
Subject: [PATCH] Correction node boundary values for Neumann BC

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

diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index d1f22b144..9c442d032 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -353,7 +353,9 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
     const NodeValue<const double> node_boundary_values = [&] {
       NodeValue<double> compute_node_boundary_values;
+      NodeValue<double> sum_mes_l;
       compute_node_boundary_values.fill(0);
+      sum_mes_l.fill(0);
       for (const auto& boundary_condition : boundary_condition_list) {
         std::visit(
           [&](auto&& bc) {
@@ -369,10 +371,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 (node_is_dirichlet[node_id] == false) {
-                    compute_node_boundary_values[node_id] += 0.5 * value_list[i_face] * mes_l[face_id];
+                    compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
+                    sum_mes_l[node_id] += mes_l[face_id];
                   }
                 }
               }
+              for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+                if (not node_is_dirichlet[node_id]) {
+                  compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
+                }
+              }
             } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
               const auto& face_list  = bc.faceList();
               const auto& value_list = bc.valueList();
@@ -535,7 +543,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                 dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
             }
           }
-        } else if (not node_is_corner[node_id]) {
+        } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id])) {
           const TinyMatrix<Dimension> I   = identity;
           const TinyVector<Dimension> n   = exterior_normal[node_id];
           const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
@@ -565,6 +573,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       CellValue<const double> fj = f->cellValues();
 
       CRSMatrix A{S.getCRSMatrix()};
+
       Vector<double> b{mesh->numberOfCells()};
       b = zero;
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-- 
GitLab