From 69f9926315934806834c19390118b890301ca58e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Thu, 25 Nov 2021 19:01:58 +0100
Subject: [PATCH] Fix pyramid node ordering when building the diamond mesh
---
src/mesh/DiamondDualConnectivityBuilder.cpp | 38 +++++++++++++--------
1 file changed, 24 insertions(+), 14 deletions(-)
diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 40e7a0dfc..bb9c9b3b6 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -86,21 +86,27 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells);
const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
+ const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
+ static_assert(Dimension > 1);
parallel_for(primal_number_of_faces, [&](FaceId face_id) {
const size_t i_cell = face_id;
const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
- if (primal_face_cell_list.size() == 1) {
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
- } else {
- Assert(primal_face_cell_list.size() == 2);
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
- }
-
- if constexpr (Dimension == 3) {
+ if constexpr (Dimension == 2) {
+ if (primal_face_cell_list.size() == 1) {
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
+ } else {
+ Assert(primal_face_cell_list.size() == 2);
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
+ }
+ } else if constexpr (Dimension == 3) {
if (primal_face_cell_list.size() == 1) {
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Pyramid;
+ if (primal_face_to_node_matrix[face_id].size() == 3) {
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Tetrahedron;
+ } else {
+ diamond_descriptor.cell_type_vector[i_cell] = CellType::Pyramid;
+ }
} else {
Assert(primal_face_cell_list.size() == 2);
diamond_descriptor.cell_type_vector[i_cell] = CellType::Diamond;
@@ -110,7 +116,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
- const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
const auto& cell_face_is_reversed = primal_connectivity.cellFaceIsReversed();
parallel_for(primal_number_of_faces, [&](FaceId face_id) {
@@ -129,12 +134,17 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size()] =
primal_number_of_nodes + cell_id;
- if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
- if constexpr (Dimension == 2) {
+ if constexpr (Dimension == 2) {
+ if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]);
-
- } else {
+ }
+ } else {
+ if (not cell_face_is_reversed(cell_id, i_face_in_cell)) {
+ // In 3D the basis of the pyramid is described in the
+ // indirect way IF the face is not reversed. In other words
+ // the "topological normal" must point to the "top" of the
+ // pyramid.
for (size_t i_node = 0; i_node < primal_face_node_list.size() / 2; ++i_node) {
std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node],
diamond_descriptor
--
GitLab