diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index f9767abd0ac47e85455b0ae2c58f08998200aa5a..f380a0fff001b0b20aee1d39bf4d3382e53cbf51 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/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) {