From b7f46037a75faa5692e193ff50f680f165301642 Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Tue, 23 Aug 2022 15:34:41 +0200
Subject: [PATCH] Fix sign of the matrix coefficients

---
 src/scheme/ScalarHybridScheme.cpp | 63 +++++++++++++++++--------------
 1 file changed, 34 insertions(+), 29 deletions(-)

diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index 10d41c948..55c38578f 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -694,6 +694,8 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
           }
 
           compute_normal_base(cell_id, i_face) = inverse(A) * (face_kappal[face_id] * face_normal);
+          // std::cout << "alpha_l = " << compute_normal_base(cell_id, i_face)[0]
+          //           << "; delta_l = " << compute_normal_base(cell_id, i_face)[1] << "\n";
         }
       }
 
@@ -1076,6 +1078,8 @@ 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_node_j));
+                // std::cout << "alpha_l = " << normal_face_base(cell_id_j, i_face)[0] << "\n";
+                // std::cout << S1(cell_id_j, cell_id_k) << "\n";
               }
             }
           }
@@ -1190,14 +1194,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             const CellId cell_id_k = face_to_cell[i_cell_k];
 
             if (cell_id_j == cell_id_k) {
-              S1(cell_id_j, cell_id_k) -=
+              S1(cell_id_j, cell_id_k) +=
                 lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
-              S2(cell_id_j, cell_id_k) +=
+              S2(cell_id_j, cell_id_k) -=
                 lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
             } else {
-              S1(cell_id_j, cell_id_k) +=
+              S1(cell_id_j, cell_id_k) -=
                 lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
-              S2(cell_id_j, cell_id_k) -=
+              S2(cell_id_j, cell_id_k) +=
                 lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
             }
           }
@@ -1214,6 +1218,8 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
     };
 
+    // std::cout << "S1 = " << S1 << "\n";
+
     for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
       S1(cell_id, cell_id) += Vj[cell_id];
       S2(cell_id, cell_id) += Vj[cell_id];
@@ -1229,13 +1235,11 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       Ph[cell_id] = Pj[cell_id];
     };
 
-    std::cout << "P = " << Ph << "\n";
+    // std::cout << "P = " << Ph << "\n";
 
     CRSMatrix A1{S1.getCRSMatrix()};
     CRSMatrix A2{S2.getCRSMatrix()};
 
-    // std::cout << A1 << "\n";
-
     Vector<double> b{mesh->numberOfCells()};
     b = zero;
     for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
@@ -1304,15 +1308,18 @@ 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) {
-                b[cell_id_j] += (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-                                face_boundary_values[face_id] *
-                                dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
-                b_prev[cell_id_j] +=
-                  (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-                  face_boundary_prev_values[face_id] *
-                  dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+              if (is_boundary_face[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] *
+                    dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+                  // std::cout << "xl - xj2 = " << xl[face_id] -
+                  b_prev[cell_id_j] +=
+                    lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+                    node_boundary_prev_values[node_id] *
+                    dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+                }
               }
             }
           }
@@ -1338,16 +1345,17 @@ 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) {
-                b[cell_id_j] +=
-                  (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-                  face_boundary_values[face_id] *
-                  dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
-                b_prev[cell_id_j] +=
-                  (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
-                  face_boundary_prev_values[face_id] *
-                  dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+              if (is_boundary_face[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] *
+                    dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+                  b_prev[cell_id_j] +=
+                    lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+                    node_boundary_prev_values[node_id] *
+                    dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+                }
               }
             }
           }
@@ -1390,9 +1398,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         B[cell_id] = dt * ((1. - cn_coeff / 2.) * b[cell_id] + cn_coeff / 2. * b_prev[cell_id]) + A2P[cell_id];
       });
 
-    std::cout << "B = " << B << "\n"
-              << "A2P = " << A2P << "\n";
-
     NodeValue<TinyVector<Dimension>> ur{mesh->connectivity()};
     ur.fill(zero);
     for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-- 
GitLab