#ifndef CONNECTIVITY_3D_HPP #define CONNECTIVITY_3D_HPP #include <Kokkos_Core.hpp> #include <PastisAssert.hpp> #include <TinyVector.hpp> #include <ConnectivityUtils.hpp> #include <vector> #include <map> #include <algorithm> #include <RefId.hpp> #include <RefNodeList.hpp> #include <RefFaceList.hpp> #include <tuple> class Connectivity3D { public: static constexpr size_t dimension = 3; private: std::vector<RefFaceList> m_ref_face_list; std::vector<RefNodeList> m_ref_node_list; const size_t m_number_of_cells; size_t m_number_of_faces; size_t m_number_of_nodes; const Kokkos::View<const unsigned short*> m_cell_nb_nodes; const Kokkos::View<const unsigned int**> m_cell_nodes; Kokkos::View<double*> m_inv_cell_nb_nodes; Kokkos::View<const unsigned short*> m_cell_nb_faces; Kokkos::View<const unsigned int**> m_cell_faces; Kokkos::View<const bool**> m_cell_faces_is_reversed; Kokkos::View<unsigned short*> m_node_nb_cells; Kokkos::View<unsigned int**> m_node_cells; Kokkos::View<unsigned short**> m_node_cell_local_node; Kokkos::View<unsigned short*> m_face_nb_cells; Kokkos::View<unsigned int**> m_face_cells; Kokkos::View<unsigned short**> m_face_cell_local_face; Kokkos::View<unsigned short*> m_face_nb_nodes; Kokkos::View<unsigned int**> m_face_nodes; Kokkos::View<unsigned short**> m_face_node_local_face; size_t m_max_nb_node_per_cell; size_t m_max_nb_face_per_cell; size_t m_max_nb_node_per_face; class Face { private: bool m_reversed; std::vector<unsigned int> m_node_id_list; public: friend std::ostream& operator<<(std::ostream& os, const Face& 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; } bool operator==(const Face& 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 Face& 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 Face& operator=(const Face&) = default; KOKKOS_INLINE_FUNCTION Face& operator=(Face&&) = default; KOKKOS_INLINE_FUNCTION Face(const Face&) = default; KOKKOS_INLINE_FUNCTION Face(Face&&) = default; KOKKOS_INLINE_FUNCTION Face(const std::vector<unsigned int>& given_node_id_list) : m_reversed(false), m_node_id_list(_sort(given_node_id_list)) { ; } KOKKOS_INLINE_FUNCTION Face() = delete; KOKKOS_INLINE_FUNCTION ~Face() = default; }; void _computeFaceCellConnectivities() { Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", m_number_of_cells); typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo; std::map<Face, std::vector<CellFaceInfo>> face_cells_map; for (unsigned int j=0; j<m_number_of_cells; ++j) { const unsigned short cell_nb_nodes = m_cell_nb_nodes[j]; switch (cell_nb_nodes) { case 4: { // tetrahedron cell_nb_faces[j] = 4; // face 0 Face f0({m_cell_nodes(j,1), m_cell_nodes(j,2), m_cell_nodes(j,3)}); face_cells_map[f0].push_back(std::make_tuple(j, 0, f0.reversed())); // face 1 Face f1({m_cell_nodes(j,0), m_cell_nodes(j,3), m_cell_nodes(j,2)}); face_cells_map[f1].push_back(std::make_tuple(j, 1, f1.reversed())); // face 2 Face f2({m_cell_nodes(j,0), m_cell_nodes(j,1), m_cell_nodes(j,3)}); face_cells_map[f2].push_back(std::make_tuple(j, 2, f2.reversed())); // face 3 Face f3({m_cell_nodes(j,0), m_cell_nodes(j,2), m_cell_nodes(j,1)}); face_cells_map[f3].push_back(std::make_tuple(j, 3, f3.reversed())); break; } case 8: { // hexahedron // face 0 Face f0({m_cell_nodes(j,3), m_cell_nodes(j,2), m_cell_nodes(j,1), m_cell_nodes(j,0)}); face_cells_map[f0].push_back(std::make_tuple(j, 0, f0.reversed())); // face 1 Face f1({m_cell_nodes(j,4), m_cell_nodes(j,5), m_cell_nodes(j,6), m_cell_nodes(j,7)}); face_cells_map[f1].push_back(std::make_tuple(j, 1, f1.reversed())); // face 2 Face f2({m_cell_nodes(j,0), m_cell_nodes(j,4), m_cell_nodes(j,7), m_cell_nodes(j,3)}); face_cells_map[f2].push_back(std::make_tuple(j, 2, f2.reversed())); // face 3 Face f3({m_cell_nodes(j,1), m_cell_nodes(j,2), m_cell_nodes(j,6), m_cell_nodes(j,5)}); face_cells_map[f3].push_back(std::make_tuple(j, 3, f3.reversed())); // face 4 Face f4({m_cell_nodes(j,0), m_cell_nodes(j,1), m_cell_nodes(j,5), m_cell_nodes(j,4)}); face_cells_map[f4].push_back(std::make_tuple(j, 4, f4.reversed())); // face 5 Face f5({m_cell_nodes(j,3), m_cell_nodes(j,7), m_cell_nodes(j,6), m_cell_nodes(j,2)}); face_cells_map[f5].push_back(std::make_tuple(j, 5, f5.reversed())); cell_nb_faces[j] = 6; break; } default: { std::cerr << "unexpected cell type!\n"; std::exit(0); } } } std::cerr << __FILE__ << ':' << __LINE__ << ':' << rang::fg::red << " only works for hexa" << rang::fg::reset << '\n'; m_cell_nb_faces = cell_nb_faces; m_number_of_faces = face_cells_map.size(); std::cerr << __FILE__ << ':' << __LINE__ << ':' << rang::fg::red << " m_max_nb_node_per_face forced to 4" << rang::fg::reset << '\n'; Kokkos::View<unsigned short*> face_nb_nodes("face_nb_nodes", m_number_of_faces); m_max_nb_node_per_face = 4; Kokkos::View<unsigned int**> face_nodes("face_nodes", m_number_of_faces, m_max_nb_node_per_face); { int l=0; for (const auto& face_cells_vector : face_cells_map) { const Face& face = face_cells_vector.first; for(size_t r=0; r<face.nodeIdList().size(); ++r) { face_nodes(l,r) = face.nodeIdList()[r]; } face_nb_nodes[l] = face.nodeIdList().size(); ++l; } } m_face_nodes = face_nodes; m_face_nb_nodes = face_nb_nodes; Kokkos::View<unsigned short*> face_nb_cells("face_nb_cells", m_number_of_faces); { int l=0; for (const auto& face_cells_vector : face_cells_map) { const auto& cells_vector = face_cells_vector.second; face_nb_cells[l] = cells_vector.size(); ++l; } } m_face_nb_cells = face_nb_cells; #warning check that the number of cell per faces is <=2 Kokkos::View<unsigned int**> face_cells("face_cells", face_cells_map.size(), 2); { int l=0; for (const auto& face_cells_vector : face_cells_map) { const auto& cells_vector = face_cells_vector.second; for (unsigned short lj=0; lj<cells_vector.size(); ++lj) { const unsigned int cell_number = std::get<0>(cells_vector[lj]); face_cells(l,lj) = cell_number; } ++l; } } m_face_cells = face_cells; std::cerr << __FILE__ << ':' << __LINE__ << ':' << rang::fg::red << " m_max_nb_face_per_cell forced to 6" << rang::fg::reset << '\n'; m_max_nb_face_per_cell = 6; Kokkos::View<unsigned int**> cell_faces("cell_faces", m_number_of_cells, m_max_nb_face_per_cell); { int l=0; for (const auto& face_cells_vector : face_cells_map) { const auto& cells_vector = face_cells_vector.second; for (unsigned short lj=0; lj<cells_vector.size(); ++lj) { unsigned int cell_number; unsigned short cell_local_face; std::tie(cell_number, cell_local_face, std::ignore) = cells_vector[lj]; cell_faces(cell_number,cell_local_face) = l; } ++l; } } m_cell_faces = cell_faces; Kokkos::View<bool**> cell_faces_is_reversed("cell_faces_is_reversed", m_number_of_cells, m_max_nb_face_per_cell); { int l=0; for (const auto& face_cells_vector : face_cells_map) { const auto& cells_vector = face_cells_vector.second; for (unsigned short lj=0; lj<cells_vector.size(); ++lj) { unsigned int cell_number; unsigned short cell_local_face; bool reversed; std::tie(cell_number, cell_local_face, reversed) = cells_vector[lj]; cell_faces_is_reversed(cell_number,cell_local_face) = reversed; } ++l; } } m_cell_faces_is_reversed = cell_faces_is_reversed; Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face", m_number_of_faces, m_max_nb_node_per_cell); { int l=0; for (const auto& face_cells_vector : face_cells_map) { const auto& cells_vector = face_cells_vector.second; for (unsigned short lj=0; lj<cells_vector.size(); ++lj) { unsigned short cell_local_face = std::get<1>(cells_vector[lj]); face_cell_local_face(l,lj) = cell_local_face; } ++l; } } m_face_cell_local_face = face_cell_local_face; } 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]; } const size_t& numberOfNodes() const { return m_number_of_nodes; } const size_t& numberOfFaces() const { return m_number_of_faces; } const size_t& numberOfCells() const { return m_number_of_cells; } const size_t& maxNbFacePerCell() const { return m_max_nb_face_per_cell; } const size_t& maxNbNodePerCell() const { return m_max_nb_node_per_cell; } const size_t& maxNbNodePerFace() const { return m_max_nb_node_per_face; } const Kokkos::View<const unsigned int**> cellNodes() const { return m_cell_nodes; } const Kokkos::View<const unsigned int**> cellFaces() const { return m_cell_faces; } const Kokkos::View<const bool**> cellFacesIsReversed() const { return m_cell_faces_is_reversed; } const Kokkos::View<const unsigned short*> nodeNbCells() const { return m_node_nb_cells; } const Kokkos::View<const unsigned short*> cellNbNodes() const { return m_cell_nb_nodes; } const Kokkos::View<const double*> invCellNbNodes() const { return m_inv_cell_nb_nodes; } const Kokkos::View<const unsigned short*> cellNbFaces() const { return m_cell_nb_faces; } const Kokkos::View<const unsigned short*> faceNbCells() const { return m_face_nb_cells; } const Kokkos::View<const unsigned short*> faceNbNodes() const { return m_face_nb_nodes; } const Kokkos::View<const unsigned int**> nodeCells() const { return m_node_cells; } const Kokkos::View<const unsigned int**> faceCells() const { return m_face_cells; } const Kokkos::View<const unsigned int**> faceNodes() const { return m_face_nodes; } const Kokkos::View<const unsigned short**> nodeCellLocalNode() const { return m_node_cell_local_node; } const Kokkos::View<const unsigned short**> faceCellLocalFace() const { return m_face_cell_local_face; } unsigned int getFaceNumber(const std::vector<unsigned int>& face_nodes) const { #warning rewrite // worst code ever const Face face(face_nodes); for (unsigned int l=0; l<m_number_of_faces; ++l) { std::vector<unsigned int> nodes(m_face_nb_nodes[l]); for (size_t r=0; r<nodes.size(); ++r) { nodes[r] = m_face_nodes(l,r); } Face f(nodes); if (f == face) { return l; } } std::cerr << "Face not found!\n"; std::exit(0); return 0; } Connectivity3D(const Connectivity3D&) = delete; Connectivity3D(const Kokkos::View<const unsigned short*> cell_nb_nodes, const Kokkos::View<const unsigned int**> cell_nodes) : m_number_of_cells (cell_nodes.extent(0)), m_cell_nb_nodes(cell_nb_nodes), m_cell_nodes (cell_nodes) { { Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells); Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int&j){ inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j]; }); m_inv_cell_nb_nodes = inv_cell_nb_nodes; } Assert(m_number_of_cells>0); ConnectivityUtils utils; utils.computeNodeCellConnectivity(m_cell_nodes, m_cell_nb_nodes, m_number_of_cells, m_max_nb_node_per_cell, m_number_of_nodes, m_node_nb_cells, m_node_cells, m_node_cell_local_node); this->_computeFaceCellConnectivities(); } ~Connectivity3D() { ; } }; #endif // CONNECTIVITY_3D_HPP