diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp index d113bcc2341201dc4070fc5c0e4c1fdd46f33289..f48d00f8acfb76fa06d03f4eb8255c6936474921 100644 --- a/src/mesh/DiamondDualMeshBuilder.cpp +++ b/src/mesh/DiamondDualMeshBuilder.cpp @@ -16,36 +16,70 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D ConnectivityDescriptor& diamond_descriptor) { const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes(); + const size_t primal_number_of_faces = primal_connectivity.numberOfFaces(); const size_t primal_number_of_cells = primal_connectivity.numberOfCells(); - const size_t diamond_number_of_nodes = primal_number_of_cells + primal_number_of_nodes; + const size_t diamond_number_of_nodes = [&]() { + if constexpr (Dimension == 1) { + return 2 * primal_number_of_nodes - primal_number_of_cells; + } else { + return primal_number_of_cells + primal_number_of_nodes; + } + }(); diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes); - const auto& primal_node_number = primal_connectivity.nodeNumber(); + if constexpr (Dimension == 1) { + const auto& primal_node_number = primal_connectivity.nodeNumber(); - for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { - diamond_descriptor.node_number_vector[primal_node_id] = primal_node_number[primal_node_id]; - } + const size_t diamond_number_of_cells = primal_number_of_faces; + const size_t number_of_kept_nodes = 2 * (diamond_number_of_nodes - diamond_number_of_cells); - const auto& primal_cell_number = primal_connectivity.cellNumber(); + const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); + size_t next_kept_node_id = 0; - const size_t max_node_number = max(primal_node_number); + for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) { + const auto& primal_node_cell_list = primal_node_to_cell_matrix[node_id]; + if (primal_node_cell_list.size() == 1) { + diamond_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[node_id]; + } + } - for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) { - diamond_descriptor.node_number_vector[primal_number_of_nodes + primal_cell_id] = - primal_cell_number[primal_cell_id] + max_node_number; - } + if (number_of_kept_nodes != next_kept_node_id) { + throw UnexpectedError("unexpected number of kept node" + std::to_string(next_kept_node_id) + + " != " + std::to_string(number_of_kept_nodes)); + } - const size_t diamond_number_of_cells = primal_connectivity.numberOfFaces(); - diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells); + } else { + const auto& primal_node_number = primal_connectivity.nodeNumber(); - const auto& primal_face_number = primal_connectivity.faceNumber(); + for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { + diamond_descriptor.node_number_vector[primal_node_id] = primal_node_number[primal_node_id]; + } + + const auto& primal_cell_number = primal_connectivity.cellNumber(); - for (FaceId i_primal_face = 0; i_primal_face < primal_connectivity.numberOfFaces(); ++i_primal_face) { - diamond_descriptor.cell_number_vector[i_primal_face] = primal_face_number[i_primal_face]; + const size_t max_node_number = max(primal_node_number); + + for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) { + diamond_descriptor.node_number_vector[primal_number_of_nodes + primal_cell_id] = + primal_cell_number[primal_cell_id] + max_node_number; + } } + const size_t diamond_number_of_cells = primal_number_of_faces; + diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells); + if constexpr (Dimension == 1) { + const auto& primal_node_number = primal_connectivity.nodeNumber(); + for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) { + diamond_descriptor.cell_number_vector[primal_node_id] = primal_node_number[primal_node_id]; + } + } else { + const auto& primal_face_number = primal_connectivity.faceNumber(); + for (FaceId i_primal_face = 0; i_primal_face < primal_number_of_faces; ++i_primal_face) { + diamond_descriptor.cell_number_vector[i_primal_face] = primal_face_number[i_primal_face]; + } + } if constexpr (Dimension == 3) { const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size(); diamond_descriptor.edge_number_vector.resize(number_of_edges); @@ -57,8 +91,6 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D } } - diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells); - diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells); const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix(); @@ -68,7 +100,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face]; if constexpr (Dimension == 1) { - throw NotImplementedError("dimension 1"); + diamond_descriptor.cell_type_vector[i_cell] = CellType::Line; } else if constexpr (Dimension == 2) { if (primal_face_cell_list.size() == 1) { diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle; @@ -90,68 +122,95 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D 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(); - for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) { - const size_t& i_diamond_cell = i_face; - const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face]; - const auto& primal_face_node_list = primal_face_to_node_matrix[i_face]; - if (primal_face_cell_list.size() == 1) { - diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 1); + if constexpr (Dimension == 1) { + const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); + const size_t number_of_kept_nodes = 2 * (diamond_number_of_nodes - diamond_number_of_cells); + size_t next_kept_node_id = 0; + + const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells(); - const CellId cell_id = primal_face_cell_list[0]; - const auto i_face_in_cell = primal_face_local_number_in_their_cells(i_face, 0); + for (NodeId i_node = 0; i_node < primal_connectivity.numberOfNodes(); ++i_node) { + const size_t& i_diamond_cell = i_node; + const auto& primal_node_cell_list = primal_node_to_cell_matrix[i_node]; + diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(2); + if (primal_node_cell_list.size() == 1) { + const CellId cell_id = primal_node_cell_list[0]; + const auto i_node_in_cell = primal_node_local_number_in_their_cells(i_node, 0); - for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) { - diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node] = primal_face_node_list[i_node]; + diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node_in_cell] = next_kept_node_id++; + diamond_descriptor.cell_to_node_vector[i_diamond_cell][1 - i_node_in_cell] = + number_of_kept_nodes + primal_node_cell_list[0]; + } else { + diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = number_of_kept_nodes + primal_node_cell_list[0]; + diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = number_of_kept_nodes + primal_node_cell_list[1]; } - diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size()] = - primal_number_of_nodes + cell_id; + } + } else { + static_assert(Dimension > 1, "invalid dimension"); + + 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(); + for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) { + const size_t& i_diamond_cell = i_face; + const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face]; + const auto& primal_face_node_list = primal_face_to_node_matrix[i_face]; + if (primal_face_cell_list.size() == 1) { + diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 1); - if (cell_face_is_reversed(cell_id, i_face_in_cell)) { - if constexpr (Dimension == 2) { - std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0], - diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]); + const CellId cell_id = primal_face_cell_list[0]; + const auto i_face_in_cell = primal_face_local_number_in_their_cells(i_face, 0); - } else { - 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 - .cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() - 1 - i_node]); - } + for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) { + diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node] = primal_face_node_list[i_node]; } - } - } else { - Assert(primal_face_cell_list.size() == 2); - diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 2); - - const CellId cell0_id = primal_face_cell_list[0]; - const CellId cell1_id = primal_face_cell_list[1]; - const auto i_face_in_cell = primal_face_local_number_in_their_cells(i_face, 0); - - if constexpr (Dimension == 2) { - Assert(primal_face_node_list.size() == 2); - diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id; - diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[0]; - diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell1_id; - diamond_descriptor.cell_to_node_vector[i_diamond_cell][3] = primal_face_node_list[1]; - - if (cell_face_is_reversed(cell0_id, i_face_in_cell)) { - std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1], - diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]); + 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) { + std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0], + diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]); + + } else { + 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 + .cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() - 1 - i_node]); + } + } } } else { - diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id; - for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) { - diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node + 1] = primal_face_node_list[i_node]; - } - diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1] = - primal_number_of_nodes + cell1_id; + Assert(primal_face_cell_list.size() == 2); + diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 2); + + const CellId cell0_id = primal_face_cell_list[0]; + const CellId cell1_id = primal_face_cell_list[1]; + const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(i_face, 0); - if (cell_face_is_reversed(cell0_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][primal_face_node_list.size() + 1]); + if constexpr (Dimension == 2) { + Assert(primal_face_node_list.size() == 2); + diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id; + diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[0]; + diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell1_id; + diamond_descriptor.cell_to_node_vector[i_diamond_cell][3] = primal_face_node_list[1]; + + if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) { + std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1], + diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]); + } + } else { + diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id; + for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) { + diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node + 1] = primal_face_node_list[i_node]; + } + diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1] = + primal_number_of_nodes + cell1_id; + + if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) { + std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0], + diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1]); + } } } } @@ -174,9 +233,11 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor); - MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor); - if constexpr (Dimension == 3) { - MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor); + if constexpr (Dimension > 1) { + MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor); + if constexpr (Dimension > 2) { + MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor); + } } { @@ -329,19 +390,44 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> diamond_descriptor.node_owner_vector.resize(diamond_descriptor.node_number_vector.size()); - const auto& primal_node_owner = primal_connectivity.nodeOwner(); - for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { - diamond_descriptor.node_owner_vector[primal_node_id] = primal_node_owner[primal_node_id]; - } - const auto& primal_cell_owner = primal_connectivity.cellOwner(); - for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) { - diamond_descriptor.node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id]; + if constexpr (Dimension == 1) { + const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); + const auto& primal_node_owner = primal_connectivity.nodeOwner(); + size_t next_kept_node_id = 0; + for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { + if (node_to_cell_matrix[primal_node_id].size() == 1) { + diamond_descriptor.node_owner_vector[next_kept_node_id++] = primal_node_owner[primal_node_id]; + } + } + const size_t number_of_kept_nodes = next_kept_node_id; + const auto& primal_cell_owner = primal_connectivity.cellOwner(); + for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) { + diamond_descriptor.node_owner_vector[number_of_kept_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id]; + } + } else { + const auto& primal_node_owner = primal_connectivity.nodeOwner(); + for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { + diamond_descriptor.node_owner_vector[primal_node_id] = primal_node_owner[primal_node_id]; + } + const auto& primal_cell_owner = primal_connectivity.cellOwner(); + for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) { + diamond_descriptor.node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id]; + } } - diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size()); - const auto& primal_face_owner = primal_connectivity.faceOwner(); - for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_cells; ++primal_face_id) { - diamond_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id]; + if constexpr (Dimension == 1) { + diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size()); + const auto& primal_node_owner = primal_connectivity.nodeOwner(); + for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) { + diamond_descriptor.cell_owner_vector[primal_node_id] = primal_node_owner[primal_node_id]; + } + } else { + diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size()); + const size_t primal_number_of_faces = primal_connectivity.numberOfFaces(); + const auto& primal_face_owner = primal_connectivity.faceOwner(); + for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_faces; ++primal_face_id) { + diamond_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id]; + } } { @@ -391,13 +477,28 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> { #warning define transfer functions - NodeId i_node = 0; - for (; i_node < primal_number_of_nodes; ++i_node) { - diamond_xr[i_node] = primal_xr[i_node]; - } + if constexpr (Dimension == 1) { + const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); + NodeId next_node_id = 0; + for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { + if (node_to_cell_matrix[primal_node_id].size() == 1) { + diamond_xr[next_node_id++] = primal_xr[primal_node_id]; + } + } + + for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) { + diamond_xr[next_node_id++] = primal_xj[i_cell]; + } + + } else { + NodeId i_node = 0; + for (; i_node < primal_number_of_nodes; ++i_node) { + diamond_xr[i_node] = primal_xr[i_node]; + } - for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) { - diamond_xr[i_node++] = primal_xj[i_cell]; + for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) { + diamond_xr[i_node++] = primal_xj[i_cell]; + } } } @@ -418,13 +519,6 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh> m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr); } -template <> -[[deprecated]] void -DiamondDualMeshBuilder::_buildDiamondMeshFrom<1>(const std::shared_ptr<const IMesh>&) -{ - m_mesh = 0; -} - DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh) { switch (p_mesh->dimension()) {