#ifndef CONNECTIVITY_HPP #define CONNECTIVITY_HPP #include <PastisMacros.hpp> #include <PastisAssert.hpp> #include <PastisOStream.hpp> #include <PastisUtils.hpp> #include <PastisTraits.hpp> #include <TinyVector.hpp> #include <ItemValue.hpp> #include <IConnectivity.hpp> #include <ConnectivityMatrix.hpp> #include <ItemToItemMatrix.hpp> #include <SubItemValuePerItem.hpp> #include <ConnectivityComputer.hpp> #include <vector> #include <algorithm> #include <CellType.hpp> #include <CSRGraph.hpp> #include <RefId.hpp> #include <ItemType.hpp> #include <RefNodeList.hpp> #include <RefFaceList.hpp> #include <SynchronizerManager.hpp> #include <tuple> #include <algorithm> #include <set> class ConnectivityDescriptor { private: std::vector<RefFaceList> m_ref_face_list; public: std::vector<std::vector<unsigned int>> cell_to_node_vector; #warning should be set in 2d std::vector<std::vector<unsigned int>> cell_to_face_vector; // std::vector<std::vector<unsigned int>> face_to_cell_vector; std::vector<std::vector<unsigned int>> face_to_node_vector; template <typename ItemOfItemT> auto& itemOfItemVector() { if constexpr (std::is_same_v<ItemOfItemT,NodeOfCell>) { return cell_to_node_vector; } else if constexpr (std::is_same_v<ItemOfItemT,FaceOfCell>) { return cell_to_face_vector; } else if constexpr (std::is_same_v<ItemOfItemT,NodeOfCell>) { return face_to_node_vector; } else { static_assert(is_false_v<ItemOfItemT>, "Unexpected item of item type"); } } std::vector<Array<bool>> cell_face_is_reversed_vector; std::vector<CellType> cell_type_vector; std::vector<int> cell_number_vector; std::vector<int> face_number_vector; std::vector<int> node_number_vector; std::vector<int> cell_owner_vector; std::vector<int> face_owner_vector; std::vector<int> node_owner_vector; const std::vector<RefFaceList>& refFaceList() const { return m_ref_face_list; } void addRefFaceList(const RefFaceList& ref_face_list) { m_ref_face_list.push_back(ref_face_list); } ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete; ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete; ConnectivityDescriptor() = default; ConnectivityDescriptor(const ConnectivityDescriptor&) = default; ConnectivityDescriptor(ConnectivityDescriptor&&) = delete; ~ConnectivityDescriptor() = default; }; template <size_t Dim> class Connectivity final : public IConnectivity { public: PASTIS_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; PASTIS_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; #warning is m_cell_global_index really necessary? should it be computed on demand instead? WeakCellValue<const int> m_cell_global_index; WeakCellValue<const int> m_cell_number; WeakFaceValue<const int> m_face_number; #warning check that m_edge_number is filled correctly 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; 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<RefFaceList> m_ref_face_list; std::vector<RefNodeList> m_ref_node_list; WeakCellValue<const double> m_inv_cell_nb_nodes; void _computeCellFaceAndFaceNodeConnectivities(); template <typename SubItemValuePerItemType> PASTIS_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; PASTIS_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: PASTIS_INLINE const auto& cellType() const { return m_cell_type; } PASTIS_INLINE const auto& cellNumber() const { return m_cell_number; } PASTIS_INLINE const auto& faceNumber() const { return m_face_number; } PASTIS_INLINE const auto& edgeNumber() const { return m_edge_number; } PASTIS_INLINE const auto& nodeNumber() const { return m_node_number; } template <ItemType item_type> PASTIS_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"); } } PASTIS_INLINE const auto& cellOwner() const { return m_cell_owner; } PASTIS_INLINE const auto& faceOwner() const { return m_face_owner; } PASTIS_INLINE const auto& edgeOwner() const { perr() << __FILE__ << ':' << __LINE__ << ": edge owner not built\n"; std::terminate(); return m_edge_owner; } PASTIS_INLINE const auto& nodeOwner() const { return m_node_owner; } template <ItemType item_type> PASTIS_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"); } } PASTIS_INLINE const auto& cellIsOwned() const { return m_cell_is_owned; } PASTIS_INLINE const auto& faceIsOwned() const { return m_face_is_owned; } PASTIS_INLINE const auto& edgeIsOwned() const { perr() << __FILE__ << ':' << __LINE__ << ": edge is owned not built\n"; std::terminate(); return m_edge_is_owned; } PASTIS_INLINE const auto& nodeIsOwned() const { return m_node_is_owned; } template <ItemType item_type> PASTIS_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"); } } PASTIS_INLINE const bool& isConnectivityMatrixBuilt(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)]; return connectivity_matrix.isBuilt(); } template <ItemType source_item_type, ItemType target_item_type> PASTIS_INLINE auto getItemToItemMatrix() const { return ItemToItemMatrix<source_item_type, target_item_type>(_getMatrix(source_item_type, target_item_type)); } PASTIS_INLINE auto cellToFaceMatrix() const { return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>(); } PASTIS_INLINE auto cellToEdgeMatrix() const { return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>(); } PASTIS_INLINE auto cellToNodeMatrix() const { return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>(); } PASTIS_INLINE auto faceToCellMatrix() const { return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>(); } PASTIS_INLINE auto faceToEdgeMatrix() const { return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>(); } PASTIS_INLINE auto faceToNodeMatrix() const { return this->template getItemToItemMatrix<ItemType::face, ItemType::node>(); } PASTIS_INLINE auto edgeToCellMatrix() const { return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>(); } PASTIS_INLINE auto edgeToFaceMatrix() const { return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>(); } PASTIS_INLINE auto edgeToNodeMatrix() const { return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>(); } PASTIS_INLINE auto nodeToCellMatrix() const { return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>(); } PASTIS_INLINE auto nodeToFaceMatrix() const { return this->template getItemToItemMatrix<ItemType::node, ItemType::face>(); } PASTIS_INLINE auto nodeToEdgeMatrix() const { return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>(); } PASTIS_INLINE const auto& cellFaceIsReversed() const { static_assert(Dimension>1, "reversed faces makes no sense in dimension 1"); return m_cell_face_is_reversed; } PASTIS_INLINE const auto& cellLocalNumbersInTheirNodes() const { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes); } PASTIS_INLINE const auto& cellLocalNumbersInTheirEdges() const { if constexpr (Dimension>2) { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_edges); } else { return cellLocalNumbersInTheirFaces(); } } PASTIS_INLINE const auto& cellLocalNumbersInTheirFaces() const { if constexpr (Dimension>1) { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_faces); } else { return cellLocalNumbersInTheirNodes(); } } PASTIS_INLINE const auto& faceLocalNumbersInTheirCells() const { if constexpr(Dimension>1) { return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_cells); } else { return nodeLocalNumbersInTheirCells(); } } PASTIS_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); } PASTIS_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); } PASTIS_INLINE const auto& edgeLocalNumbersInTheirCells() const { if constexpr (Dimension>2) { return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_cells); } else { return faceLocalNumbersInTheirCells(); } } PASTIS_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); } PASTIS_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); } PASTIS_INLINE const auto& nodeLocalNumbersInTheirCells() const { return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells); } PASTIS_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); } PASTIS_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); } size_t numberOfRefFaceList() const { return m_ref_face_list.size(); } const RefFaceList& refFaceList(const size_t& i) const { return m_ref_face_list[i]; } void addRefNodeList(const RefNodeList& ref_node_list) { m_ref_node_list.push_back(ref_node_list); } size_t numberOfRefNodeList() const { return m_ref_node_list.size(); } const RefNodeList& refNodeList(const size_t& i) const { return m_ref_node_list[i]; } PASTIS_INLINE CSRGraph cellToCellGraph() const { std::vector<std::set<int>> cell_cells(this->numberOfCells()); if constexpr (true) { 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); } } } else { const auto& node_to_cell_matrix = this->nodeToCellMatrix(); for (NodeId l=0; l<this->numberOfNodes(); ++l) { const auto& node_to_cell = node_to_cell_matrix[l]; for (size_t i_cell=0; i_cell<node_to_cell.size(); ++i_cell) { const CellId cell_0 = node_to_cell[i_cell]; for (size_t j_cell=0; j_cell<i_cell; ++j_cell) { const CellId cell_1 = node_to_cell[j_cell]; 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 CSRGraph(entries, neighbors); } PASTIS_INLINE size_t numberOfNodes() const final { const auto& node_to_cell_matrix = this->_getMatrix(ItemType::node,ItemType::cell); return node_to_cell_matrix.numRows(); } PASTIS_INLINE size_t numberOfEdges() const final { const auto& edge_to_node_matrix = this->_getMatrix(ItemType::edge,ItemType::node); return edge_to_node_matrix.numRows(); } PASTIS_INLINE size_t numberOfFaces() const final { const auto& face_to_node_matrix = this->_getMatrix(ItemType::face,ItemType::cell); return face_to_node_matrix.numRows(); } PASTIS_INLINE size_t numberOfCells() const final { const auto& cell_to_node_matrix = this->_getMatrix(ItemType::cell,ItemType::node); return cell_to_node_matrix.numRows(); } CellValue<const double> invCellNbNodes() const { #warning add calculation on demand when variables will be defined return m_inv_cell_nb_nodes; } Connectivity(const Connectivity&) = delete; private: Connectivity(); void _buildFrom(const ConnectivityDescriptor& descriptor); public: ~Connectivity() { auto& manager = SynchronizerManager::instance(); manager.deleteConnectivitySynchronizer(this); } }; template <size_t Dim> PASTIS_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