#ifndef CONNECTIVITY_HPP #define CONNECTIVITY_HPP #include <PastisAssert.hpp> #include <TinyVector.hpp> #include <Kokkos_Core.hpp> #include <ItemValue.hpp> #include <IConnectivity.hpp> #include <ConnectivityMatrix.hpp> #include <SubItemValuePerItem.hpp> #include <ConnectivityComputer.hpp> #include <vector> #include <unordered_map> #include <algorithm> #include <CellType.hpp> #include <RefId.hpp> #include <ItemType.hpp> #include <RefNodeList.hpp> #include <RefFaceList.hpp> #include <tuple> #include <algorithm> template <size_t Dimension> class Connectivity; template <size_t Dimension> class ConnectivityFace; template<> class ConnectivityFace<1> { public: friend struct Hash; struct Hash { size_t operator()(const ConnectivityFace& f) const; }; }; template<> class ConnectivityFace<2> { public: friend struct Hash; struct Hash { size_t operator()(const ConnectivityFace& f) const { size_t hash = 0; hash ^= std::hash<unsigned int>()(f.m_node0_id); hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1; return hash; } }; unsigned int m_node0_id; unsigned int m_node1_id; friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f) { os << f.m_node0_id << ' ' << f.m_node1_id << ' '; return os; } KOKKOS_INLINE_FUNCTION bool operator==(const ConnectivityFace& f) const { return ((m_node0_id == f.m_node0_id) and (m_node1_id == f.m_node1_id)); } KOKKOS_INLINE_FUNCTION bool operator<(const ConnectivityFace& f) const { return ((m_node0_id<f.m_node0_id) or ((m_node0_id == f.m_node0_id) and (m_node1_id<f.m_node1_id))); } KOKKOS_INLINE_FUNCTION ConnectivityFace& operator=(const ConnectivityFace&) = default; KOKKOS_INLINE_FUNCTION ConnectivityFace& operator=(ConnectivityFace&&) = default; KOKKOS_INLINE_FUNCTION ConnectivityFace(const std::vector<unsigned int>& given_node_id_list) { Assert(given_node_id_list.size()==2); #warning rework this dirty constructor const auto& [min, max] = std::minmax(given_node_id_list[0], given_node_id_list[1]); m_node0_id = min; m_node1_id = max; } KOKKOS_INLINE_FUNCTION ConnectivityFace(const ConnectivityFace&) = default; KOKKOS_INLINE_FUNCTION ConnectivityFace(ConnectivityFace&&) = default; KOKKOS_INLINE_FUNCTION ~ConnectivityFace() = default; }; template <> class ConnectivityFace<3> { private: friend class Connectivity<3>; friend struct Hash; struct Hash { size_t operator()(const ConnectivityFace& f) const { size_t hash = 0; for (size_t i=0; i<f.m_node_id_list.size(); ++i) { hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i; } return hash; } }; bool m_reversed; std::vector<unsigned int> m_node_id_list; friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f) { for (auto id : f.m_node_id_list) { std::cout << id << ' '; } return os; } KOKKOS_INLINE_FUNCTION const bool& reversed() const { return m_reversed; } KOKKOS_INLINE_FUNCTION const std::vector<unsigned int>& nodeIdList() const { return m_node_id_list; } KOKKOS_INLINE_FUNCTION std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list) { const auto min_id = std::min_element(node_list.begin(), node_list.end()); const int shift = std::distance(node_list.begin(), min_id); std::vector<unsigned int> rotated_node_list(node_list.size()); if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) { for (size_t i=0; i<node_list.size(); ++i) { rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()]; m_reversed = true; } } else { for (size_t i=0; i<node_list.size(); ++i) { rotated_node_list[i] = node_list[(shift+i)%node_list.size()]; } } return rotated_node_list; } KOKKOS_INLINE_FUNCTION ConnectivityFace(const std::vector<unsigned int>& given_node_id_list) : m_reversed(false), m_node_id_list(_sort(given_node_id_list)) { ; } public: bool operator==(const ConnectivityFace& f) const { if (m_node_id_list.size() == f.nodeIdList().size()) { for (size_t j=0; j<m_node_id_list.size(); ++j) { if (m_node_id_list[j] != f.nodeIdList()[j]) { return false; } } return true; } return false; } KOKKOS_INLINE_FUNCTION bool operator<(const ConnectivityFace& f) const { const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size()); for (size_t i=0; i<min_nb_nodes; ++i) { if (m_node_id_list[i] < f.m_node_id_list[i]) return true; if (m_node_id_list[i] != f.m_node_id_list[i]) return false; } return m_node_id_list.size() < f.m_node_id_list.size(); } KOKKOS_INLINE_FUNCTION ConnectivityFace& operator=(const ConnectivityFace&) = default; KOKKOS_INLINE_FUNCTION ConnectivityFace& operator=(ConnectivityFace&&) = default; KOKKOS_INLINE_FUNCTION ConnectivityFace(const ConnectivityFace&) = default; KOKKOS_INLINE_FUNCTION ConnectivityFace(ConnectivityFace&&) = default; KOKKOS_INLINE_FUNCTION ConnectivityFace() = delete; KOKKOS_INLINE_FUNCTION ~ConnectivityFace() = default; }; template <size_t Dimension> class Connectivity final : public IConnectivity { private: constexpr static auto& itemTId = ItemTypeId<Dimension>::itemTId; public: static constexpr size_t dimension = Dimension; private: ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1]; CellValue<const CellType> m_cell_type; FaceValuePerCell<const bool> m_cell_face_is_reversed; NodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes; EdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges; FaceValuePerCell<const unsigned short> m_cell_local_numbers_in_their_faces; CellValuePerFace<const unsigned short> m_face_local_numbers_in_their_cells; EdgeValuePerFace<const unsigned short> m_face_local_numbers_in_their_edges; NodeValuePerFace<const unsigned short> m_face_local_numbers_in_their_nodes; CellValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_cells; FaceValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_faces; NodeValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_nodes; CellValuePerNode<const unsigned short> m_node_local_numbers_in_their_cells; EdgeValuePerNode<const unsigned short> m_node_local_numbers_in_their_edges; FaceValuePerNode<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; CellValue<const double> m_inv_cell_nb_nodes; using Face = ConnectivityFace<Dimension>; std::unordered_map<Face, unsigned int, typename Face::Hash> m_face_number_map; void _computeCellFaceAndFaceNodeConnectivities(); template <typename SubItemValuePerItemType> KOKKOS_INLINE_FUNCTION const SubItemValuePerItemType& _lazzyBuildItemNumberInTheirChild(const SubItemValuePerItemType& sub_item_value_per_item) const { if (not sub_item_value_per_item.isBuilt()) { const_cast<SubItemValuePerItemType&>(sub_item_value_per_item) = m_connectivity_computer . computeLocalItemNumberInChildItem<SubItemValuePerItemType::item_t, SubItemValuePerItemType::sub_item_t>(*this); } return sub_item_value_per_item; } public: KOKKOS_INLINE_FUNCTION 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(); } KOKKOS_INLINE_FUNCTION 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; } KOKKOS_INLINE_FUNCTION const auto& cellFaceIsReversed() const { static_assert(dimension>1, "reversed faces makes no sense in dimension 1"); return m_cell_face_is_reversed; } KOKKOS_INLINE_FUNCTION const auto& cellLocalNumbersInTheirNodes() const { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes); } KOKKOS_INLINE_FUNCTION const auto& cellLocalNumbersInTheirEdges() const { if constexpr (dimension>2) { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_edges); } else { return cellLocalNumbersInTheirFaces(); } } KOKKOS_INLINE_FUNCTION const auto& cellLocalNumbersInTheirFaces() const { if constexpr (dimension>1) { return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_faces); } else { return cellLocalNumbersInTheirNodes(); } } KOKKOS_INLINE_FUNCTION const auto& faceLocalNumbersInTheirCells() const { if constexpr(dimension>1) { return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_cells); } else { return nodeLocalNumbersInTheirCells(); } } KOKKOS_INLINE_FUNCTION 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); } KOKKOS_INLINE_FUNCTION const auto& faceLocalNumbersInTheirNodes() const { static_assert(dimension>1,"this function has no meaning in 1d"); return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_nodes); } KOKKOS_INLINE_FUNCTION const auto& edgeLocalNumbersInTheirCells() const { if constexpr (dimension>2) { return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_cells); } else { return faceLocalNumbersInTheirCells(); } } KOKKOS_INLINE_FUNCTION 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); } KOKKOS_INLINE_FUNCTION 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); } KOKKOS_INLINE_FUNCTION const auto& nodeLocalNumbersInTheirCells() const { return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells); } KOKKOS_INLINE_FUNCTION 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); } KOKKOS_INLINE_FUNCTION const auto& nodeLocalNumbersInTheirFaces() const { static_assert(dimension>1,"this function has no meaning in 1d"); return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces); } void addRefFaceList(const RefFaceList& ref_face_list) { m_ref_face_list.push_back(ref_face_list); } 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]; } KOKKOS_INLINE_FUNCTION size_t numberOfNodes() const final { const auto& node_to_cell_matrix = this->getMatrix(ItemType::node,ItemType::cell); return node_to_cell_matrix.numRows(); } KOKKOS_INLINE_FUNCTION size_t numberOfEdges() const final { const auto& edge_to_node_matrix = this->getMatrix(ItemType::edge,ItemType::node); return edge_to_node_matrix.numRows(); } KOKKOS_INLINE_FUNCTION size_t numberOfFaces() const final { const auto& face_to_node_matrix = this->getMatrix(ItemType::face,ItemType::cell); return face_to_node_matrix.numRows(); } KOKKOS_INLINE_FUNCTION size_t numberOfCells() const final { const auto& cell_to_node_matrix = this->getMatrix(ItemType::cell,ItemType::node); return cell_to_node_matrix.numRows(); } const CellValue<const double>& invCellNbNodes() const { #warning add calculation on demand when variables will be defined return m_inv_cell_nb_nodes; } unsigned int getFaceNumber(const std::vector<unsigned int>& face_nodes) const { const Face face(face_nodes); auto i_face = m_face_number_map.find(face); if (i_face == m_face_number_map.end()) { std::cerr << "Face " << face << " not found!\n"; throw std::exception(); std::exit(0); } return i_face->second; } Connectivity(const Connectivity&) = delete; Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector, const std::vector<CellType>& cell_type_vector); ~Connectivity() { ; } }; using Connectivity3D = Connectivity<3>; using Connectivity2D = Connectivity<2>; using Connectivity1D = Connectivity<1>; #endif // CONNECTIVITY_HPP