From 7bc8d0dba628ec5027c334015721c0079832ac3e Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Mon, 18 Jul 2022 10:19:48 +0200
Subject: [PATCH] Change of the angles calculation for Neumann BC
---
src/scheme/ScalarNodalScheme.cpp | 46 +++++++++++++++++++++++++++++---
1 file changed, 43 insertions(+), 3 deletions(-)
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 9367249f7..33baf960f 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -313,6 +313,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
+ const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
@@ -320,6 +321,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const auto& node_local_numbers_in_their_cells = mesh->connectivity().nodeLocalNumbersInTheirCells();
const CellValue<const double> Vj = mesh_data.Vj();
const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& nl = mesh_data.nl();
const NodeValue<const bool> node_is_neumann = [&] {
NodeValue<bool> compute_node_is_neumann{mesh->connectivity()};
@@ -409,6 +411,35 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return compute_node_is_corner;
}();
+ const NodeValue<const bool> node_is_angle = [&] {
+ NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
+ compute_node_is_angle.fill(false);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (is_boundary_node[node_id]) {
+ const auto& node_to_face = node_to_face_matrix[node_id];
+
+ TinyVector<Dimension> n1 = zero;
+ TinyVector<Dimension> n2 = zero;
+
+ for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+ FaceId face_id = node_to_face[i_face];
+ if (is_boundary_face[face_id]) {
+ if (l2Norm(n1) == 0) {
+ n1 = nl[face_id];
+ } else {
+ n2 = nl[face_id];
+ }
+ }
+ if (n1 != n2) {
+ compute_node_is_angle[node_id] = true;
+ }
+ }
+ }
+ });
+ return compute_node_is_angle;
+ }();
+
const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] {
NodeValue<TinyVector<Dimension>> compute_exterior_normal{mesh->connectivity()};
compute_exterior_normal.fill(zero);
@@ -429,6 +460,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return compute_exterior_normal;
}();
+ // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ // std::cout << node_id << " " << exterior_normal[node_id] << "\n";
+ //};
+
const FaceValue<const double> mes_l = [&] {
if constexpr (Dimension == 1) {
FaceValue<double> compute_mes_l{mesh->connectivity()};
@@ -459,8 +494,13 @@ 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 (not node_is_dirichlet[node_id]) {
- compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
- sum_mes_l[node_id] += mes_l[face_id];
+ if (node_is_angle[node_id]) {
+ compute_node_boundary_values[node_id] +=
+ value_list[i_face] * std::abs(dot(nl[face_id], exterior_normal[node_id]));
+ } else {
+ compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
+ sum_mes_l[node_id] += mes_l[face_id];
+ }
}
}
}
@@ -487,7 +527,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
boundary_condition);
}
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
- if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id])) {
+ if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
} else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
--
GitLab