From 24c5be377cc73b23ad9adf67052c06a0d49c5169 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com> Date: Tue, 4 Jan 2022 20:13:52 +0100 Subject: [PATCH] Add dual 1d connectivity This is shared by direct 1d dual, 1d diamond dual and 1d median dual. --- src/mesh/CMakeLists.txt | 1 + src/mesh/DiamondDualConnectivityBuilder.cpp | 160 +----------- src/mesh/DiamondDualMeshBuilder.cpp | 1 + src/mesh/Dual1DConnectivityBuilder.cpp | 247 ++++++++++++++++++ src/mesh/Dual1DConnectivityBuilder.hpp | 39 +++ src/mesh/DualConnectivityManager.cpp | 61 ++++- src/mesh/DualConnectivityManager.hpp | 11 + src/mesh/DualMeshType.hpp | 10 +- src/mesh/MedianDualConnectivityBuilder.cpp | 6 +- ...malToDiamondDualConnectivityDataMapper.hpp | 3 + .../PrimalToDual1DConnectivityDataMapper.hpp | 140 ++++++++++ ...imalToMedianDualConnectivityDataMapper.hpp | 3 + tests/CMakeLists.txt | 1 + ...malToDiamondDualConnectivityDataMapper.cpp | 139 +--------- ...t_PrimalToDual1DConnectivityDataMapper.cpp | 160 ++++++++++++ ...imalToMedianDualConnectivityDataMapper.cpp | 28 ++ 16 files changed, 711 insertions(+), 299 deletions(-) create mode 100644 src/mesh/Dual1DConnectivityBuilder.cpp create mode 100644 src/mesh/Dual1DConnectivityBuilder.hpp create mode 100644 src/mesh/PrimalToDual1DConnectivityDataMapper.hpp create mode 100644 tests/test_PrimalToDual1DConnectivityDataMapper.cpp diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt index 873a94dd3..b85297c51 100644 --- a/src/mesh/CMakeLists.txt +++ b/src/mesh/CMakeLists.txt @@ -9,6 +9,7 @@ add_library( ConnectivityDispatcher.cpp DiamondDualConnectivityBuilder.cpp DiamondDualMeshBuilder.cpp + Dual1DConnectivityBuilder.cpp DualConnectivityManager.cpp DualMeshManager.cpp GmshReader.cpp diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp index ec1620fab..3b8d6c836 100644 --- a/src/mesh/DiamondDualConnectivityBuilder.cpp +++ b/src/mesh/DiamondDualConnectivityBuilder.cpp @@ -17,6 +17,8 @@ void DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity, ConnectivityDescriptor& diamond_descriptor) { + static_assert((Dimension == 2) or (Dimension == 3), "invalid connectivity dimension"); + 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(); @@ -188,93 +190,11 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec }); } -template <> -void -DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Connectivity<1>& primal_connectivity, - 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 = 2 * primal_number_of_nodes - primal_number_of_cells; - - 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_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); - size_t next_kept_node_id = 0; - - diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes); - - const auto& primal_node_number = primal_connectivity.nodeNumber(); - - for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { - const auto& primal_node_cell_list = primal_node_to_cell_matrix[primal_node_id]; - if (primal_node_cell_list.size() == 1) { - diamond_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[primal_node_id]; - } - } - - const size_t primal_number_of_kept_nodes = next_kept_node_id; - - const auto& primal_cell_number = primal_connectivity.cellNumber(); - - const size_t cell_number_shift = max(primal_node_number) + 1; - - for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) { - diamond_descriptor.node_number_vector[primal_number_of_kept_nodes + primal_cell_id] = - primal_cell_number[primal_cell_id] + cell_number_shift; - } - - 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)); - } - - diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells); - - 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]; - } - - diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells); - - for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) { - const size_t i_cell = node_id; - - diamond_descriptor.cell_type_vector[i_cell] = CellType::Line; - } - - diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells); - - const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells(); - - { - size_t next_kept_node_id = 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 auto i_node_in_cell = primal_node_local_number_in_their_cells(i_node, 0); - - 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]; - } - } - } -} - template <size_t Dimension> void DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivity& i_primal_connectivity) { - static_assert(Dimension <= 3, "invalid connectivity dimension"); + static_assert((Dimension == 2) or (Dimension == 3), "invalid connectivity dimension"); using ConnectivityType = Connectivity<Dimension>; const ConnectivityType& primal_connectivity = dynamic_cast<const ConnectivityType&>(i_primal_connectivity); @@ -283,11 +203,9 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor); - if constexpr (Dimension > 1) { - ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor); - if constexpr (Dimension > 2) { - ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor); - } + ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor); + if constexpr (Dimension == 3) { + ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor); } { @@ -328,7 +246,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit } } - if constexpr (Dimension > 1) { + { const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix(); using Face = ConnectivityFace<Dimension>; @@ -381,7 +299,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit } } - if constexpr (Dimension > 2) { + if constexpr (Dimension == 3) { const auto& primal_edge_to_node_matrix = primal_connectivity.edgeToNodeMatrix(); using Edge = ConnectivityFace<2>; @@ -437,21 +355,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit diamond_descriptor.node_owner_vector.resize(diamond_descriptor.node_number_vector.size()); - 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]; @@ -462,13 +366,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit } } - 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(); @@ -515,40 +413,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit m_connectivity = ConnectivityType::build(diamond_descriptor); - if constexpr (Dimension == 1) { - const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); - - NodeId dual_node_id = 0; - m_primal_node_to_dual_node_map = [&]() { - std::vector<std::pair<NodeId, NodeId>> primal_node_to_dual_node_vector; - 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) { - primal_node_to_dual_node_vector.push_back(std::make_pair(primal_node_id, dual_node_id++)); - } - } - return convert_to_array(primal_node_to_dual_node_vector); - }(); - - m_primal_cell_to_dual_node_map = [&]() { - CellIdToNodeIdMap primal_cell_to_dual_node_map{primal_number_of_cells}; - for (CellId primal_cell_id = 0; primal_cell_id < primal_cell_to_dual_node_map.size(); ++primal_cell_id) { - primal_cell_to_dual_node_map[primal_cell_id] = std::make_pair(primal_cell_id, dual_node_id++); - } - return primal_cell_to_dual_node_map; - }(); - - m_primal_face_to_dual_cell_map = [&]() { - FaceIdToCellIdMap primal_face_to_dual_cell_map{primal_connectivity.numberOfFaces()}; - for (size_t id = 0; id < primal_face_to_dual_cell_map.size(); ++id) { - const CellId dual_cell_id = id; - const FaceId primal_face_id = id; - - primal_face_to_dual_cell_map[id] = std::make_pair(primal_face_id, dual_cell_id); - } - return primal_face_to_dual_cell_map; - }(); - } - m_mapper = std::make_shared<PrimalToDiamondDualConnectivityDataMapper<Dimension>>(primal_connectivity, dynamic_cast<const ConnectivityType&>( @@ -564,10 +428,6 @@ DiamondDualConnectivityBuilder::DiamondDualConnectivityBuilder(const IConnectivi throw NotImplementedError("Construction of diamond dual mesh is not implemented in parallel"); } switch (connectivity.dimension()) { - case 1: { - this->_buildDiamondConnectivityFrom<1>(connectivity); - break; - } case 2: { this->_buildDiamondConnectivityFrom<2>(connectivity); break; diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp index 48a4646be..5db07c20f 100644 --- a/src/mesh/DiamondDualMeshBuilder.cpp +++ b/src/mesh/DiamondDualMeshBuilder.cpp @@ -7,6 +7,7 @@ #include <mesh/MeshData.hpp> #include <mesh/MeshDataManager.hpp> #include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp> +#include <mesh/PrimalToDual1DConnectivityDataMapper.hpp> template <size_t Dimension> void diff --git a/src/mesh/Dual1DConnectivityBuilder.cpp b/src/mesh/Dual1DConnectivityBuilder.cpp new file mode 100644 index 000000000..aa7dcdd0f --- /dev/null +++ b/src/mesh/Dual1DConnectivityBuilder.cpp @@ -0,0 +1,247 @@ +#include <mesh/Dual1DConnectivityBuilder.hpp> + +#include <mesh/Connectivity.hpp> +#include <mesh/ConnectivityDescriptor.hpp> +#include <mesh/ConnectivityDispatcher.hpp> +#include <mesh/ItemValueUtils.hpp> +#include <mesh/Mesh.hpp> +#include <mesh/PrimalToDual1DConnectivityDataMapper.hpp> +#include <mesh/RefId.hpp> +#include <utils/Array.hpp> +#include <utils/Messenger.hpp> + +#include <vector> + +void +Dual1DConnectivityBuilder::_buildConnectivityDescriptor(const Connectivity<1>& primal_connectivity, + ConnectivityDescriptor& dual_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 dual_number_of_nodes = 2 * primal_number_of_nodes - primal_number_of_cells; + + const size_t dual_number_of_cells = primal_number_of_faces; + const size_t number_of_kept_nodes = 2 * (dual_number_of_nodes - dual_number_of_cells); + + const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); + size_t next_kept_node_id = 0; + + dual_descriptor.node_number_vector.resize(dual_number_of_nodes); + + const auto& primal_node_number = primal_connectivity.nodeNumber(); + + for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) { + const auto& primal_node_cell_list = primal_node_to_cell_matrix[primal_node_id]; + if (primal_node_cell_list.size() == 1) { + dual_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[primal_node_id]; + } + } + + const size_t primal_number_of_kept_nodes = next_kept_node_id; + + const auto& primal_cell_number = primal_connectivity.cellNumber(); + + const size_t cell_number_shift = max(primal_node_number) + 1; + + for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) { + dual_descriptor.node_number_vector[primal_number_of_kept_nodes + primal_cell_id] = + primal_cell_number[primal_cell_id] + cell_number_shift; + } + + 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)); + } + + dual_descriptor.cell_number_vector.resize(dual_number_of_cells); + + for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) { + dual_descriptor.cell_number_vector[primal_node_id] = primal_node_number[primal_node_id]; + } + + dual_descriptor.cell_type_vector.resize(dual_number_of_cells); + + for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) { + const size_t i_cell = node_id; + + dual_descriptor.cell_type_vector[i_cell] = CellType::Line; + } + + dual_descriptor.cell_to_node_vector.resize(dual_number_of_cells); + + const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells(); + + { + size_t next_kept_node_id = 0; + for (NodeId i_node = 0; i_node < primal_connectivity.numberOfNodes(); ++i_node) { + const size_t& i_dual_cell = i_node; + const auto& primal_node_cell_list = primal_node_to_cell_matrix[i_node]; + dual_descriptor.cell_to_node_vector[i_dual_cell].resize(2); + if (primal_node_cell_list.size() == 1) { + const auto i_node_in_cell = primal_node_local_number_in_their_cells(i_node, 0); + + dual_descriptor.cell_to_node_vector[i_dual_cell][i_node_in_cell] = next_kept_node_id++; + dual_descriptor.cell_to_node_vector[i_dual_cell][1 - i_node_in_cell] = + number_of_kept_nodes + primal_node_cell_list[0]; + } else { + dual_descriptor.cell_to_node_vector[i_dual_cell][0] = number_of_kept_nodes + primal_node_cell_list[0]; + dual_descriptor.cell_to_node_vector[i_dual_cell][1] = number_of_kept_nodes + primal_node_cell_list[1]; + } + } + } +} + +void +Dual1DConnectivityBuilder::_buildConnectivityFrom(const IConnectivity& i_primal_connectivity) +{ + using ConnectivityType = Connectivity<1>; + + const ConnectivityType& primal_connectivity = dynamic_cast<const ConnectivityType&>(i_primal_connectivity); + + ConnectivityDescriptor dual_descriptor; + + this->_buildConnectivityDescriptor(primal_connectivity, dual_descriptor); + + { + const std::unordered_map<unsigned int, NodeId> node_to_id_map = [&] { + std::unordered_map<unsigned int, NodeId> node_to_id_map; + for (size_t i_node = 0; i_node < dual_descriptor.node_number_vector.size(); ++i_node) { + node_to_id_map[dual_descriptor.node_number_vector[i_node]] = i_node; + } + return node_to_id_map; + }(); + + for (size_t i_node_list = 0; i_node_list < primal_connectivity.template numberOfRefItemList<ItemType::node>(); + ++i_node_list) { + const auto& primal_ref_node_list = primal_connectivity.template refItemList<ItemType::node>(i_node_list); + const auto& primal_node_list = primal_ref_node_list.list(); + + const std::vector<NodeId> dual_node_list = [&]() { + std::vector<NodeId> dual_node_list; + + for (size_t i_primal_node = 0; i_primal_node < primal_node_list.size(); ++i_primal_node) { + auto primal_node_id = primal_node_list[i_primal_node]; + const auto i_dual_node = node_to_id_map.find(primal_connectivity.nodeNumber()[primal_node_id]); + if (i_dual_node != node_to_id_map.end()) { + dual_node_list.push_back(i_dual_node->second); + } + } + + return dual_node_list; + }(); + + if (parallel::allReduceOr(dual_node_list.size() > 0)) { + Array<NodeId> node_array(dual_node_list.size()); + for (size_t i = 0; i < dual_node_list.size(); ++i) { + node_array[i] = dual_node_list[i]; + } + dual_descriptor.addRefItemList(RefNodeList{primal_ref_node_list.refId(), node_array}); + } + } + } + + const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes(); + const size_t primal_number_of_cells = primal_connectivity.numberOfCells(); + + dual_descriptor.node_owner_vector.resize(dual_descriptor.node_number_vector.size()); + + { + 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) { + dual_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) { + dual_descriptor.node_owner_vector[number_of_kept_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id]; + } + } + + { + dual_descriptor.cell_owner_vector.resize(dual_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) { + dual_descriptor.cell_owner_vector[primal_node_id] = primal_node_owner[primal_node_id]; + } + } + + { + std::vector<int> face_cell_owner(dual_descriptor.face_number_vector.size()); + std::fill(std::begin(face_cell_owner), std::end(face_cell_owner), parallel::size()); + + for (size_t i_cell = 0; i_cell < dual_descriptor.cell_to_face_vector.size(); ++i_cell) { + const auto& cell_face_list = dual_descriptor.cell_to_face_vector[i_cell]; + for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) { + const size_t face_id = cell_face_list[i_face]; + face_cell_owner[face_id] = std::min(face_cell_owner[face_id], dual_descriptor.cell_number_vector[i_cell]); + } + } + + dual_descriptor.face_owner_vector.resize(face_cell_owner.size()); + for (size_t i_face = 0; i_face < face_cell_owner.size(); ++i_face) { + dual_descriptor.face_owner_vector[i_face] = dual_descriptor.cell_owner_vector[face_cell_owner[i_face]]; + } + } + + m_connectivity = ConnectivityType::build(dual_descriptor); + + { + const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); + + NodeId dual_node_id = 0; + m_primal_node_to_dual_node_map = [&]() { + std::vector<std::pair<NodeId, NodeId>> primal_node_to_dual_node_vector; + 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) { + primal_node_to_dual_node_vector.push_back(std::make_pair(primal_node_id, dual_node_id++)); + } + } + return convert_to_array(primal_node_to_dual_node_vector); + }(); + + m_primal_cell_to_dual_node_map = [&]() { + CellIdToNodeIdMap primal_cell_to_dual_node_map{primal_number_of_cells}; + for (CellId primal_cell_id = 0; primal_cell_id < primal_cell_to_dual_node_map.size(); ++primal_cell_id) { + primal_cell_to_dual_node_map[primal_cell_id] = std::make_pair(primal_cell_id, dual_node_id++); + } + return primal_cell_to_dual_node_map; + }(); + + m_primal_face_to_dual_cell_map = [&]() { + FaceIdToCellIdMap primal_face_to_dual_cell_map{primal_connectivity.numberOfFaces()}; + for (size_t id = 0; id < primal_face_to_dual_cell_map.size(); ++id) { + const CellId dual_cell_id = id; + const FaceId primal_face_id = id; + + primal_face_to_dual_cell_map[id] = std::make_pair(primal_face_id, dual_cell_id); + } + return primal_face_to_dual_cell_map; + }(); + } + + m_mapper = + std::make_shared<PrimalToDual1DConnectivityDataMapper>(primal_connectivity, + dynamic_cast<const ConnectivityType&>(*m_connectivity), + m_primal_node_to_dual_node_map, + m_primal_cell_to_dual_node_map, + m_primal_face_to_dual_cell_map); +} + +Dual1DConnectivityBuilder::Dual1DConnectivityBuilder(const IConnectivity& connectivity) +{ + if (parallel::size() > 1) { + throw NotImplementedError("Construction of diamond dual mesh is not implemented in parallel"); + } + + if (connectivity.dimension() == 1) { + this->_buildConnectivityFrom(connectivity); + } else { + throw UnexpectedError("invalid connectivity dimension: " + std::to_string(connectivity.dimension())); + } +} diff --git a/src/mesh/Dual1DConnectivityBuilder.hpp b/src/mesh/Dual1DConnectivityBuilder.hpp new file mode 100644 index 000000000..dc534d0f4 --- /dev/null +++ b/src/mesh/Dual1DConnectivityBuilder.hpp @@ -0,0 +1,39 @@ +#ifndef DUAL_1D_CONNECTIVITY_BUILDER_HPP +#define DUAL_1D_CONNECTIVITY_BUILDER_HPP + +#include <mesh/ConnectivityBuilderBase.hpp> +#include <mesh/IPrimalToDualConnectivityDataMapper.hpp> +#include <mesh/ItemIdToItemIdMap.hpp> + +#include <memory> + +template <size_t> +class Connectivity; +class ConnectivityDescriptor; + +class Dual1DConnectivityBuilder : public ConnectivityBuilderBase +{ + NodeIdToNodeIdMap m_primal_node_to_dual_node_map; + CellIdToNodeIdMap m_primal_cell_to_dual_node_map; + FaceIdToCellIdMap m_primal_face_to_dual_cell_map; + + std::shared_ptr<IPrimalToDualConnectivityDataMapper> m_mapper; + + void _buildConnectivityDescriptor(const Connectivity<1>&, ConnectivityDescriptor&); + + void _buildConnectivityFrom(const IConnectivity&); + + friend class DualConnectivityManager; + Dual1DConnectivityBuilder(const IConnectivity&); + + public: + std::shared_ptr<IPrimalToDualConnectivityDataMapper> + mapper() const + { + return m_mapper; + } + + ~Dual1DConnectivityBuilder() = default; +}; + +#endif // DUAL_1D_CONNECTIVITY_BUILDER_HPP diff --git a/src/mesh/DualConnectivityManager.cpp b/src/mesh/DualConnectivityManager.cpp index 5a09a8214..fefc4bde5 100644 --- a/src/mesh/DualConnectivityManager.cpp +++ b/src/mesh/DualConnectivityManager.cpp @@ -2,9 +2,11 @@ #include <mesh/Connectivity.hpp> #include <mesh/DiamondDualConnectivityBuilder.hpp> +#include <mesh/Dual1DConnectivityBuilder.hpp> #include <mesh/DualConnectivityManager.hpp> #include <mesh/MedianDualConnectivityBuilder.hpp> #include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp> +#include <mesh/PrimalToDual1DConnectivityDataMapper.hpp> #include <mesh/PrimalToMedianDualConnectivityDataMapper.hpp> #include <utils/Exceptions.hpp> @@ -78,11 +80,14 @@ DualConnectivityManager::_getDualConnectivityInfo(const DualMeshType& type, cons return i_connectivity->second; } else { switch (type) { - case DualMeshType::Median: { - return this->_buildDualConnectivity<MedianDualConnectivityBuilder>(key, connectivity); - } case DualMeshType::Diamond: { return this->_buildDualConnectivity<DiamondDualConnectivityBuilder>(key, connectivity); + } + case DualMeshType::Dual1D: { + return this->_buildDualConnectivity<Dual1DConnectivityBuilder>(key, connectivity); + } + case DualMeshType::Median: { + return this->_buildDualConnectivity<MedianDualConnectivityBuilder>(key, connectivity); } // LCOV_EXCL_START default: { @@ -93,18 +98,50 @@ DualConnectivityManager::_getDualConnectivityInfo(const DualMeshType& type, cons } } +std::shared_ptr<const Connectivity<1>> +DualConnectivityManager::getDual1DConnectivity(const Connectivity<1>& connectivity) +{ + return std::dynamic_pointer_cast<const Connectivity<1>>( + this->_getDualConnectivityInfo(DualMeshType::Dual1D, connectivity).dualConnectivity()); +} + +std::shared_ptr<const PrimalToDual1DConnectivityDataMapper> +DualConnectivityManager::getPrimalToDual1DConnectivityDataMapper(const Connectivity<1>& connectivity) +{ + auto i_data_mapper = + this->_getDualConnectivityInfo(DualMeshType::Dual1D, connectivity).connectivityToDualConnectivityDataMapper(); + auto dual_1d_data_mapper = std::dynamic_pointer_cast<const PrimalToDual1DConnectivityDataMapper>(i_data_mapper); + + if (dual_1d_data_mapper.use_count() > 0) { + return dual_1d_data_mapper; + } else { + // LCOV_EXCL_START + throw UnexpectedError("invalid connectivity data mapper type"); + // LCOV_EXCL_STOP + } +} + template <size_t Dimension> std::shared_ptr<const Connectivity<Dimension>> DualConnectivityManager::getDiamondDualConnectivity(const Connectivity<Dimension>& connectivity) { + auto dual_mesh_type = (connectivity.dimension() == 1) ? DualMeshType::Dual1D : DualMeshType::Diamond; + return std::dynamic_pointer_cast<const Connectivity<Dimension>>( - this->_getDualConnectivityInfo(DualMeshType::Diamond, connectivity).dualConnectivity()); + this->_getDualConnectivityInfo(dual_mesh_type, connectivity).dualConnectivity()); +} + +std::shared_ptr<const PrimalToDual1DConnectivityDataMapper> +DualConnectivityManager::getPrimalToDiamondDualConnectivityDataMapper(const Connectivity<1>& connectivity) +{ + return this->getPrimalToDual1DConnectivityDataMapper(connectivity); } template <size_t Dimension> std::shared_ptr<const PrimalToDiamondDualConnectivityDataMapper<Dimension>> DualConnectivityManager::getPrimalToDiamondDualConnectivityDataMapper(const Connectivity<Dimension>& connectivity) { + static_assert((Dimension == 2) or (Dimension == 3)); auto i_data_mapper = this->_getDualConnectivityInfo(DualMeshType::Diamond, connectivity).connectivityToDualConnectivityDataMapper(); auto diamond_data_mapper = @@ -123,8 +160,16 @@ template <size_t Dimension> std::shared_ptr<const Connectivity<Dimension>> DualConnectivityManager::getMedianDualConnectivity(const Connectivity<Dimension>& connectivity) { + auto dual_mesh_type = (connectivity.dimension() == 1) ? DualMeshType::Dual1D : DualMeshType::Median; + return std::dynamic_pointer_cast<const Connectivity<Dimension>>( - this->_getDualConnectivityInfo(DualMeshType::Median, connectivity).dualConnectivity()); + this->_getDualConnectivityInfo(dual_mesh_type, connectivity).dualConnectivity()); +} + +std::shared_ptr<const PrimalToDual1DConnectivityDataMapper> +DualConnectivityManager::getPrimalToMedianDualConnectivityDataMapper(const Connectivity<1>& connectivity) +{ + return this->getPrimalToDual1DConnectivityDataMapper(connectivity); } template <size_t Dimension> @@ -154,9 +199,6 @@ template std::shared_ptr<const Connectivity<2>> DualConnectivityManager::getDiam template std::shared_ptr<const Connectivity<3>> DualConnectivityManager::getDiamondDualConnectivity( const Connectivity<3>& connectivity); -template std::shared_ptr<const PrimalToDiamondDualConnectivityDataMapper<1>> -DualConnectivityManager::getPrimalToDiamondDualConnectivityDataMapper(const Connectivity<1>&); - template std::shared_ptr<const PrimalToDiamondDualConnectivityDataMapper<2>> DualConnectivityManager::getPrimalToDiamondDualConnectivityDataMapper(const Connectivity<2>&); @@ -172,9 +214,6 @@ template std::shared_ptr<const Connectivity<2>> DualConnectivityManager::getMedi template std::shared_ptr<const Connectivity<3>> DualConnectivityManager::getMedianDualConnectivity( const Connectivity<3>& connectivity); -template std::shared_ptr<const PrimalToMedianDualConnectivityDataMapper<1>> -DualConnectivityManager::getPrimalToMedianDualConnectivityDataMapper(const Connectivity<1>&); - template std::shared_ptr<const PrimalToMedianDualConnectivityDataMapper<2>> DualConnectivityManager::getPrimalToMedianDualConnectivityDataMapper(const Connectivity<2>&); diff --git a/src/mesh/DualConnectivityManager.hpp b/src/mesh/DualConnectivityManager.hpp index a45e2226d..6389c5d45 100644 --- a/src/mesh/DualConnectivityManager.hpp +++ b/src/mesh/DualConnectivityManager.hpp @@ -17,6 +17,8 @@ class PrimalToDiamondDualConnectivityDataMapper; template <size_t Dimension> class PrimalToMedianDualConnectivityDataMapper; +class PrimalToDual1DConnectivityDataMapper; + class DualConnectivityManager { private: @@ -97,9 +99,16 @@ class DualConnectivityManager void deleteConnectivity(const IConnectivity*); + std::shared_ptr<const Connectivity<1>> getDual1DConnectivity(const Connectivity<1>&); + + std::shared_ptr<const PrimalToDual1DConnectivityDataMapper> getPrimalToDual1DConnectivityDataMapper( + const Connectivity<1>&); + template <size_t Dimension> std::shared_ptr<const Connectivity<Dimension>> getMedianDualConnectivity(const Connectivity<Dimension>&); + std::shared_ptr<const PrimalToDual1DConnectivityDataMapper> // special 1D case + getPrimalToMedianDualConnectivityDataMapper(const Connectivity<1>&); template <size_t Dimension> std::shared_ptr<const PrimalToMedianDualConnectivityDataMapper<Dimension>> getPrimalToMedianDualConnectivityDataMapper(const Connectivity<Dimension>&); @@ -107,6 +116,8 @@ class DualConnectivityManager template <size_t Dimension> std::shared_ptr<const Connectivity<Dimension>> getDiamondDualConnectivity(const Connectivity<Dimension>&); + std::shared_ptr<const PrimalToDual1DConnectivityDataMapper> // special 1D case + getPrimalToDiamondDualConnectivityDataMapper(const Connectivity<1>&); template <size_t Dimension> std::shared_ptr<const PrimalToDiamondDualConnectivityDataMapper<Dimension>> getPrimalToDiamondDualConnectivityDataMapper(const Connectivity<Dimension>&); diff --git a/src/mesh/DualMeshType.hpp b/src/mesh/DualMeshType.hpp index 16c370c7e..967d6a453 100644 --- a/src/mesh/DualMeshType.hpp +++ b/src/mesh/DualMeshType.hpp @@ -9,6 +9,7 @@ enum class DualMeshType { Diamond, + Dual1D, Median }; @@ -17,12 +18,15 @@ std::string name(DualMeshType type) { switch (type) { - case DualMeshType::Median: { - return "median"; - } case DualMeshType::Diamond: { return "diamond"; } + case DualMeshType::Dual1D: { + return "dual 1d"; + } + case DualMeshType::Median: { + return "median"; + } default: { throw UnexpectedError("unexpected dual mesh type"); } diff --git a/src/mesh/MedianDualConnectivityBuilder.cpp b/src/mesh/MedianDualConnectivityBuilder.cpp index ee345aa72..e260e2e15 100644 --- a/src/mesh/MedianDualConnectivityBuilder.cpp +++ b/src/mesh/MedianDualConnectivityBuilder.cpp @@ -392,13 +392,9 @@ MedianDualConnectivityBuilder::_buildConnectivityFrom<2>(const IConnectivity& i_ MedianDualConnectivityBuilder::MedianDualConnectivityBuilder(const IConnectivity& connectivity) { if (parallel::size() > 1) { - throw NotImplementedError("Construction of diamond dual mesh is not implemented in parallel"); + throw NotImplementedError("Construction of median dual mesh is not implemented in parallel"); } switch (connectivity.dimension()) { - case 1: { - throw NotImplementedError("must be treated as a special common case with diamond dual connectivities"); - break; - } case 2: { this->_buildConnectivityFrom<2>(connectivity); break; diff --git a/src/mesh/PrimalToDiamondDualConnectivityDataMapper.hpp b/src/mesh/PrimalToDiamondDualConnectivityDataMapper.hpp index a8e7d0e1a..b80c26d18 100644 --- a/src/mesh/PrimalToDiamondDualConnectivityDataMapper.hpp +++ b/src/mesh/PrimalToDiamondDualConnectivityDataMapper.hpp @@ -11,6 +11,9 @@ template <size_t Dimension> class PrimalToDiamondDualConnectivityDataMapper : public IPrimalToDualConnectivityDataMapper { + static_assert(Dimension == 2 or Dimension == 3, + "primal to diamond dual connectivity mapper is defined in dimension 2 or 3"); + private: const IConnectivity* m_primal_connectivity; const IConnectivity* m_dual_connectivity; diff --git a/src/mesh/PrimalToDual1DConnectivityDataMapper.hpp b/src/mesh/PrimalToDual1DConnectivityDataMapper.hpp new file mode 100644 index 000000000..1cad07fb1 --- /dev/null +++ b/src/mesh/PrimalToDual1DConnectivityDataMapper.hpp @@ -0,0 +1,140 @@ +#ifndef PRIMAL_TO_DUAL_1D_CONNECTIVITY_DATA_MAPPER_HPP +#define PRIMAL_TO_DUAL_1D_CONNECTIVITY_DATA_MAPPER_HPP + +#include <mesh/Connectivity.hpp> +#include <mesh/IPrimalToDualConnectivityDataMapper.hpp> +#include <mesh/ItemIdToItemIdMap.hpp> +#include <mesh/ItemValue.hpp> +#include <utils/Array.hpp> +#include <utils/PugsAssert.hpp> + +class PrimalToDual1DConnectivityDataMapper : public IPrimalToDualConnectivityDataMapper +{ + private: + const IConnectivity* m_primal_connectivity; + const IConnectivity* m_dual_connectivity; + + ConstNodeIdToNodeIdMap m_primal_node_to_dual_node_map; + ConstCellIdToNodeIdMap m_primal_cell_to_dual_node_map; + ConstFaceIdToCellIdMap m_primal_face_to_dual_cell_map; + + public: + template <typename OriginDataType1, typename OriginDataType2, typename DestinationDataType> + void + toDualNode(const NodeValue<OriginDataType1>& primal_node_value, + const CellValue<OriginDataType2>& primal_cell_value, + const NodeValue<DestinationDataType>& dual_node_value) const + { + static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant"); + static_assert(std::is_same_v<std::remove_const_t<OriginDataType1>, DestinationDataType>, "incompatible types"); + static_assert(std::is_same_v<std::remove_const_t<OriginDataType2>, DestinationDataType>, "incompatible types"); + + Assert(m_primal_connectivity == primal_cell_value.connectivity_ptr().get(), + "unexpected connectivity for primal CellValue"); + Assert(m_primal_connectivity == primal_node_value.connectivity_ptr().get(), + "unexpected connectivity for primal NodeValue"); + Assert(m_dual_connectivity == dual_node_value.connectivity_ptr().get(), + "unexpected connectivity for dual NodeValue"); + + parallel_for( + m_primal_node_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) { + const auto [primal_node_id, dual_node_id] = m_primal_node_to_dual_node_map[i]; + + dual_node_value[dual_node_id] = primal_node_value[primal_node_id]; + }); + + parallel_for( + m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) { + const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i]; + dual_node_value[dual_node_id] = primal_cell_value[primal_cell_id]; + }); + } + + template <typename OriginDataType, typename DestinationDataType1, typename DestinationDataType2> + void + fromDualNode(const NodeValue<OriginDataType>& dual_node_value, + const NodeValue<DestinationDataType1>& primal_node_value, + const CellValue<DestinationDataType2>& primal_cell_value) const + { + static_assert(not std::is_const_v<DestinationDataType1>, "destination data type must not be constant"); + static_assert(not std::is_const_v<DestinationDataType2>, "destination data type must not be constant"); + static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType1>, "incompatible types"); + static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType2>, "incompatible types"); + + Assert(m_primal_connectivity == primal_cell_value.connectivity_ptr().get(), + "unexpected connectivity for primal CellValue"); + Assert(m_primal_connectivity == primal_node_value.connectivity_ptr().get(), + "unexpected connectivity for primal NodeValue"); + Assert(m_dual_connectivity == dual_node_value.connectivity_ptr().get(), + "unexpected connectivity for dual NodeValue"); + + parallel_for( + m_primal_node_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) { + const auto [primal_node_id, dual_node_id] = m_primal_node_to_dual_node_map[i]; + + primal_node_value[primal_node_id] = dual_node_value[dual_node_id]; + }); + + parallel_for( + m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) { + const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i]; + primal_cell_value[primal_cell_id] = dual_node_value[dual_node_id]; + }); + } + + template <typename OriginDataType, typename DestinationDataType> + void + toDualCell(const FaceValue<OriginDataType>& primal_face_value, + const CellValue<DestinationDataType>& dual_cell_value) const + { + static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant"); + static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType>, "incompatible types"); + + Assert(m_primal_connectivity == primal_face_value.connectivity_ptr().get(), + "unexpected connectivity for primal FaceValue"); + Assert(m_dual_connectivity == dual_cell_value.connectivity_ptr().get(), + "unexpected connectivity for dual CellValue"); + + parallel_for( + m_primal_face_to_dual_cell_map.size(), PUGS_LAMBDA(size_t i) { + const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i]; + + dual_cell_value[dual_cell_id] = primal_face_value[primal_face_id]; + }); + } + + template <typename OriginDataType, typename DestinationDataType> + void + fromDualCell(const CellValue<DestinationDataType>& dual_cell_value, + const FaceValue<OriginDataType>& primal_face_value) const + { + static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant"); + static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType>, "incompatible types"); + + Assert(m_primal_connectivity == primal_face_value.connectivity_ptr().get(), + "unexpected connectivity for primal FaceValue"); + Assert(m_dual_connectivity == dual_cell_value.connectivity_ptr().get(), + "unexpected connectivity for dual CellValue"); + + parallel_for( + m_primal_face_to_dual_cell_map.size(), PUGS_LAMBDA(size_t i) { + const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i]; + + primal_face_value[primal_face_id] = dual_cell_value[dual_cell_id]; + }); + } + + PrimalToDual1DConnectivityDataMapper(const Connectivity<1>& primal_connectivity, + const Connectivity<1>& dual_connectivity, + const ConstNodeIdToNodeIdMap& primal_node_to_dual_node_map, + const ConstCellIdToNodeIdMap& primal_cell_to_dual_node_map, + const ConstFaceIdToCellIdMap& primal_face_to_dual_cell_map) + : m_primal_connectivity{&primal_connectivity}, + m_dual_connectivity{&dual_connectivity}, + m_primal_node_to_dual_node_map{primal_node_to_dual_node_map}, + m_primal_cell_to_dual_node_map{primal_cell_to_dual_node_map}, + m_primal_face_to_dual_cell_map{primal_face_to_dual_cell_map} + {} +}; + +#endif // PRIMAL_TO_DUAL_1D_CONNECTIVITY_DATA_MAPPER_HPP diff --git a/src/mesh/PrimalToMedianDualConnectivityDataMapper.hpp b/src/mesh/PrimalToMedianDualConnectivityDataMapper.hpp index c511930eb..15b846ac6 100644 --- a/src/mesh/PrimalToMedianDualConnectivityDataMapper.hpp +++ b/src/mesh/PrimalToMedianDualConnectivityDataMapper.hpp @@ -11,6 +11,9 @@ template <size_t Dimension> class PrimalToMedianDualConnectivityDataMapper : public IPrimalToDualConnectivityDataMapper { + static_assert(Dimension == 2 or Dimension == 3, + "primal to median dual connectivity mapper is defined in dimension 2 or 3"); + private: const IConnectivity* m_primal_connectivity; const IConnectivity* m_dual_connectivity; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6edfaeb58..6985be9cb 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -103,6 +103,7 @@ add_executable (unit_tests test_ParseError.cpp test_PETScUtils.cpp test_PrimalToDiamondDualConnectivityDataMapper.cpp + test_PrimalToDual1DConnectivityDataMapper.cpp test_PrimalToMedianDualConnectivityDataMapper.cpp test_PrismGaussQuadrature.cpp test_PrismTransformation.cpp diff --git a/tests/test_PrimalToDiamondDualConnectivityDataMapper.cpp b/tests/test_PrimalToDiamondDualConnectivityDataMapper.cpp index ec628b753..d34c29d41 100644 --- a/tests/test_PrimalToDiamondDualConnectivityDataMapper.cpp +++ b/tests/test_PrimalToDiamondDualConnectivityDataMapper.cpp @@ -4,6 +4,7 @@ #include <MeshDataBaseForTests.hpp> #include <mesh/DualConnectivityManager.hpp> #include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp> +#include <mesh/PrimalToDual1DConnectivityDataMapper.hpp> #include <mesh/Connectivity.hpp> #include <mesh/Mesh.hpp> @@ -24,141 +25,19 @@ TEST_CASE("PrimalToDiamondDualConnectivityDataMapper", "[mesh]") std::shared_ptr p_diamond_dual_connectivity = DualConnectivityManager::instance().getDiamondDualConnectivity(primal_connectivity); - const ConnectivityType& dual_connectivity = *p_diamond_dual_connectivity; - - std::shared_ptr p_primal_to_dual_connectivity_mapper = - DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(primal_connectivity); - const auto& primal_to_dual_connectivity_mapper = *p_primal_to_dual_connectivity_mapper; - - SECTION("check data transfer between primal faces to dual cells") - { - FaceValue<size_t> primal_face_value{primal_connectivity}; - parallel_for( - primal_face_value.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) { primal_face_value[face_id] = face_id; }); - - CellValue<size_t> dual_from_primal_face_value{dual_connectivity}; - - primal_to_dual_connectivity_mapper.toDualCell(primal_face_value, dual_from_primal_face_value); - - { - bool correct_value = true; - for (CellId cell_id = 0; cell_id < dual_connectivity.numberOfCells(); ++cell_id) { - correct_value &= (cell_id == dual_from_primal_face_value[cell_id]); - } - REQUIRE(correct_value); - } - - { - bool is_same = true; - FaceValue<size_t> primal_from_dual_cell_value{primal_connectivity}; - primal_to_dual_connectivity_mapper.fromDualCell(dual_from_primal_face_value, primal_from_dual_cell_value); - for (FaceId face_id = 0; face_id < primal_connectivity.numberOfFaces(); ++face_id) { - is_same &= (primal_from_dual_cell_value[face_id] == primal_face_value[face_id]); - } - REQUIRE(is_same); - } -#ifndef NDEBUG - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualCell(FaceValue<size_t>{dual_connectivity}, - dual_from_primal_face_value), - "unexpected connectivity for primal FaceValue"); - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualCell(primal_face_value, - CellValue<size_t>{primal_connectivity}), - "unexpected connectivity for dual CellValue"); + std::shared_ptr p_dual_1d_connectivity = + DualConnectivityManager::instance().getDual1DConnectivity(primal_connectivity); - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualCell(dual_from_primal_face_value, - FaceValue<size_t>{dual_connectivity}), - "unexpected connectivity for primal FaceValue"); - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualCell(CellValue<size_t>{primal_connectivity}, - primal_face_value), - "unexpected connectivity for dual CellValue"); + REQUIRE(p_diamond_dual_connectivity == p_dual_1d_connectivity); -#endif // NDEBUG - } - - SECTION("check data transfer between primal nodes and cells to dual nodes") - { - const auto& primal_cell_to_node_matrix = primal_connectivity.cellToNodeMatrix(); - CellValue<size_t> primal_cell_value{primal_connectivity}; - parallel_for( - primal_cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { - // 2 is the number of boundary nodes - primal_cell_value[cell_id] = primal_cell_to_node_matrix[cell_id].size() + 2 + cell_id; - }); - - const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); - NodeValue<size_t> primal_node_value{primal_connectivity}; - primal_node_value.fill(0); - size_t count = 0; - for (NodeId node_id = 0; node_id < primal_node_value.numberOfItems(); ++node_id) { - if (primal_node_to_cell_matrix[node_id].size() == 1) { - primal_node_value[node_id] = primal_node_to_cell_matrix[node_id].size() + count; - count++; - } - } - NodeValue<size_t> dual_from_primal_cell_and_node_value{dual_connectivity}; - - primal_to_dual_connectivity_mapper.toDualNode(primal_node_value, primal_cell_value, - dual_from_primal_cell_and_node_value); - - { - const auto& dual_node_to_cell_matrix = dual_connectivity.nodeToCellMatrix(); - bool correct_dual_value = true; - for (NodeId node_id = 0; node_id < dual_connectivity.numberOfNodes(); ++node_id) { - correct_dual_value &= - (dual_from_primal_cell_and_node_value[node_id] == dual_node_to_cell_matrix[node_id].size() + node_id); - - CHECK(dual_from_primal_cell_and_node_value[node_id] == dual_node_to_cell_matrix[node_id].size() + node_id); - } - REQUIRE(correct_dual_value); - } - - { - bool is_same = true; - NodeValue<size_t> primal_from_dual_node_value{primal_connectivity}; - CellValue<size_t> primal_from_dual_cell_value{primal_connectivity}; - primal_to_dual_connectivity_mapper.fromDualNode(dual_from_primal_cell_and_node_value, - primal_from_dual_node_value, primal_from_dual_cell_value); - for (CellId cell_id = 0; cell_id < primal_connectivity.numberOfCells(); ++cell_id) { - is_same &= (primal_cell_value[cell_id] == primal_from_dual_cell_value[cell_id]); - } - - for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) { - if (primal_node_to_cell_matrix[node_id].size() == 1) { - is_same &= (primal_node_value[node_id] == primal_from_dual_node_value[node_id]); - } - } - - REQUIRE(is_same); - } - -#ifndef NDEBUG - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualNode(NodeValue<size_t>{dual_connectivity}, - primal_cell_value, - dual_from_primal_cell_and_node_value), - "unexpected connectivity for primal NodeValue"); - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualNode(primal_node_value, - CellValue<size_t>{dual_connectivity}, - dual_from_primal_cell_and_node_value), - "unexpected connectivity for primal CellValue"); - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualNode(primal_node_value, primal_cell_value, - NodeValue<size_t>{primal_connectivity}), - "unexpected connectivity for dual NodeValue"); + std::shared_ptr p_primal_to_diamond_dual_connectivity_mapper = + DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(primal_connectivity); - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualNode(dual_from_primal_cell_and_node_value, - NodeValue<size_t>{dual_connectivity}, - primal_cell_value), - "unexpected connectivity for primal NodeValue"); - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualNode(dual_from_primal_cell_and_node_value, - primal_node_value, - CellValue<size_t>{dual_connectivity}), - "unexpected connectivity for primal CellValue"); - REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualNode(NodeValue<size_t>{primal_connectivity}, - primal_node_value, primal_cell_value), - "unexpected connectivity for dual NodeValue"); + std::shared_ptr p_primal_to_dual_1d_connectivity_mapper = + DualConnectivityManager::instance().getPrimalToDual1DConnectivityDataMapper(primal_connectivity); -#endif // NDEBUG - } + REQUIRE(p_primal_to_diamond_dual_connectivity_mapper == p_primal_to_dual_1d_connectivity_mapper); } SECTION("2D") diff --git a/tests/test_PrimalToDual1DConnectivityDataMapper.cpp b/tests/test_PrimalToDual1DConnectivityDataMapper.cpp new file mode 100644 index 000000000..fef9680c8 --- /dev/null +++ b/tests/test_PrimalToDual1DConnectivityDataMapper.cpp @@ -0,0 +1,160 @@ +#include <catch2/catch_test_macros.hpp> +#include <catch2/matchers/catch_matchers_all.hpp> + +#include <MeshDataBaseForTests.hpp> +#include <mesh/DualConnectivityManager.hpp> +#include <mesh/PrimalToDual1DConnectivityDataMapper.hpp> + +#include <mesh/Connectivity.hpp> +#include <mesh/Mesh.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("PrimalToDual1DConnectivityDataMapper", "[mesh]") +{ + constexpr static size_t Dimension = 1; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + std::shared_ptr<const MeshType> mesh = MeshDataBaseForTests::get().unordered1DMesh(); + const ConnectivityType& primal_connectivity = mesh->connectivity(); + + std::shared_ptr p_dual_1d_connectivity = + DualConnectivityManager::instance().getDual1DConnectivity(primal_connectivity); + const ConnectivityType& dual_connectivity = *p_dual_1d_connectivity; + + std::shared_ptr p_primal_to_dual_connectivity_mapper = + DualConnectivityManager::instance().getPrimalToDual1DConnectivityDataMapper(primal_connectivity); + const auto& primal_to_dual_connectivity_mapper = *p_primal_to_dual_connectivity_mapper; + + SECTION("check data transfer between primal faces to dual cells") + { + FaceValue<size_t> primal_face_value{primal_connectivity}; + parallel_for( + primal_face_value.numberOfItems(), PUGS_LAMBDA(const FaceId face_id) { primal_face_value[face_id] = face_id; }); + + CellValue<size_t> dual_from_primal_face_value{dual_connectivity}; + + primal_to_dual_connectivity_mapper.toDualCell(primal_face_value, dual_from_primal_face_value); + + { + bool correct_value = true; + for (CellId cell_id = 0; cell_id < dual_connectivity.numberOfCells(); ++cell_id) { + correct_value &= (cell_id == dual_from_primal_face_value[cell_id]); + } + REQUIRE(correct_value); + } + + { + bool is_same = true; + FaceValue<size_t> primal_from_dual_cell_value{primal_connectivity}; + primal_to_dual_connectivity_mapper.fromDualCell(dual_from_primal_face_value, primal_from_dual_cell_value); + for (FaceId face_id = 0; face_id < primal_connectivity.numberOfFaces(); ++face_id) { + is_same &= (primal_from_dual_cell_value[face_id] == primal_face_value[face_id]); + } + REQUIRE(is_same); + } + +#ifndef NDEBUG + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualCell(FaceValue<size_t>{dual_connectivity}, + dual_from_primal_face_value), + "unexpected connectivity for primal FaceValue"); + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualCell(primal_face_value, + CellValue<size_t>{primal_connectivity}), + "unexpected connectivity for dual CellValue"); + + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualCell(dual_from_primal_face_value, + FaceValue<size_t>{dual_connectivity}), + "unexpected connectivity for primal FaceValue"); + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualCell(CellValue<size_t>{primal_connectivity}, + primal_face_value), + "unexpected connectivity for dual CellValue"); + +#endif // NDEBUG + } + + SECTION("check data transfer between primal nodes and cells to dual nodes") + { + const auto& primal_cell_to_node_matrix = primal_connectivity.cellToNodeMatrix(); + CellValue<size_t> primal_cell_value{primal_connectivity}; + parallel_for( + primal_cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { + // 2 is the number of boundary nodes + primal_cell_value[cell_id] = primal_cell_to_node_matrix[cell_id].size() + 2 + cell_id; + }); + + const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix(); + NodeValue<size_t> primal_node_value{primal_connectivity}; + primal_node_value.fill(0); + size_t count = 0; + for (NodeId node_id = 0; node_id < primal_node_value.numberOfItems(); ++node_id) { + if (primal_node_to_cell_matrix[node_id].size() == 1) { + primal_node_value[node_id] = primal_node_to_cell_matrix[node_id].size() + count; + count++; + } + } + NodeValue<size_t> dual_from_primal_cell_and_node_value{dual_connectivity}; + + primal_to_dual_connectivity_mapper.toDualNode(primal_node_value, primal_cell_value, + dual_from_primal_cell_and_node_value); + + { + const auto& dual_node_to_cell_matrix = dual_connectivity.nodeToCellMatrix(); + bool correct_dual_value = true; + for (NodeId node_id = 0; node_id < dual_connectivity.numberOfNodes(); ++node_id) { + correct_dual_value &= + (dual_from_primal_cell_and_node_value[node_id] == dual_node_to_cell_matrix[node_id].size() + node_id); + + CHECK(dual_from_primal_cell_and_node_value[node_id] == dual_node_to_cell_matrix[node_id].size() + node_id); + } + REQUIRE(correct_dual_value); + } + + { + bool is_same = true; + NodeValue<size_t> primal_from_dual_node_value{primal_connectivity}; + CellValue<size_t> primal_from_dual_cell_value{primal_connectivity}; + primal_to_dual_connectivity_mapper.fromDualNode(dual_from_primal_cell_and_node_value, primal_from_dual_node_value, + primal_from_dual_cell_value); + for (CellId cell_id = 0; cell_id < primal_connectivity.numberOfCells(); ++cell_id) { + is_same &= (primal_cell_value[cell_id] == primal_from_dual_cell_value[cell_id]); + } + + for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) { + if (primal_node_to_cell_matrix[node_id].size() == 1) { + is_same &= (primal_node_value[node_id] == primal_from_dual_node_value[node_id]); + } + } + + REQUIRE(is_same); + } + +#ifndef NDEBUG + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualNode(NodeValue<size_t>{dual_connectivity}, + primal_cell_value, + dual_from_primal_cell_and_node_value), + "unexpected connectivity for primal NodeValue"); + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualNode(primal_node_value, + CellValue<size_t>{dual_connectivity}, + dual_from_primal_cell_and_node_value), + "unexpected connectivity for primal CellValue"); + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.toDualNode(primal_node_value, primal_cell_value, + NodeValue<size_t>{primal_connectivity}), + "unexpected connectivity for dual NodeValue"); + + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualNode(dual_from_primal_cell_and_node_value, + NodeValue<size_t>{dual_connectivity}, + primal_cell_value), + "unexpected connectivity for primal NodeValue"); + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualNode(dual_from_primal_cell_and_node_value, + primal_node_value, + CellValue<size_t>{dual_connectivity}), + "unexpected connectivity for primal CellValue"); + REQUIRE_THROWS_WITH(primal_to_dual_connectivity_mapper.fromDualNode(NodeValue<size_t>{primal_connectivity}, + primal_node_value, primal_cell_value), + "unexpected connectivity for dual NodeValue"); + +#endif // NDEBUG + } +} diff --git a/tests/test_PrimalToMedianDualConnectivityDataMapper.cpp b/tests/test_PrimalToMedianDualConnectivityDataMapper.cpp index 2b75b8840..437557704 100644 --- a/tests/test_PrimalToMedianDualConnectivityDataMapper.cpp +++ b/tests/test_PrimalToMedianDualConnectivityDataMapper.cpp @@ -3,6 +3,7 @@ #include <MeshDataBaseForTests.hpp> #include <mesh/DualConnectivityManager.hpp> +#include <mesh/PrimalToDual1DConnectivityDataMapper.hpp> #include <mesh/PrimalToMedianDualConnectivityDataMapper.hpp> #include <mesh/Connectivity.hpp> @@ -12,6 +13,33 @@ TEST_CASE("PrimalToMedianDualConnectivityDataMapper", "[mesh]") { + SECTION("1D") + { + constexpr static size_t Dimension = 1; + + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + + std::shared_ptr<const MeshType> mesh = MeshDataBaseForTests::get().unordered1DMesh(); + const ConnectivityType& primal_connectivity = mesh->connectivity(); + + std::shared_ptr p_median_dual_connectivity = + DualConnectivityManager::instance().getMedianDualConnectivity(primal_connectivity); + + std::shared_ptr p_dual_1d_connectivity = + DualConnectivityManager::instance().getDual1DConnectivity(primal_connectivity); + + REQUIRE(p_median_dual_connectivity == p_dual_1d_connectivity); + + std::shared_ptr p_primal_to_median_dual_connectivity_mapper = + DualConnectivityManager::instance().getPrimalToMedianDualConnectivityDataMapper(primal_connectivity); + + std::shared_ptr p_primal_to_dual_1d_connectivity_mapper = + DualConnectivityManager::instance().getPrimalToDual1DConnectivityDataMapper(primal_connectivity); + + REQUIRE(p_primal_to_median_dual_connectivity_mapper == p_primal_to_dual_1d_connectivity_mapper); + } + SECTION("2D") { constexpr static size_t Dimension = 2; -- GitLab