From 2f3941ac6c4b2c5d7a67e46166a8ea8e30ed58c5 Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Fri, 16 Sep 2022 16:46:19 +0200
Subject: [PATCH] Hybrid scheme for 3D with Cartesian corners ad straight
borders
---
src/scheme/ScalarHybridScheme.cpp | 139 ++++++++++++++++++------------
1 file changed, 85 insertions(+), 54 deletions(-)
diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index 84fed063f..a9bafc3e2 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -456,41 +456,60 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (is_boundary_node[node_id]) {
- const auto& node_to_cell = node_to_cell_matrix[node_id];
- const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
- for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
- const unsigned int i_node = node_local_number_in_its_cell[i_cell];
- const CellId cell_id = node_to_cell[i_cell];
- const auto& cell_nodes = cell_to_node_matrix[cell_id];
- if (Dimension == 2) {
- NodeId i_node_prev;
- NodeId i_node_next;
- if (i_node == 0) {
- i_node_prev = cell_nodes.size() - 1;
- } else {
- i_node_prev = i_node - 1;
- }
- if (i_node == cell_nodes.size() - 1) {
- i_node_next = 0;
- } else {
- i_node_next = i_node + 1;
- }
- const NodeId prev_node_id = cell_to_node_matrix[cell_id][i_node_prev];
- const NodeId next_node_id = cell_to_node_matrix[cell_id][i_node_next];
- if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
- compute_node_is_corner[node_id] = true;
- }
- } else if (Dimension == 3) {
- compute_node_is_corner[node_id] = true;
- } else {
- throw NotImplementedError("The scheme is not implemented in this dimension.");
- }
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ // const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ if (node_to_cell.size() == 1) {
+ compute_node_is_corner[node_id] = true;
+ //// std::cout << "node is corner : " << node_id << "; " << compute_node_is_corner[node_id] << "\n";
+ // for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ // const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+ // const CellId cell_id = node_to_cell[i_cell];
+ // const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ // NodeId i_node_prev;
+ // NodeId i_node_next;
+ // if (i_node == 0) {
+ // i_node_prev = cell_nodes.size() - 1;
+ // } else {
+ // i_node_prev = i_node - 1;
+ // }
+ // if (i_node == cell_nodes.size() - 1) {
+ // i_node_next = 0;
+ // } else {
+ // i_node_next = i_node + 1;
+ // }
+ // const NodeId prev_node_id = cell_to_node_matrix[cell_id][i_node_prev];
+ // const NodeId next_node_id = cell_to_node_matrix[cell_id][i_node_next];
+ // if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
+ // compute_node_is_corner[node_id] = true;
+ // }
+ // }
+ } else if ((Dimension == 3) && (node_to_cell.size() == 2)) {
+ compute_node_is_corner[node_id] = true;
+ //// std::cout << " node is corner droit " << node_id << "; " << compute_node_is_corner[node_id] << "\n";
}
}
});
return compute_node_is_corner;
}();
- ////std::cout << "Test 1.9 \n";
+
+ const NodeValue<const bool> node_is_straight_border = [&] {
+ NodeValue<bool> compute_node_is_straight_border{mesh->connectivity()};
+ compute_node_is_straight_border.fill(false);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (is_boundary_node[node_id]) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ if ((Dimension == 3) && (node_to_cell.size() == 2)) {
+ compute_node_is_straight_border[node_id] = true;
+ //// std::cout << "node_is_straight_border : " << node_id << "; " <<
+ // compute_node_is_straight_border[node_id]
+ // << "\n";
+ }
+ }
+ });
+ return compute_node_is_straight_border;
+ }();
+ // std::cout << "Test 1.9 \n";
const NodeValue<const bool> node_is_angle = [&] {
NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
@@ -498,7 +517,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
// for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- if (is_boundary_node[node_id]) {
+ if (is_boundary_node[node_id] && Dimension == 2) {
const auto& node_to_face = node_to_face_matrix[node_id];
TinyVector<Dimension> n1 = zero;
@@ -517,7 +536,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
if (l2Norm(n1 - n2) > (1.E-15) and l2Norm(n1 + n2) > (1.E-15) and not node_is_corner[node_id]) {
compute_node_is_angle[node_id] = true;
}
- // std::cout << node_id << " " << n1 << " " << n2 << " " << compute_node_is_angle[node_id] << "\n";
+ //// std::cout << node_id << " " << n1 << " " << n2 << " " << compute_node_is_angle[node_id] << "\n";
}
});
return compute_node_is_angle;
@@ -538,7 +557,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
const double norm_exterior_normal = l2Norm(compute_exterior_normal[node_id]);
compute_exterior_normal[node_id] *= 1. / norm_exterior_normal;
- // std::cout << node_id << " " << compute_exterior_normal[node_id] << "\n";
+ //// std::cout << node_id << " " << compute_exterior_normal[node_id] << "\n";
}
}
return compute_exterior_normal;
@@ -546,7 +565,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
// std::cout << "Test 1.11 \n";
// for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
- // std::cout << node_id << " " << exterior_normal[node_id] << " " << node_is_angle[node_id] << "\n";
+ // // std::cout << node_id << " " << exterior_normal[node_id] << " " << node_is_angle[node_id] << "\n";
//};
const FaceValue<const double> mes_l = [&] {
@@ -677,7 +696,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
compute_xj1_xj2_orth.fill(zero);
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
const auto& cell_to_face = cell_to_face_matrix[cell_id];
- // std::cout << "Test 1.16.1 \n";
+ //// std::cout << "Test 1.16.1 \n";
for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
const FaceId face_id = cell_to_face[i_face];
@@ -695,15 +714,16 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
const auto& cell_to_face = cell_to_face_matrix[cell_id];
for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
- const FaceId face_id = cell_to_face[i_face];
- TinyVector<Dimension> v;
- v[0] = xj1_xj2(cell_id, i_face)[0] + 1;
- v[1] = xj1_xj2(cell_id, i_face)[1];
- v[2] = 0;
- if (pow(dot(v, xj1_xj2(cell_id, i_face)), 2) == pow(l2Norm(v) * l2Norm(xj1_xj2(cell_id, i_face)), 2)) {
- v[2] = 1;
- }
+ const FaceId face_id = cell_to_face[i_face];
+ TinyVector<Dimension> v = zero;
if (Dimension == 3) {
+ v[0] = xj1_xj2(cell_id, i_face)[0] + 1;
+ v[1] = xj1_xj2(cell_id, i_face)[1];
+ v[2] = 0;
+ if (pow(dot(v, xj1_xj2(cell_id, i_face)), 2) == pow(l2Norm(v) * l2Norm(xj1_xj2(cell_id, i_face)), 2)) {
+ v[2] = 1;
+ }
+
// compute_w1(cell_id, i_face) = crossProduct(v, xj1_xj2(cell_id, i_face));
compute_w1(cell_id, i_face)[0] = v[1] * xj1_xj2(cell_id, i_face)[2] - v[2] * xj1_xj2(cell_id, i_face)[1];
compute_w1(cell_id, i_face)[1] = v[2] * xj1_xj2(cell_id, i_face)[0] - v[0] * xj1_xj2(cell_id, i_face)[2];
@@ -715,6 +735,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_w1;
}();
+ // std::cout << "Test 1.17.1 \n";
const FaceValuePerCell<const TinyVector<Dimension>> w2 = [&] {
FaceValuePerCell<TinyVector<Dimension>> compute_w2{mesh->connectivity()};
@@ -739,14 +760,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_w2;
}();
+ // std::cout << "Test 1.17.2 \n";
const FaceValuePerCell<const TinyVector<Dimension>> normal_face_base = [&] {
FaceValuePerCell<TinyVector<Dimension>> compute_normal_base{mesh->connectivity()};
compute_normal_base.fill(zero);
- // std::cout << "Test 1.17.01 \n";
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
const auto& cell_to_face = cell_to_face_matrix[cell_id];
- // std::cout << "Test 1.17.1 \n";
+
for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
const FaceId face_id = cell_to_face[i_face];
// TinyMatrix<Dimension> A;
@@ -756,7 +777,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
// A(0, 1) = xj1_xj2(cell_id, i_face)[0];
// A(1, 0) = xj1_xj2_orth(cell_id, i_face)[1];
// A(1, 1) = xj1_xj2(cell_id, i_face)[1];
- // std::cout << "Test 1.17.2 \n";
+
if (dot(nl[face_id], xj1_xj2(cell_id, i_face)) > 0) {
face_normal = nl[face_id];
} else {
@@ -764,7 +785,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
// compute_normal_base(cell_id, i_face) = inverse(A) * (face_kappal[face_id] * face_normal);
- // std::cout << "Test 1.17.3 \n";
+
// std::cout << face_kappal[face_id] << "\n";
// std::cout << xj1_xj2_orth(cell_id, i_face) << "\n";
// std::cout << face_normal << "\n";
@@ -785,7 +806,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
dot(face_kappal[face_id] * face_normal, xj1_xj2(cell_id, i_face)) /
pow(l2Norm(xj1_xj2(cell_id, i_face)), 2);
}
- // std::cout << "Test 1.17.4 \n";
+
// std::cout << cell_id << " " << face_id << "; alpha_l = " << compute_normal_base(cell_id, i_face)[0]
// << "; delta_l = " << compute_normal_base(cell_id, i_face)[1] << "\n";
}
@@ -938,7 +959,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
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];
}
- // std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n";
+ //// std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n";
}
return compute_node_boundary_values;
}();
@@ -1023,7 +1044,9 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ // std::cout << node_id << " " << node_is_corner[node_id] << "\n";
if (not node_is_corner[node_id]) {
+ // std::cout << node_id << " " << node_betar[node_id];
kappa_invBeta[node_id] = node_kappar[node_id] * inverse(node_betar[node_id]);
}
});
@@ -1038,10 +1061,11 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
if (node_is_corner[node_id]) {
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
const auto& node_to_cell = node_to_cell_matrix[node_id];
-
- size_t i_cell = 0;
- const CellId cell_id = node_to_cell[i_cell];
+ beta[node_id] = zero;
+ size_t i_cell = 0;
+ const CellId cell_id = node_to_cell[i_cell];
if (Dimension == 2) {
+ // std::cout << "Test 1.24.1 \n";
const unsigned int i_node = node_local_number_in_its_cell[i_cell];
const auto& cell_nodes = cell_to_node_matrix[cell_id];
@@ -1064,11 +1088,18 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
Cjr2[0] = -0.5 * (xrp2[1] - xrp1[1]);
Cjr2[1] = 0.5 * (xrp2[0] - xrp1[0]);
beta[node_id] = tensorProduct(Cjr1, (xr[node_id] - xj1)) + tensorProduct(Cjr2, (xr[node_id] - xj2));
- } else {
+ } else if (Dimension == 3 && not node_is_straight_border[node_id]) {
+ // std::cout << "Test 1.24.2 \n";
+ // std::cout << "coin : " << node_id << " " << Vj[cell_id] << "\n";
const TinyMatrix<Dimension> I = identity;
beta[node_id] = 0.125 * Vj[cell_id] * I;
// throw NotImplementedError("The scheme is not implemented in this dimension.");
+ } else if (Dimension == 3 && node_is_straight_border[node_id]) {
+ // std::cout << "Test 1.24.3 \n";
+ // std::cout << " bord droit : " << node_id << " " << Vj[cell_id] << "\n";
+ const TinyMatrix<Dimension> I = identity;
+ beta[node_id] += 0.125 * Vj[cell_id] * I;
}
}
});
@@ -1383,7 +1414,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
};
- // std::cout << "S1 = " << S1 << "\n";
+ //// std::cout << "S1 = " << S1 << "\n";
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
S1(cell_id, cell_id) += Vj[cell_id];
--
GitLab