diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp index 40e7a0dfc21e1b093a560d92ca3848c99e63307f..bb9c9b3b643c87c4a71de31cf6a2a5c60e53b0e2 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