Test matrix with the flux at the faces

src/scheme/ScalarNodalScheme.cpp

@@ -313,6 +313,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
     const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
 
+    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();
     const auto& cell_to_node_matrix               = mesh->connectivity().cellToNodeMatrix();
@@ -736,6 +737,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       non_zeros.fill(4);   // Modif pour optimiser
       CRSMatrixDescriptor<double> S1(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
       CRSMatrixDescriptor<double> S2(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
+      CRSMatrixDescriptor<double> C(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
+
       for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
         const auto& node_to_cell = node_to_cell_matrix[node_id];
         for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
@@ -744,6 +747,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             const CellId cell_id_k   = node_to_cell[i_cell_k];
             S1(cell_id_j, cell_id_k) = 0;
             S2(cell_id_j, cell_id_k) = 0;
+            C(cell_id_j, cell_id_k)  = 0;
           }
         }
       }
@@ -831,6 +835,32 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         }
       }
 
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        if (not is_boundary_face[face_id]) {
+          const auto& face_to_cell = face_to_cell_matrix[face_id];
+          CellId cell_id_j         = face_to_cell[0];
+          CellId cell_id_k         = face_to_cell[1];
+          C(cell_id_j, cell_id_j) += dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xj[cell_id_k]);
+          C(cell_id_k, cell_id_k) += dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xj[cell_id_k]);
+          C(cell_id_j, cell_id_k) -= dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xj[cell_id_k]);
+          C(cell_id_k, cell_id_j) -= dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xj[cell_id_k]);
+        } else {
+          const auto& face_to_cell = face_to_cell_matrix[face_id];
+          CellId cell_id_j         = face_to_cell[0];
+          C(cell_id_j, cell_id_j) += dt * mes_l[face_id] / l2Norm(xj[cell_id_j] - xl[face_id]);
+        }
+      }
+
+      const double epsilon = 0.5;
+
+      for (CellId cell_id_j = 0; cell_id_j < mesh->numberOfCells(); ++cell_id_j) {
+        for (CellId cell_id_k = 0; cell_id_k < mesh->numberOfCells(); ++cell_id_k) {
+          if (S1(cell_id_j, cell_id_k) != 0 or C(cell_id_j, cell_id_k) != 0) {
+            S1(cell_id_j, cell_id_k) = (1. - epsilon) * S1(cell_id_j, cell_id_k) + epsilon * C(cell_id_j, cell_id_k);
+          }
+        }
+      }
+
       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];
@@ -849,6 +879,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       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) {