From f0c17b59be66e8f7b777d1c93ee10ed3a60791a3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 16 Apr 2021 14:36:26 +0200
Subject: [PATCH] Fix mass term discretization
---
src/scheme/ScalarDiamondScheme.cpp | 17 +++++++++--------
1 file changed, 9 insertions(+), 8 deletions(-)
diff --git a/src/scheme/ScalarDiamondScheme.cpp b/src/scheme/ScalarDiamondScheme.cpp
index 889a1054a..464d94384 100644
--- a/src/scheme/ScalarDiamondScheme.cpp
+++ b/src/scheme/ScalarDiamondScheme.cpp
@@ -604,9 +604,6 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
return computed_alpha_l;
}();
- // double lambda = (Tf == 0) ? 0 : 1;
- double time_factor = 1; //(lambda == 0) ? 1 : dt;
-
const CellValue<const double> primal_Vj = mesh_data.Vj();
SparseMatrixDescriptor S{number_of_dof};
@@ -618,18 +615,22 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
const CellId cell_id2 = primal_face_to_cell[j_cell];
if (i_cell == j_cell) {
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) +=
- (time_factor * beta_l + (*alpha)[cell_id1] * primal_Vj[cell_id1]);
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= beta_l;
}
} else {
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= time_factor * beta_l;
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
}
}
}
}
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ const size_t j = cell_dof_number[cell_id];
+ S(j, j) += (*alpha)[cell_id] * primal_Vj[cell_id];
+ }
+
const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
@@ -662,7 +663,7 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
- S(cell_dof_number[i_id], cell_dof_number[j_id]) -= time_factor * w_rj(node_id, j_cell) * a_ir;
+ S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * a_ir;
}
@@ -671,7 +672,7 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
FaceId l_id = node_to_face_matrix[node_id][l_face];
if (primal_face_is_on_boundary[l_id]) {
- S(cell_dof_number[i_id], face_dof_number[l_id]) -= time_factor * w_rl(node_id, l_face) * a_ir;
+ S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * a_ir;
}
--
GitLab