From e6906815f84b5f2721cfbb3a255ebe3ea6e39e4f Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Mon, 30 May 2022 09:56:56 +0200
Subject: [PATCH] Add Beta_r at the corners for Dirichlet BC.

---
 src/scheme/ScalarNodalScheme.cpp | 109 ++++++++++++++++++++++++++++---
 1 file changed, 100 insertions(+), 9 deletions(-)

diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 8246e2121..c349f4b5e 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -420,10 +420,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             kappa[node_id]           = zero;
             const auto& node_to_cell = node_to_cell_matrix[node_id];
             double weight            = 0;
-            for (size_t j = 0; j < node_to_cell.size(); ++j) {
-              const CellId J      = node_to_cell[j];
-              double local_weight = 1. / l2Norm(xr[node_id] - xj[J]);
-              kappa[node_id] += local_weight * cell_kappaj[J];
+            for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+              const CellId cell_id = node_to_cell[i_cell];
+              double local_weight  = 1. / l2Norm(xr[node_id] - xj[cell_id]);
+              kappa[node_id] += local_weight * cell_kappaj[cell_id];
               weight += local_weight;
             }
             kappa[node_id] = 1. / weight * kappa[node_id];
@@ -438,10 +438,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             kappa[node_id]           = zero;
             const auto& node_to_cell = node_to_cell_matrix[node_id];
             double weight            = 0;
-            for (size_t j = 0; j < node_to_cell.size(); ++j) {
-              const CellId J      = node_to_cell[j];
-              double local_weight = 1. / l2Norm(xr[node_id] - xj[J]);
-              kappa[node_id] += local_weight * cell_kappajb[J];
+            for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+              const CellId cell_id = node_to_cell[i_cell];
+              double local_weight  = 1. / l2Norm(xr[node_id] - xj[cell_id]);
+              kappa[node_id] += local_weight * cell_kappajb[cell_id];
               weight += local_weight;
             }
             kappa[node_id] = 1. / weight * kappa[node_id];
@@ -487,6 +487,71 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         return kappa_invBeta;
       }();
 
+      const NodeValue<const TinyMatrix<Dimension>> corner_betar = [&] {
+        NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            if (node_is_corner[node_id]) {
+              const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+              const auto& node_to_cell                  = node_to_cell_matrix[node_id];
+
+              size_t i_cell             = node_to_cell[0];
+              const CellId cell_id      = node_to_cell[i_cell];
+              const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+
+              const auto& cell_nodes  = cell_to_node_matrix[cell_id];
+              const NodeId node_id_p1 = cell_to_node_matrix[cell_id][(i_node + 1) % cell_nodes.size()];
+              const NodeId node_id_p2 = cell_to_node_matrix[cell_id][(i_node + 2) % cell_nodes.size()];
+              const NodeId node_id_m1 = cell_to_node_matrix[cell_id][(i_node - 1) % cell_nodes.size()];
+
+              const TinyVector<Dimension> xj1 = 1. / 3. * (xr[node_id] + xr[node_id_m1] + xr[node_id_p2]);
+              const TinyVector<Dimension> xj2 = 1. / 3. * (xr[node_id] + xr[node_id_p2] + xr[node_id_p1]);
+
+              const TinyVector<Dimension> xrm1 = xr[node_id_m1];
+              const TinyVector<Dimension> xrp1 = xr[node_id_p1];
+              const TinyVector<Dimension> xrp2 = xr[node_id_p2];
+
+              TinyVector<Dimension> Cjr1;
+              TinyVector<Dimension> Cjr2;
+
+              if (Dimension == 2) {
+                Cjr1[0] = -0.5 * (xrm1[1] - xrp2[1]);
+                Cjr1[1] = 0.5 * (xrm1[0] - xrp2[0]);
+                Cjr2[0] = -0.5 * (xrp2[1] - xrp1[1]);
+                Cjr2[1] = 0.5 * (xrp2[0] - xrp1[0]);
+              } else {
+                std::cout << "Dimension " << Dimension << " not implemented. \n";
+                throw NotImplementedError("The scheme is not implemented in this dimension.");
+              }
+
+              beta[node_id] = tensorProduct(Cjr1, (xr[node_id] - xj1)) + tensorProduct(Cjr2, (xr[node_id] - xj2));
+            }
+          });
+        return beta;
+      }();
+
+      const NodeValue<const TinyMatrix<Dimension>> corner_kappar_invBetar = [&] {
+        NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            if (node_is_corner[node_id]) {
+              kappa_invBeta[node_id] = node_kappar[node_id] * inverse(corner_betar[node_id]);
+            }
+          });
+        return kappa_invBeta;
+      }();
+
+      const NodeValue<const TinyMatrix<Dimension>> corner_kapparb_invBetar = [&] {
+        NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            if (node_is_corner[node_id]) {
+              kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(corner_betar[node_id]);
+            }
+          });
+        return kappa_invBeta;
+      }();
+
       const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
         NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()};
         parallel_for(
@@ -564,7 +629,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 ((node_is_neumann[node_id]) && (not node_is_corner[node_id])) {
+        } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id]) && (not node_is_dirichlet[node_id])) {
           const TinyMatrix<Dimension> I   = identity;
           const TinyVector<Dimension> n   = exterior_normal[node_id];
           const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
@@ -584,6 +649,32 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                     P * Cjr(cell_id_k, i_node_j));
             }
           }
+        } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) {
+          for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+            const CellId cell_id_j = node_to_cell[i_cell_j];
+            const size_t i_node_j  = node_local_number_in_its_cells[i_cell_j];
+
+            for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+              const CellId cell_id_k = node_to_cell[i_cell_k];
+              const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
+
+              S(cell_id_j, cell_id_k) +=
+                dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+            }
+          }
+        } else if (node_is_dirichlet[node_id] && node_is_corner[node_id]) {
+          for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+            const CellId cell_id_j = node_to_cell[i_cell_j];
+            const size_t i_node_j  = node_local_number_in_its_cells[i_cell_j];
+
+            for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+              const CellId cell_id_k = node_to_cell[i_cell_k];
+              const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
+
+              S(cell_id_j, cell_id_k) +=
+                dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+            }
+          }
         }
       }
 
-- 
GitLab