Fix mass term discretization

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,16 +615,20 @@ 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();
@@ -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;
                         }