#ifndef CONNECTIVITY_HPP #define CONNECTIVITY_HPP #include <PastisAssert.hpp> #include <TinyVector.hpp> #include <Kokkos_Core.hpp> #include <ConnectivityMatrix.hpp> #include <ConnectivityComputer.hpp> #include <SubItemValuePerItem.hpp> #include <vector> #include <unordered_map> #include <algorithm> #include <RefId.hpp> #include <TypeOfItem.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 { private: constexpr static auto& itemId = ItemId<Dimension>::itemId; public: static constexpr size_t dimension = Dimension; template <TypeOfItem item_type_0, TypeOfItem item_type_1> const ConnectivityMatrix& itemToItemMatrix() const { return m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)]; }; const auto& cellFaceIsReversed() const { return m_cell_face_is_reversed; } const auto& cellToNodeLocalCell() const { return m_cell_to_node_local_cell; } const auto& nodeToCellLocalNode() const { return m_node_to_cell_local_node; } const auto& faceToCellLocalFace() const { if constexpr(dimension>1) { return m_face_to_cell_local_face; } else { return m_node_to_cell_local_node; } } const auto& nodeToFaceLocalNode() const { static_assert("this function has no meaning in 1d"); return m_node_to_face_local_node; } private: ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1]; FaceValuePerCell<const bool> m_cell_face_is_reversed; NodeValuePerCell<const unsigned short> m_cell_to_node_local_cell; CellValuePerFace<const unsigned short> m_face_to_cell_local_face; CellValuePerNode<const unsigned short> m_node_to_cell_local_node; // not plugged ... #warning remaining def NodeValuePerFace<unsigned short> m_node_to_face_local_node; FaceValuePerNode<unsigned short> m_face_to_node_local_face; // ... not plugged ConnectivityComputer m_connectivity_computer; std::vector<RefFaceList> m_ref_face_list; std::vector<RefNodeList> m_ref_node_list; Kokkos::View<double*> m_inv_cell_nb_nodes; Kokkos::View<const unsigned short*> m_node_nb_faces; Kokkos::View<const unsigned int**> m_node_faces; using Face = ConnectivityFace<Dimension>; std::unordered_map<Face, unsigned int, typename Face::Hash> m_face_number_map; void _computeFaceCellConnectivities(); public: 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 { const auto& node_to_cell_matrix = m_item_to_item_matrix[itemId(TypeOfItem::node)][itemId(TypeOfItem::cell)]; return node_to_cell_matrix.numRows(); } KOKKOS_INLINE_FUNCTION size_t numberOfFaces() const { const auto& face_to_node_matrix = m_item_to_item_matrix[itemId(TypeOfItem::face)][itemId(TypeOfItem::node)]; return face_to_node_matrix.numRows(); } KOKKOS_INLINE_FUNCTION size_t numberOfCells() const { const auto& cell_to_node_matrix = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)]; return cell_to_node_matrix.numRows(); } const Kokkos::View<const double*> invCellNbNodes() const { 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) { auto& cell_to_node_matrix = m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)]; cell_to_node_matrix = cell_by_node_vector; Assert(this->numberOfCells()>0); { Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells()); Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){ const auto& cell_nodes = cell_to_node_matrix.rowConst(j); inv_cell_nb_nodes[j] = 1./cell_nodes.length; }); m_inv_cell_nb_nodes = inv_cell_nb_nodes; } auto& node_to_cell_matrix = m_item_to_item_matrix[itemId(TypeOfItem::node)][itemId(TypeOfItem::cell)]; m_connectivity_computer.computeInverseConnectivityMatrix(cell_to_node_matrix, node_to_cell_matrix); CellValuePerNode<unsigned short> node_to_cell_local_node(*this); m_connectivity_computer.computeLocalChildItemNumberInItem(cell_to_node_matrix, node_to_cell_matrix, node_to_cell_local_node); m_node_to_cell_local_node = node_to_cell_local_node; NodeValuePerCell<unsigned short> cell_to_node_local_cell(*this); m_connectivity_computer.computeLocalChildItemNumberInItem(node_to_cell_matrix, cell_to_node_matrix, cell_to_node_local_cell); m_cell_to_node_local_cell = cell_to_node_local_cell; if constexpr (Dimension>1) { this->_computeFaceCellConnectivities(); } } ~Connectivity() { ; } }; using Connectivity3D = Connectivity<3>; using Connectivity2D = Connectivity<2>; using Connectivity1D = Connectivity<1>; #endif // CONNECTIVITY_HPP