diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 104c89351e05b286ae5fb5b660a6dd50866fa603..7fd220f64f08b536255d45ec89fbb6bfa46fbbd0 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/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; }