diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index d1f22b1442754f35a06f8e7bae2e1a67d149bdf6..9c442d0329900c6bf49f52ca050db4dc6cef511d 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -353,7 +353,9 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const NodeValue<const double> node_boundary_values = [&] { NodeValue<double> compute_node_boundary_values; + NodeValue<double> sum_mes_l; compute_node_boundary_values.fill(0); + sum_mes_l.fill(0); for (const auto& boundary_condition : boundary_condition_list) { std::visit( [&](auto&& bc) { @@ -369,10 +371,16 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) { const NodeId node_id = face_nodes[i_node]; if (node_is_dirichlet[node_id] == false) { - compute_node_boundary_values[node_id] += 0.5 * value_list[i_face] * mes_l[face_id]; + compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id]; + sum_mes_l[node_id] += mes_l[face_id]; } } } + for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { + if (not node_is_dirichlet[node_id]) { + compute_node_boundary_values[node_id] /= sum_mes_l[node_id]; + } + } } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) { const auto& face_list = bc.faceList(); const auto& value_list = bc.valueList(); @@ -535,7 +543,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } - } else if (not node_is_corner[node_id]) { + } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id])) { const TinyMatrix<Dimension> I = identity; const TinyVector<Dimension> n = exterior_normal[node_id]; const TinyMatrix<Dimension> nxn = tensorProduct(n, n); @@ -565,6 +573,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand CellValue<const double> fj = f->cellValues(); CRSMatrix A{S.getCRSMatrix()}; + Vector<double> b{mesh->numberOfCells()}; b = zero; for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {