diff --git a/src/mesh/MeshBoundary.hpp b/src/mesh/MeshBoundary.hpp index 8e3961e52f5fc4b7cfc1107254e1e5a265f19122..46a5f55d0d7e96a19aeaec58770bf37c672b1a0b 100644 --- a/src/mesh/MeshBoundary.hpp +++ b/src/mesh/MeshBoundary.hpp @@ -10,6 +10,8 @@ template <size_t dimension> class MeshBoundary { + protected: + Kokkos::View<const unsigned int*> m_node_list; public: MeshBoundary& operator=(const MeshBoundary&) = default; MeshBoundary& operator=(MeshBoundary&&) = default; @@ -27,6 +29,28 @@ class MeshBoundary std::exit(1); } }); + + const Kokkos::View<const unsigned short*> face_nb_nodes = mesh.connectivity().faceNbNodes(); + const Kokkos::View<const unsigned int**> face_nodes = mesh.connectivity().faceNodes(); + + std::vector<unsigned int> node_ids; + // not enough but should reduce significantly the number of resizing + node_ids.reserve(dimension*face_list.extent(0)); + for (size_t l=0; l<face_list.extent(0); ++l) { + const size_t face_number = face_list[l]; + for (size_t r=0; r<face_nb_nodes[face_number]; ++r) { + node_ids.push_back(face_nodes(face_number,r)); + } + } + std::sort(node_ids.begin(), node_ids.end()); + auto last = std::unique(node_ids.begin(), node_ids.end()); + node_ids.resize(std::distance(node_ids.begin(), last)); + + Kokkos::View<unsigned int*> node_list("node_list", node_ids.size()); + Kokkos::parallel_for(node_ids.size(), KOKKOS_LAMBDA(const int& r){ + node_list[r] = node_ids[r]; + }); + m_node_list = node_list; } MeshBoundary() = default; @@ -45,8 +69,7 @@ class MeshFlatBoundary const Rd m_outgoing_normal; template <typename MeshType> - inline Rd _getOutgoingNormal(const MeshType& mesh, - const Kokkos::View<const unsigned int*>& face_list); + inline Rd _getOutgoingNormal(const MeshType& mesh); public: MeshFlatBoundary& operator=(const MeshFlatBoundary&) = default; MeshFlatBoundary& operator=(MeshFlatBoundary&&) = default; @@ -55,7 +78,7 @@ class MeshFlatBoundary MeshFlatBoundary(const MeshType& mesh, const Kokkos::View<const unsigned int*>& face_list) : MeshBoundary<dimension>(mesh, face_list), - m_outgoing_normal(_getOutgoingNormal(mesh, face_list)) + m_outgoing_normal(_getOutgoingNormal(mesh)) { ; } @@ -71,14 +94,11 @@ template <> template <typename MeshType> inline TinyVector<2,double> MeshFlatBoundary<2>:: -_getOutgoingNormal(const MeshType& mesh, - const Kokkos::View<const unsigned int*>& face_list) +_getOutgoingNormal(const MeshType& mesh) { static_assert(MeshType::dimension == 2); typedef TinyVector<2,double> R2; - const Kokkos::View<const unsigned short*> face_nb_nodes = mesh.connectivity().faceNbNodes(); - const Kokkos::View<const unsigned int**> face_nodes = mesh.connectivity().faceNodes(); const Kokkos::View<const R2*> xr = mesh.xr(); R2 xmin(std::numeric_limits<double>::max(), @@ -87,21 +107,13 @@ _getOutgoingNormal(const MeshType& mesh, R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()); - if (face_list.extent(0)==0) { - std::cerr << "face_list is empty! Cannot compute normal!\n"; - std::exit(1); - } - - for (size_t l=0; l<face_list.extent(0); ++l) { - const unsigned int face_number = face_list[l]; - for (int r=0; r<face_nb_nodes(face_number); ++r) { - const R2& x = xr(face_nodes(face_number,r)); - if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) { - xmin = x; - } - if ((x[0]>xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) { - xmax = x; - } + for (size_t r=0; r<m_node_list.extent(0); ++r) { + const R2& x = xr[m_node_list[r]]; + if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) { + xmin = x; + } + if ((x[0]>xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) { + xmax = x; } }