From 85547a23cea42228ac6da0beae34dc9fd8bb9d0c Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Fri, 26 Aug 2022 10:19:05 +0200
Subject: [PATCH] Fix matrix indices for Dirichlet BC

---
 src/scheme/ScalarHybridScheme.cpp | 67 +++++++++----------------------
 1 file changed, 20 insertions(+), 47 deletions(-)

diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index 61285fab5..ecf6d84ed 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -1085,16 +1085,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
                 S2(cell_id_j, cell_id_k) -=
                   lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
                   dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
-                // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //       normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) !=
-                //     nl[face_id]) {
-                //   std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n"
-                //             << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //                  normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face)
-                //             << " = " << nl[face_id] << "\n"
-                //             << "\n";
-                // }
-                // std::cout << S1(cell_id_j, cell_id_k) << "\n";
               }
             }
           }
@@ -1146,23 +1136,15 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
               const FaceId face_id     = cell_to_face[i_face];
               const auto& face_to_node = face_to_node_matrix[face_id];
-
-              if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
-                S1(cell_id_j, cell_id_k) +=
-                  lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-                  dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
-                S2(cell_id_j, cell_id_k) -=
-                  lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-                  dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
-                // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //       normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) !=
-                //     nl[face_id]) {
-                //   std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n"
-                //             << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //                  normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face)
-                //             << " = " << nl[face_id] << "\n"
-                //             << "\n";
-                // }
+              if (face_is_dirichlet[face_id]) {
+                if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+                  S1(cell_id_j, cell_id_k) +=
+                    lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+                    dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
+                  S2(cell_id_j, cell_id_k) -=
+                    lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+                    dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
+                }
               }
             }
           }
@@ -1188,25 +1170,16 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
               const FaceId face_id     = cell_to_face[i_face];
               const auto& face_to_node = face_to_node_matrix[face_id];
+              if (face_is_dirichlet[face_id]) {
+                if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+                  S1(cell_id_j, cell_id_k) +=
+                    lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+                    dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
 
-              if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
-                S1(cell_id_j, cell_id_k) +=
-                  lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-                  dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
-
-                S2(cell_id_j, cell_id_k) -=
-                  lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-                  dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
-
-                // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //       normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) !=
-                //     nl[face_id]) {
-                //   std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n"
-                //             << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //                  normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face)
-                //             << " = " << nl[face_id] << "\n"
-                //             << "\n";
-                // }
+                  S2(cell_id_j, cell_id_k) -=
+                    lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+                    dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
+                }
               }
             }
           }
@@ -1340,7 +1313,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
               const FaceId face_id     = cell_to_face[i_face];
               const auto& face_to_node = face_to_node_matrix[face_id];
-              if (is_boundary_face[face_id]) {
+              if (face_is_dirichlet[face_id]) {
                 if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
                   b[cell_id_j] += lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
                                   node_boundary_values[node_id] *
@@ -1376,7 +1349,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
               const FaceId face_id     = cell_to_face[i_face];
               const auto& face_to_node = face_to_node_matrix[face_id];
-              if (is_boundary_face[face_id]) {
+              if (face_is_dirichlet[face_id]) {
                 if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
                   b[cell_id_j] +=
                     lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-- 
GitLab