#ifndef CONNECTIVITY_HPP #define CONNECTIVITY_HPP #include <algebra/TinyVector.hpp> #include <mesh/CellType.hpp> #include <mesh/ConnectivityComputer.hpp> #include <mesh/ConnectivityMatrix.hpp> #include <mesh/IConnectivity.hpp> #include <mesh/ItemToItemMatrix.hpp> #include <mesh/ItemType.hpp> #include <mesh/ItemValue.hpp> #include <mesh/RefId.hpp> #include <mesh/RefItemList.hpp> #include <mesh/SubItemValuePerItem.hpp> #include <utils/CRSGraph.hpp> #include <utils/Exceptions.hpp> #include <utils/PugsAssert.hpp> #include <utils/PugsMacros.hpp> #include <utils/PugsTraits.hpp> #include <utils/PugsUtils.hpp> #include <algorithm> #include <iostream> #include <set> #include <vector> class ConnectivityDescriptor; template <size_t Dim> class Connectivity final : public IConnectivity { public: PUGS_INLINE static std::shared_ptr<Connectivity<Dim>> build(const ConnectivityDescriptor&); private: constexpr static auto& itemTId = ItemTypeId<Dim>::itemTId; public: static constexpr size_t Dimension = Dim; PUGS_INLINE size_t dimension() const final { return Dimension; } private: ConnectivityMatrix m_item_to_item_matrix[Dimension + 1][Dimension + 1]; WeakCellValue<const CellType> m_cell_type; WeakCellValue<const int> m_cell_global_index; WeakCellValue<const int> m_cell_number; WeakFaceValue<const int> m_face_number; WeakEdgeValue<const int> m_edge_number; WeakNodeValue<const int> m_node_number; WeakCellValue<const int> m_cell_owner; WeakCellValue<const bool> m_cell_is_owned; WeakFaceValue<const int> m_face_owner; WeakFaceValue<const bool> m_face_is_owned; WeakEdgeValue<const int> m_edge_owner; WeakEdgeValue<const bool> m_edge_is_owned; WeakNodeValue<const int> m_node_owner; WeakNodeValue<const bool> m_node_is_owned; WeakFaceValuePerCell<const bool> m_cell_face_is_reversed; WeakEdgeValuePerFace<const bool> m_face_edge_is_reversed; WeakNodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes; WeakEdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges; WeakFaceValuePerCell<const unsigned short> m_cell_local_numbers_in_their_faces; WeakCellValuePerFace<const unsigned short> m_face_local_numbers_in_their_cells; WeakEdgeValuePerFace<const unsigned short> m_face_local_numbers_in_their_edges; WeakNodeValuePerFace<const unsigned short> m_face_local_numbers_in_their_nodes; WeakCellValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_cells; WeakFaceValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_faces; WeakNodeValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_nodes; WeakCellValuePerNode<const unsigned short> m_node_local_numbers_in_their_cells; WeakEdgeValuePerNode<const unsigned short> m_node_local_numbers_in_their_edges; WeakFaceValuePerNode<const unsigned short> m_node_local_numbers_in_their_faces; ConnectivityComputer m_connectivity_computer; std::vector<RefCellList> m_ref_cell_list_vector; std::vector<RefFaceList> m_ref_face_list_vector; std::vector<RefEdgeList> m_ref_edge_list_vector; std::vector<RefNodeList> m_ref_node_list_vector; WeakCellValue<const double> m_inv_cell_nb_nodes; void _computeCellFaceAndFaceNodeConnectivities(); template <typename SubItemValuePerItemType> PUGS_INLINE const SubItemValuePerItemType& _lazzyBuildItemNumberInTheirChild(const SubItemValuePerItemType& sub_item_value_per_item) const { using ReversedItemOfItem = typename SubItemValuePerItemType::ItemOfItemType::Reversed; if (not sub_item_value_per_item.isBuilt()) { const_cast<SubItemValuePerItemType&>(sub_item_value_per_item) = m_connectivity_computer.computeLocalItemNumberInChildItem<ReversedItemOfItem>(*this); } return sub_item_value_per_item; } friend class ConnectivityComputer; PUGS_INLINE const ConnectivityMatrix& _getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const { const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)]; if (not connectivity_matrix.isBuilt()) { const_cast<ConnectivityMatrix&>(connectivity_matrix) = m_connectivity_computer.computeConnectivityMatrix(*this, item_type_0, item_type_1); } return connectivity_matrix; } public: PUGS_INLINE const auto& cellType() const { return m_cell_type; } PUGS_INLINE const auto& cellNumber() const { return m_cell_number; } PUGS_INLINE const auto& faceNumber() const { return m_face_number; } PUGS_INLINE const auto& edgeNumber() const { return m_edge_number; } PUGS_INLINE const auto& nodeNumber() const { return m_node_number; } template <ItemType item_type> PUGS_INLINE const auto& number() const { if constexpr (item_type == ItemType::cell) { return m_cell_number; } else if constexpr (item_type == ItemType::face) { return m_face_number; } else if constexpr (item_type == ItemType::edge) { return m_edge_number; } else if constexpr (item_type == ItemType::node) { return m_node_number; } else { static_assert(is_false_item_type_v<item_type>, "unknown ItemType"); } } PUGS_INLINE const auto& cellOwner() const { return m_cell_owner; } PUGS_INLINE const auto& faceOwner() const { return m_face_owner; } PUGS_INLINE const auto& edgeOwner() const { return m_edge_owner; } PUGS_INLINE const auto& nodeOwner() const { return m_node_owner; } template <ItemType item_type> PUGS_INLINE const auto& owner() const { if constexpr (item_type == ItemType::cell) { return m_cell_owner; } else if constexpr (item_type == ItemType::face) { return m_face_owner; } else if constexpr (item_type == ItemType::edge) { return m_edge_owner; } else if constexpr (item_type == ItemType::node) { return m_node_owner; } else { static_assert(is_false_item_type_v<item_type>, "unknown ItemType"); } } PUGS_INLINE const auto& cellIsOwned() const { return m_cell_is_owned; } PUGS_INLINE const auto& faceIsOwned() const { return m_face_is_owned; } PUGS_INLINE const auto& edgeIsOwned() const { return m_edge_is_owned; } PUGS_INLINE const auto& nodeIsOwned() const { return m_node_is_owned; } template <ItemType item_type> PUGS_INLINE const auto& isOwned() const { if constexpr (item_type == ItemType::cell) { return m_cell_is_owned; } else if constexpr (item_type == ItemType::face) { return m_face_is_owned; } else if constexpr (item_type == ItemType::edge) { return m_edge_is_owned; } else if constexpr (item_type == ItemType::node) { return m_node_is_owned; } else { static_assert(is_false_item_type_v<item_type>, "unknown ItemType"); } } PUGS_INLINE const bool& isConnectivityMatrixBuilt(ItemType item_type_0, ItemType item_type_1) const { const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)]; return connectivity_matrix.isBuilt(); } template <ItemType source_item_type, ItemType target_item_type> PUGS_INLINE auto getItemToItemMatrix() const { return ItemToItemMatrix<source_item_type, target_item_type>(_getMatrix(source_item_type, target_item_type)); } PUGS_INLINE auto cellToFaceMatrix() const { return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>(); } PUGS_INLINE auto cellToEdgeMatrix() const { return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>(); } PUGS_INLINE auto cellToNodeMatrix() const { return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>(); } PUGS_INLINE auto faceToCellMatrix() const { return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>(); } PUGS_INLINE auto faceToEdgeMatrix() const { return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>(); } PUGS_INLINE auto faceToNodeMatrix() const { return this->template getItemToItemMatrix<ItemType::face, ItemType::node>(); } PUGS_INLINE auto edgeToCellMatrix() const { return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>(); } PUGS_INLINE auto edgeToFaceMatrix() const { return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>(); } PUGS_INLINE auto edgeToNodeMatrix() const { return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>(); } PUGS_INLINE auto nodeToCellMatrix() const { return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>(); } PUGS_INLINE auto nodeToFaceMatrix() const { return this->template getItemToItemMatrix<ItemType::node, ItemType::face>(); } PUGS_INLINE auto nodeToEdgeMatrix() const { return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>(); } PUGS_INLINE const auto& cellFaceIsReversed() const { static_assert(Dimension > 1, "reversed faces makes no sense in dimension 1"); return m_cell_face_is_reversed; } PUGS_INLINE const auto& faceEdgeIsReversed() const { static_assert(Dimension > 2, "reversed edges makes no sense in dimension 1 or 2"); return m_face_edge_is_reversed; } PUGS_INLINE const auto& cellLocalNumbersInTheirNodes() const { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes); } PUGS_INLINE const auto& cellLocalNumbersInTheirEdges() const { if constexpr (Dimension > 2) { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_edges); } else { return cellLocalNumbersInTheirFaces(); } } PUGS_INLINE const auto& cellLocalNumbersInTheirFaces() const { if constexpr (Dimension > 1) { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_faces); } else { return cellLocalNumbersInTheirNodes(); } } PUGS_INLINE const auto& faceLocalNumbersInTheirCells() const { if constexpr (Dimension > 1) { return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_cells); } else { return nodeLocalNumbersInTheirCells(); } } PUGS_INLINE const auto& faceLocalNumbersInTheirEdges() const { static_assert(Dimension > 2, "this function has no meaning in 1d or 2d"); return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_edges); } PUGS_INLINE const auto& faceLocalNumbersInTheirNodes() const { static_assert(Dimension > 1, "this function has no meaning in 1d"); return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_nodes); } PUGS_INLINE const auto& edgeLocalNumbersInTheirCells() const { if constexpr (Dimension > 2) { return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_cells); } else { return faceLocalNumbersInTheirCells(); } } PUGS_INLINE const auto& edgeLocalNumbersInTheirFaces() const { static_assert(Dimension > 2, "this function has no meaning in 1d or 2d"); return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_faces); } PUGS_INLINE const auto& edgeLocalNumbersInTheirNodes() const { static_assert(Dimension > 2, "this function has no meaning in 1d or 2d"); return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_nodes); } PUGS_INLINE const auto& nodeLocalNumbersInTheirCells() const { return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells); } PUGS_INLINE const auto& nodeLocalNumbersInTheirEdges() const { static_assert(Dimension > 2, "this function has no meaning in 1d or 2d"); return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_edges); } PUGS_INLINE const auto& nodeLocalNumbersInTheirFaces() const { static_assert(Dimension > 1, "this function has no meaning in 1d"); return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces); } template <ItemType item_type> size_t numberOfRefItemList() const { if constexpr (item_type == ItemType::cell) { return m_ref_cell_list_vector.size(); } else if constexpr (item_type == ItemType::face) { return m_ref_face_list_vector.size(); } else if constexpr (item_type == ItemType::edge) { return m_ref_edge_list_vector.size(); } else if constexpr (item_type == ItemType::node) { return m_ref_node_list_vector.size(); } else { static_assert(is_false_item_type_v<item_type>, "Unexpected item type"); } } template <ItemType item_type> const RefItemList<item_type>& refItemList(size_t i) const { if constexpr (item_type == ItemType::cell) { return m_ref_cell_list_vector[i]; } else if constexpr (item_type == ItemType::face) { return m_ref_face_list_vector[i]; } else if constexpr (item_type == ItemType::edge) { return m_ref_edge_list_vector[i]; } else if constexpr (item_type == ItemType::node) { return m_ref_node_list_vector[i]; } else { static_assert(is_false_item_type_v<item_type>, "Unexpected item type"); } } template <ItemType item_type> void addRefItemList(const RefItemList<item_type>& ref_item_list) { if constexpr (item_type == ItemType::cell) { m_ref_cell_list_vector.push_back(ref_item_list); } else if constexpr (item_type == ItemType::face) { m_ref_face_list_vector.push_back(ref_item_list); } else if constexpr (item_type == ItemType::edge) { m_ref_edge_list_vector.push_back(ref_item_list); } else if constexpr (item_type == ItemType::node) { m_ref_node_list_vector.push_back(ref_item_list); } else { static_assert(is_false_item_type_v<item_type>, "Unexpected item type"); } } PUGS_INLINE CRSGraph cellToCellGraph() const { std::vector<std::set<int>> cell_cells(this->numberOfCells()); const auto& face_to_cell_matrix = this->faceToCellMatrix(); for (FaceId l = 0; l < this->numberOfFaces(); ++l) { const auto& face_to_cell = face_to_cell_matrix[l]; if (face_to_cell.size() > 1) { const CellId cell_0 = face_to_cell[0]; const CellId cell_1 = face_to_cell[1]; cell_cells[cell_0].insert(cell_1); cell_cells[cell_1].insert(cell_0); } } Array<int> entries(this->numberOfCells() + 1); entries[0] = 0; for (size_t j = 0; j < this->numberOfCells(); ++j) { entries[j + 1] = entries[j] + cell_cells[j].size(); } Array<int> neighbors(entries[this->numberOfCells()]); { size_t k = 0; for (size_t j = 0; j < this->numberOfCells(); ++j) { for (CellId cell_id : cell_cells[j]) { neighbors[k] = m_cell_global_index[cell_id]; ++k; } } } return CRSGraph(entries, neighbors); } PUGS_INLINE size_t numberOfNodes() const final { const auto& node_to_cell_matrix = this->_getMatrix(ItemType::node, ItemType::cell); return node_to_cell_matrix.numRows(); } PUGS_INLINE size_t numberOfEdges() const final { if constexpr (Dimension == 1) { return this->numberOfNodes(); } else if constexpr (Dimension == 2) { return this->numberOfFaces(); } else { const auto& edge_to_node_matrix = this->_getMatrix(ItemType::edge, ItemType::node); return edge_to_node_matrix.numRows(); } } PUGS_INLINE size_t numberOfFaces() const final { const auto& face_to_node_matrix = this->_getMatrix(ItemType::face, ItemType::cell); return face_to_node_matrix.numRows(); } PUGS_INLINE size_t numberOfCells() const final { const auto& cell_to_node_matrix = this->_getMatrix(ItemType::cell, ItemType::node); return cell_to_node_matrix.numRows(); } template <ItemType item_type> PUGS_INLINE size_t numberOf() const { if constexpr (item_type == ItemType::cell) { return this->numberOfCells(); } else if constexpr (item_type == ItemType::face) { return this->numberOfFaces(); } else if constexpr (item_type == ItemType::edge) { return this->numberOfEdges(); } else if constexpr (item_type == ItemType::node) { return this->numberOfNodes(); } else { static_assert(is_false_item_type_v<item_type>, "unknown ItemType"); } } CellValue<const double> invCellNbNodes() const { if (not m_inv_cell_nb_nodes.isBuilt()) { const auto& cell_to_node_matrix = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)]; WeakCellValue<double> inv_cell_nb_nodes(*this); parallel_for( this->numberOfCells(), PUGS_LAMBDA(CellId j) { const auto& cell_nodes = cell_to_node_matrix.rowConst(j); inv_cell_nb_nodes[j] = 1. / cell_nodes.length; }); const_cast<WeakCellValue<const double>&>(m_inv_cell_nb_nodes) = inv_cell_nb_nodes; } return m_inv_cell_nb_nodes; } Connectivity(const Connectivity&) = delete; private: Connectivity(); void _buildFrom(const ConnectivityDescriptor& descriptor); public: ~Connectivity() = default; }; template <size_t Dim> PUGS_INLINE std::shared_ptr<Connectivity<Dim>> Connectivity<Dim>::build(const ConnectivityDescriptor& descriptor) { std::shared_ptr<Connectivity<Dim>> connectivity_ptr(new Connectivity<Dim>); connectivity_ptr->_buildFrom(descriptor); return connectivity_ptr; } using Connectivity3D = Connectivity<3>; using Connectivity2D = Connectivity<2>; using Connectivity1D = Connectivity<1>; #endif // CONNECTIVITY_HPP