From 89a5ada3b68a5e60e88fbfd8325c976f73745a5c Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Wed, 14 Sep 2022 16:38:53 +0200
Subject: [PATCH] Change of the base xj1_xj2_orth for 3D

---
 src/scheme/ScalarHybridScheme.cpp | 433 ++++++++++++++++++++----------
 1 file changed, 291 insertions(+), 142 deletions(-)

diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index f0c4414e9..84fed063f 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -332,7 +332,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         throw NormalError(error_msg.str());
       }
     }
-
+    // std::cout << "Test 1.2 \n";
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
 
     const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
@@ -346,6 +346,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
     const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
     const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
 
+    // std::cout << "Test 1.3 \n";
     const auto& face_to_cell_matrix               = mesh->connectivity().faceToCellMatrix();
     const auto& node_to_face_matrix               = mesh->connectivity().nodeToFaceMatrix();
     const auto& face_to_node_matrix               = mesh->connectivity().faceToNodeMatrix();
@@ -358,6 +359,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
     const auto& node_to_cell_matrix               = mesh->connectivity().nodeToCellMatrix();
     const auto& nl                                = mesh_data.nl();
 
+    // std::cout << "Test 1.4 \n";
     const NodeValue<const bool> node_is_neumann = [&] {
       NodeValue<bool> compute_node_is_neumann{mesh->connectivity()};
       compute_node_is_neumann.fill(false);
@@ -381,6 +383,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_node_is_neumann;
     }();
+    // std::cout << "Test 1.5 \n";
 
     const FaceValue<const bool> face_is_neumann = [&] {
       FaceValue<bool> compute_face_is_neumann{mesh->connectivity()};
@@ -401,6 +404,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_face_is_neumann;
     }();
+    // std::cout << "Test 1.6 \n";
 
     const NodeValue<const bool> node_is_dirichlet = [&] {
       NodeValue<bool> compute_node_is_dirichlet{mesh->connectivity()};
@@ -423,6 +427,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_node_is_dirichlet;
     }();
+    // std::cout << "Test 1.7 \n";
 
     const FaceValue<const bool> face_is_dirichlet = [&] {
       FaceValue<bool> compute_face_is_dirichlet{mesh->connectivity()};
@@ -443,6 +448,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_face_is_dirichlet;
     }();
+    // std::cout << "Test 1.8 \n";
 
     const NodeValue<const bool> node_is_corner = [&] {
       NodeValue<bool> compute_node_is_corner{mesh->connectivity()};
@@ -456,28 +462,35 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
               const unsigned int i_node = node_local_number_in_its_cell[i_cell];
               const CellId cell_id      = node_to_cell[i_cell];
               const auto& cell_nodes    = cell_to_node_matrix[cell_id];
-              NodeId i_node_prev;
-              NodeId i_node_next;
-              if (i_node == 0) {
-                i_node_prev = cell_nodes.size() - 1;
-              } else {
-                i_node_prev = i_node - 1;
-              }
-              if (i_node == cell_nodes.size() - 1) {
-                i_node_next = 0;
-              } else {
-                i_node_next = i_node + 1;
-              }
-              const NodeId prev_node_id = cell_to_node_matrix[cell_id][i_node_prev];
-              const NodeId next_node_id = cell_to_node_matrix[cell_id][i_node_next];
-              if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
+              if (Dimension == 2) {
+                NodeId i_node_prev;
+                NodeId i_node_next;
+                if (i_node == 0) {
+                  i_node_prev = cell_nodes.size() - 1;
+                } else {
+                  i_node_prev = i_node - 1;
+                }
+                if (i_node == cell_nodes.size() - 1) {
+                  i_node_next = 0;
+                } else {
+                  i_node_next = i_node + 1;
+                }
+                const NodeId prev_node_id = cell_to_node_matrix[cell_id][i_node_prev];
+                const NodeId next_node_id = cell_to_node_matrix[cell_id][i_node_next];
+                if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
+                  compute_node_is_corner[node_id] = true;
+                }
+              } else if (Dimension == 3) {
                 compute_node_is_corner[node_id] = true;
+              } else {
+                throw NotImplementedError("The scheme is not implemented in this dimension.");
               }
             }
           }
         });
       return compute_node_is_corner;
     }();
+    ////std::cout << "Test 1.9 \n";
 
     const NodeValue<const bool> node_is_angle = [&] {
       NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
@@ -509,6 +522,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return compute_node_is_angle;
     }();
+    // std::cout << "Test 1.10 \n";
 
     const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] {
       NodeValue<TinyVector<Dimension>> compute_exterior_normal{mesh->connectivity()};
@@ -529,6 +543,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_exterior_normal;
     }();
+    // std::cout << "Test 1.11 \n";
 
     // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
     //  std::cout << node_id << "  " << exterior_normal[node_id] << "  " << node_is_angle[node_id] << "\n";
@@ -543,6 +558,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         return mesh_data.ll();
       }
     }();
+    // std::cout << "Test 1.12 \n";
 
     const NodeValue<const TinyVector<Dimension>> normal_angle_base = [&] {
       NodeValue<TinyVector<Dimension>> compute_normal_base{mesh->connectivity()};
@@ -582,6 +598,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_normal_base;
     }();
+    // std::cout << "Test 1.13 \n";
 
     const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
       NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
@@ -601,6 +618,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return kappa;
     }();
+    // std::cout << "Test 1.14 \n";
 
     const FaceValue<const TinyMatrix<Dimension>> face_kappal = [&] {
       FaceValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
@@ -608,10 +626,15 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       parallel_for(
         mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
           const auto& face_to_node = face_to_node_matrix[face_id];
-          kappa[face_id]           = 0.5 * (node_kappar[face_to_node[0]] + node_kappar[face_to_node[1]]);
+          const size_t n           = face_to_node.size();
+          for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) {
+            kappa[face_id] += node_kappar[face_to_node[i_node]];
+          }
+          kappa[face_id] = 1. / n * kappa[face_id];
         });
       return kappa;
     }();
+    // std::cout << "Test 1.15 \n";
 
     const FaceValuePerCell<const TinyVector<Dimension>> xj1_xj2 = [&] {
       FaceValuePerCell<TinyVector<Dimension>> compute_xj1_xj2{mesh->connectivity()};
@@ -647,13 +670,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_xj1_xj2;
     }();
+    // std::cout << "Test 1.16 \n";
 
     const FaceValuePerCell<const TinyVector<Dimension>> xj1_xj2_orth = [&] {
       FaceValuePerCell<TinyVector<Dimension>> compute_xj1_xj2_orth{mesh->connectivity()};
       compute_xj1_xj2_orth.fill(zero);
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
         const auto& cell_to_face = cell_to_face_matrix[cell_id];
-
+        // std::cout << "Test 1.16.1 \n";
         for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
           const FaceId face_id = cell_to_face[i_face];
 
@@ -663,44 +687,133 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_xj1_xj2_orth;
     }();
+    // std::cout << "Test 1.17 \n";
+
+    const FaceValuePerCell<const TinyVector<Dimension>> w1 = [&] {
+      FaceValuePerCell<TinyVector<Dimension>> compute_w1{mesh->connectivity()};
+      compute_w1.fill(zero);
+      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+        const auto& cell_to_face = cell_to_face_matrix[cell_id];
+        for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+          const FaceId face_id = cell_to_face[i_face];
+          TinyVector<Dimension> v;
+          v[0] = xj1_xj2(cell_id, i_face)[0] + 1;
+          v[1] = xj1_xj2(cell_id, i_face)[1];
+          v[2] = 0;
+          if (pow(dot(v, xj1_xj2(cell_id, i_face)), 2) == pow(l2Norm(v) * l2Norm(xj1_xj2(cell_id, i_face)), 2)) {
+            v[2] = 1;
+          }
+          if (Dimension == 3) {
+            // compute_w1(cell_id, i_face) = crossProduct(v, xj1_xj2(cell_id, i_face));
+            compute_w1(cell_id, i_face)[0] = v[1] * xj1_xj2(cell_id, i_face)[2] - v[2] * xj1_xj2(cell_id, i_face)[1];
+            compute_w1(cell_id, i_face)[1] = v[2] * xj1_xj2(cell_id, i_face)[0] - v[0] * xj1_xj2(cell_id, i_face)[2];
+            compute_w1(cell_id, i_face)[2] = v[0] * xj1_xj2(cell_id, i_face)[1] - v[1] * xj1_xj2(cell_id, i_face)[0];
+          } else {
+            compute_w1(cell_id, i_face) = zero;
+          }
+        }
+      }
+      return compute_w1;
+    }();
+
+    const FaceValuePerCell<const TinyVector<Dimension>> w2 = [&] {
+      FaceValuePerCell<TinyVector<Dimension>> compute_w2{mesh->connectivity()};
+      compute_w2.fill(zero);
+      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+        const auto& cell_to_face = cell_to_face_matrix[cell_id];
+        for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+          const FaceId face_id = cell_to_face[i_face];
+
+          // compute_w2(cell_id, i_face) = crossProduct(xj1_xj2(cell_id, i_face), w1(cell_id, i_face));
+          if (Dimension == 3) {
+            compute_w2(cell_id, i_face)[0] = w1(cell_id, i_face)[1] * xj1_xj2(cell_id, i_face)[2] -
+                                             w1(cell_id, i_face)[2] * xj1_xj2(cell_id, i_face)[1];
+            compute_w2(cell_id, i_face)[1] = w1(cell_id, i_face)[2] * xj1_xj2(cell_id, i_face)[0] -
+                                             w1(cell_id, i_face)[0] * xj1_xj2(cell_id, i_face)[2];
+            compute_w2(cell_id, i_face)[2] = w1(cell_id, i_face)[0] * xj1_xj2(cell_id, i_face)[1] -
+                                             w1(cell_id, i_face)[1] * xj1_xj2(cell_id, i_face)[0];
+          } else {
+            compute_w2(cell_id, i_face) = zero;
+          }
+        }
+      }
+      return compute_w2;
+    }();
 
     const FaceValuePerCell<const TinyVector<Dimension>> normal_face_base = [&] {
       FaceValuePerCell<TinyVector<Dimension>> compute_normal_base{mesh->connectivity()};
       compute_normal_base.fill(zero);
+      // std::cout << "Test 1.17.01 \n";
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
         const auto& cell_to_face = cell_to_face_matrix[cell_id];
-
+        // std::cout << "Test 1.17.1 \n";
         for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
           const FaceId face_id = cell_to_face[i_face];
-          TinyMatrix<Dimension> A;
+          // TinyMatrix<Dimension> A;
           TinyVector<Dimension> face_normal;
 
-          A(0, 0) = xj1_xj2_orth(cell_id, i_face)[0];
-          A(0, 1) = xj1_xj2(cell_id, i_face)[0];
-          A(1, 0) = xj1_xj2_orth(cell_id, i_face)[1];
-          A(1, 1) = xj1_xj2(cell_id, i_face)[1];
-
+          // A(0, 0) = xj1_xj2_orth(cell_id, i_face)[0];
+          // A(0, 1) = xj1_xj2(cell_id, i_face)[0];
+          // A(1, 0) = xj1_xj2_orth(cell_id, i_face)[1];
+          // A(1, 1) = xj1_xj2(cell_id, i_face)[1];
+          // std::cout << "Test 1.17.2 \n";
           if (dot(nl[face_id], xj1_xj2(cell_id, i_face)) > 0) {
             face_normal = nl[face_id];
           } else {
             face_normal = -nl[face_id];
           }
 
-          // compute_normal_base(cell_id, i_face)[0] =
-          //   dot(face_kappal[face_id] * face_normal, 1. / l2Norm(xj1_xj2(cell_id, i_face)) * xj1_xj2(cell_id,
-          //   i_face));
-          // compute_normal_base(cell_id, i_face)[1] =
-          //   dot(face_kappal[face_id] * face_normal,
-          //       1. / l2Norm(xj1_xj2_orth(cell_id, i_face)) * xj1_xj2_orth(cell_id, i_face));
-          compute_normal_base(cell_id, i_face) = inverse(A) * (face_kappal[face_id] * face_normal);
-
-          // std::cout << cell_id << " " << face_id << "; alpha_l = " << compute_normal_base(cell_id, i_face)[0]
-          //           << "; delta_l = " << compute_normal_base(cell_id, i_face)[1] << "\n";
+          // compute_normal_base(cell_id, i_face) = inverse(A) * (face_kappal[face_id] * face_normal);
+          // std::cout << "Test 1.17.3 \n";
+          // std::cout << face_kappal[face_id] << "\n";
+          // std::cout << xj1_xj2_orth(cell_id, i_face) << "\n";
+          // std::cout << face_normal << "\n";
+
+          if (Dimension == 2) {
+            compute_normal_base(cell_id, i_face)[1] =
+              dot(face_kappal[face_id] * face_normal, xj1_xj2_orth(cell_id, i_face)) /
+              pow(l2Norm(xj1_xj2_orth(cell_id, i_face)), 2);
+            compute_normal_base(cell_id, i_face)[0] =
+              dot(face_kappal[face_id] * face_normal, xj1_xj2(cell_id, i_face)) /
+              pow(l2Norm(xj1_xj2(cell_id, i_face)), 2);
+          } else if (Dimension == 3) {
+            compute_normal_base(cell_id, i_face)[1] =
+              dot(face_kappal[face_id] * face_normal, w1(cell_id, i_face)) / pow(l2Norm(w1(cell_id, i_face)), 2);
+            compute_normal_base(cell_id, i_face)[2] =
+              dot(face_kappal[face_id] * face_normal, w2(cell_id, i_face)) / pow(l2Norm(w2(cell_id, i_face)), 2);
+            compute_normal_base(cell_id, i_face)[0] =
+              dot(face_kappal[face_id] * face_normal, xj1_xj2(cell_id, i_face)) /
+              pow(l2Norm(xj1_xj2(cell_id, i_face)), 2);
+          }
+          // std::cout << "Test 1.17.4 \n";
+          //  std::cout << cell_id << " " << face_id << "; alpha_l = " << compute_normal_base(cell_id, i_face)[0]
+          //            << "; delta_l = " << compute_normal_base(cell_id, i_face)[1] << "\n";
         }
       }
 
       return compute_normal_base;
     }();
+    // std::cout << "Test 1.18 \n";
+
+    const FaceValuePerCell<const TinyVector<Dimension>> base_xj1_xj2_orth = [&] {
+      FaceValuePerCell<TinyVector<Dimension>> compute_w{mesh->connectivity()};
+      compute_w.fill(zero);
+      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+        const auto& cell_to_face = cell_to_face_matrix[cell_id];
+        for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+          const FaceId face_id = cell_to_face[i_face];
+          if (Dimension == 2) {
+            compute_w(cell_id, i_face) = normal_face_base(cell_id, i_face)[1] * xj1_xj2_orth(cell_id, i_face);
+          } else if (Dimension == 3) {
+            compute_w(cell_id, i_face) = normal_face_base(cell_id, i_face)[1] * w1(cell_id, i_face) +
+                                         normal_face_base(cell_id, i_face)[2] * w2(cell_id, i_face);
+          } else {
+            throw NotImplementedError("The scheme is not implemented in this dimension.");
+          }
+        }
+      }
+      return compute_w;
+    }();
 
     const FaceValue<const double> face_boundary_values = [&] {
       FaceValue<double> compute_face_boundary_values{mesh->connectivity()};
@@ -735,6 +848,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_face_boundary_values;
     }();
+    // std::cout << "Test 1.19 \n";
 
     const FaceValue<const double> face_boundary_prev_values = [&] {
       FaceValue<double> compute_face_boundary_values{mesh->connectivity()};
@@ -769,6 +883,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_face_boundary_values;
     }();
+    // std::cout << "Test 1.20 \n";
 
     const NodeValue<const double> node_boundary_values = [&] {
       NodeValue<double> compute_node_boundary_values{mesh->connectivity()};
@@ -827,6 +942,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_node_boundary_values;
     }();
+    // std::cout << "Test 1.21 \n";
 
     const NodeValue<const double> node_boundary_prev_values = [&] {
       NodeValue<double> compute_node_boundary_values{mesh->connectivity()};
@@ -883,6 +999,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
       }
       return compute_node_boundary_values;
     }();
+    // std::cout << "Test 1.22 \n";
 
     const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
       NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
@@ -900,6 +1017,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return beta;
     }();
+    // std::cout << "Test 1.23 \n";
 
     const NodeValue<const TinyMatrix<Dimension>> kappar_invBetar = [&] {
       NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
@@ -911,6 +1029,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return kappa_invBeta;
     }();
+    // std::cout << "Test 1.24 \n";
 
     const NodeValue<const TinyMatrix<Dimension>> corner_betar = [&] {
       NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
@@ -920,39 +1039,42 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             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             = 0;
-            const CellId cell_id      = node_to_cell[i_cell];
-            const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+            size_t i_cell        = 0;
+            const CellId cell_id = node_to_cell[i_cell];
+            if (Dimension == 2) {
+              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 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> 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];
+              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;
+              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]);
+              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]);
+              beta[node_id] = tensorProduct(Cjr1, (xr[node_id] - xj1)) + tensorProduct(Cjr2, (xr[node_id] - xj2));
             } else {
-              throw NotImplementedError("The scheme is not implemented in this dimension.");
-            }
+              const TinyMatrix<Dimension> I = identity;
 
-            beta[node_id] = tensorProduct(Cjr1, (xr[node_id] - xj1)) + tensorProduct(Cjr2, (xr[node_id] - xj2));
+              beta[node_id] = 0.125 * Vj[cell_id] * I;
+              // throw NotImplementedError("The scheme is not implemented in this dimension.");
+            }
           }
         });
       return beta;
     }();
+    // std::cout << "Test 1.25 \n";
 
     const NodeValue<const TinyMatrix<Dimension>> corner_kappar_invBetar = [&] {
       NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
@@ -964,6 +1086,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return kappa_invBeta;
     }();
+    // std::cout << "Test 1.26 \n";
 
     const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
       NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()};
@@ -981,6 +1104,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return compute_sum_Cjr;
     }();
+    // std::cout << "Test 1.27 \n";
 
     const NodeValue<const TinyVector<Dimension>> v = [&] {
       NodeValue<TinyVector<Dimension>> compute_v{mesh->connectivity()};
@@ -993,6 +1117,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return compute_v;
     }();
+    // std::cout << "Test 1.28 \n";
 
     const NodeValuePerCell<const double> theta = [&] {
       NodeValuePerCell<double> compute_theta{mesh->connectivity()};
@@ -1008,6 +1133,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return compute_theta;
     }();
+    // std::cout << "Test 1.29 \n";
 
     const NodeValue<const double> sum_theta = [&] {
       NodeValue<double> compute_sum_theta{mesh->connectivity()};
@@ -1027,7 +1153,9 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
         });
       return compute_sum_theta;
     }();
+    // std::cout << "Test 1.30 \n";
 
+    // std::cout << "Test 1.1 \n";
     const Array<int> non_zeros{mesh->numberOfCells()};
     non_zeros.fill(4);   // Modif pour optimiser
     CRSMatrixDescriptor<double> S1(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
@@ -1068,16 +1196,20 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
 
             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];
+              const FaceId face_id              = cell_to_face[i_face];
+              const auto& face_to_node          = face_to_node_matrix[face_id];
+              const size_t numberOfNodesPerFace = face_to_node.size();
               if (not is_boundary_face[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));
+                for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) {
+                  if (face_to_node[i_face_to_node] == node_id) {
+                    S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] /
+                                                numberOfNodesPerFace *
+                                                dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k),
+                                                    base_xj1_xj2_orth(cell_id_j, i_face));
+                    S2(cell_id_j, cell_id_k) -= lambda * dt * (cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace *
+                                                dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k),
+                                                    base_xj1_xj2_orth(cell_id_j, i_face));
+                  }
                 }
               }
             }
@@ -1111,22 +1243,25 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
 
             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];
+              const FaceId face_id              = cell_to_face[i_face];
+              const auto& face_to_node          = face_to_node_matrix[face_id];
+              const size_t numberOfNodesPerFace = face_to_node.size();
               if (not is_boundary_face[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) -
-                           theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
-                        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) -
-                           theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
-                        xj1_xj2_orth(cell_id_j, i_face));
+                for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) {
+                  if (face_to_node[i_face_to_node] == node_id) {
+                    S1(cell_id_j, cell_id_k) +=
+                      lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace *
+                      dot(inverse(node_betar[node_id]) *
+                            (Cjr(cell_id_k, i_node_k) -
+                             theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
+                          base_xj1_xj2_orth(cell_id_j, i_face));
+                    S2(cell_id_j, cell_id_k) -=
+                      lambda * dt * (cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace *
+                      dot(inverse(node_betar[node_id]) *
+                            (Cjr(cell_id_k, i_node_k) -
+                             theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
+                          base_xj1_xj2_orth(cell_id_j, i_face));
+                  }
                 }
               }
             }
@@ -1151,16 +1286,20 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
 
             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];
+              const FaceId face_id              = cell_to_face[i_face];
+              const auto& face_to_node          = face_to_node_matrix[face_id];
+              const size_t numberOfNodesPerFace = face_to_node.size();
               if (face_is_dirichlet[face_id] or not is_boundary_face[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));
+                for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) {
+                  if (face_to_node[i_face_to_node] == node_id) {
+                    S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] /
+                                                numberOfNodesPerFace *
+                                                dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k),
+                                                    base_xj1_xj2_orth(cell_id_j, i_face));
+                    S2(cell_id_j, cell_id_k) -= lambda * dt * (cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace *
+                                                dot(inverse(node_betar[node_id]) * Cjr(cell_id_k, i_node_k),
+                                                    base_xj1_xj2_orth(cell_id_j, i_face));
+                  }
                 }
               }
             }
@@ -1185,17 +1324,21 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
 
             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];
+              const FaceId face_id              = cell_to_face[i_face];
+              const auto& face_to_node          = face_to_node_matrix[face_id];
+              const size_t numberOfNodesPerFace = face_to_node.size();
               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));
-
-                  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));
+                for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) {
+                  if (face_to_node[i_face_to_node] == node_id) {
+                    S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] /
+                                                numberOfNodesPerFace *
+                                                dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k),
+                                                    base_xj1_xj2_orth(cell_id_j, i_face));
+
+                    S2(cell_id_j, cell_id_k) -= lambda * dt * (cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace *
+                                                dot(inverse(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k),
+                                                    base_xj1_xj2_orth(cell_id_j, i_face));
+                  }
                 }
               }
             }
@@ -1217,14 +1360,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
 
             if (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];
+                lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0];
               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];
+                lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0];
             } else {
               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];
+                lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0];
               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];
+                lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0];
             }
           }
         }
@@ -1233,9 +1376,9 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
           const CellId cell_id_j = face_to_cell[i_cell_j];
           const size_t i_face    = face_local_number_in_its_cells[i_cell_j];
           S1(cell_id_j, cell_id_j) +=
-            lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
+            lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0];
           S2(cell_id_j, cell_id_j) -=
-            lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
+            lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0];
         }
       }
     };
@@ -1299,19 +1442,22 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             const auto& cell_to_face = cell_to_face_matrix[cell_id];
 
             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];
+              const FaceId face_id              = cell_to_face[i_face];
+              const auto& face_to_node          = face_to_node_matrix[face_id];
+              const size_t numberOfNodesPerFace = face_to_node.size();
               if (not is_boundary_face[face_id]) {
-                if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
-                  b[cell_id] += lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id, i_face)[0] *
-                                node_boundary_values[node_id] /
-                                (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
-                                dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id, i_face));
-                  // std::cout << "xl - xj2 = " << xl[face_id] -
-                  b_prev[cell_id] +=
-                    lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id, i_face)[0] *
-                    node_boundary_prev_values[node_id] / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
-                    dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id, i_face));
+                for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) {
+                  if (face_to_node[i_face_to_node] == node_id) {
+                    b[cell_id] +=
+                      lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_values[node_id] /
+                      (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
+                      dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id, i_face));
+                    // std::cout << "xl - xj2 = " << xl[face_id] -
+                    b_prev[cell_id] +=
+                      lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_prev_values[node_id] /
+                      (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
+                      dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id, i_face));
+                  }
                 }
               }
             }
@@ -1349,18 +1495,20 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
 
             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];
+              const FaceId face_id              = cell_to_face[i_face];
+              const auto& face_to_node          = face_to_node_matrix[face_id];
+              const size_t numberOfNodesPerFace = face_to_node.size();
               if (face_is_dirichlet[face_id] or not 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_face));
-                  // 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_face));
+                for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) {
+                  if (face_to_node[i_face_to_node] == node_id) {
+                    b[cell_id_j] +=
+                      lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_values[node_id] *
+                      dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id_j, i_face));
+                    // std::cout << "xl - xj2 = " << xl[face_id] -
+                    b_prev[cell_id_j] +=
+                      lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_prev_values[node_id] *
+                      dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id_j, i_face));
+                  }
                 }
               }
             }
@@ -1385,18 +1533,19 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
 
             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];
+              const FaceId face_id              = cell_to_face[i_face];
+              const auto& face_to_node          = face_to_node_matrix[face_id];
+              const size_t numberOfNodesPerFace = face_to_node.size();
               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] *
-                    dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_face));
-                  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_face));
+                for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) {
+                  if (face_to_node[i_face_to_node] == node_id) {
+                    b[cell_id_j] +=
+                      lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_values[node_id] *
+                      dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id_j, i_face));
+                    b_prev[cell_id_j] +=
+                      lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_prev_values[node_id] *
+                      dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id_j, i_face));
+                  }
                 }
               }
             }
@@ -1414,9 +1563,9 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
           const size_t i_face    = face_local_number_in_its_cells[i_cell_j];
 
           b[cell_id_j] +=
-            lambda * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1] * face_boundary_values[face_id];
+            lambda * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * face_boundary_values[face_id];
           b_prev[cell_id_j] +=
-            lambda * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1] * face_boundary_prev_values[face_id];
+            lambda * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] * face_boundary_prev_values[face_id];
         }
       } else if (face_is_neumann[face_id]) {
         for (size_t i_cell_j = 0; i_cell_j < face_to_cell.size(); ++i_cell_j) {
-- 
GitLab