Merge branch 'feature/Nodal_diffusion' of gitlab.delpinux.fr:code/pugs into feature/Nodal_diffusion

src/scheme/ScalarNodalScheme.cpp

@@ -586,8 +586,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               const CellId cell_id_k = node_to_cell[i_cell_k];
               const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
-              S(cell_id_j, cell_id_k) -=
-                dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+              S(cell_id_j, cell_id_k) +=
+                dt / Vj[cell_id_j] * dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
             }
           }
         } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id]) && (not node_is_dirichlet[node_id])) {
@@ -604,7 +604,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               const CellId cell_id_k = node_to_cell[i_cell_k];
               const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
-              S(cell_id_j, cell_id_k) -=
+              S(cell_id_j, cell_id_k) +=
+                dt / Vj[cell_id_j] *
                 dot(kappar_invBetar[node_id] *
                       (Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
                     Q * Cjr(cell_id_j, i_node_j));
@@ -619,8 +620,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               const CellId cell_id_k = node_to_cell[i_cell_k];
               const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
-              S(cell_id_j, cell_id_k) -=
-                dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+              S(cell_id_j, cell_id_k) +=
+                dt / Vj[cell_id_j] * dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
             }
           }
         } else if (node_is_dirichlet[node_id] && node_is_corner[node_id]) {
@@ -632,19 +633,14 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               const CellId cell_id_k = node_to_cell[i_cell_k];
               const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
-              S(cell_id_j, cell_id_k) -=
+              S(cell_id_j, cell_id_k) +=
+                dt / Vj[cell_id_j] *
                 dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
             }
           }
         }
       }
 
-      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) {
-          S(cell_id_j, cell_id_k) = -dt / Vj[cell_id_j] * S(cell_id_j, cell_id_k);
-        }
-      };
-
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
         S(cell_id, cell_id) += 1;
       }