diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 72936fe9222f2a9ac273f2190f37d79cf1513b39..104c89351e05b286ae5fb5b660a6dd50866fa603 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -586,7 +586,7 @@ 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) -= dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } @@ -604,10 +604,10 @@ 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) -= 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_k, i_node_j)); + Q * Cjr(cell_id_j, i_node_j)); } } } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) { @@ -619,7 +619,7 @@ 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) -= dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } @@ -632,7 +632,7 @@ 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) -= dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } @@ -641,7 +641,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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) = -dt / Vj[cell_id_j] * S(cell_id_j, cell_id_k); } }; @@ -720,23 +720,23 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand } } } - }; + } + }; - Vector<double> T{mesh->numberOfCells()}; - T = zero; + Vector<double> T{mesh->numberOfCells()}; + T = zero; - Vector<double> B{mesh->numberOfCells()}; - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; }); + Vector<double> B{mesh->numberOfCells()}; + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; }); - LinearSolver solver; - solver.solveLocalSystem(A, T, B); + LinearSolver solver; + solver.solveLocalSystem(A, T, B); - m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh); - auto& solution = *m_solution; - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; }); - } + m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh); + auto& solution = *m_solution; + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; }); } } };