Skip to content
Snippets Groups Projects
Commit 1a66c65f authored by Axelle Drouard's avatar Axelle Drouard
Browse files

Correction node boundary values for Neumann BC

parent 6a6be095
No related branches found
No related tags found
No related merge requests found
......@@ -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,9 +371,15 @@ 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();
......@@ -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) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment