#ifndef MESH_NODE_BOUNDARY_HPP #define MESH_NODE_BOUNDARY_HPP #include <Kokkos_Core.hpp> #include <TinyVector.hpp> #include <iostream> template <size_t dimension> class MeshNodeBoundary { protected: Kokkos::View<const unsigned int*> m_node_list; public: MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default; MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default; template <typename MeshType> MeshNodeBoundary(const MeshType& mesh, const Kokkos::View<const unsigned int*>& face_list) { static_assert(dimension == MeshType::dimension); const Kokkos::View<const unsigned short*> face_nb_cells = mesh.connectivity().faceNbCells(); Kokkos::parallel_for(face_list.extent(0), KOKKOS_LAMBDA(const int& l){ if (face_nb_cells[face_list[l]]>1) { std::cerr << "internal faces cannot be used to define mesh boundaries\n"; 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; } MeshNodeBoundary() = default; MeshNodeBoundary(const MeshNodeBoundary&) = default; MeshNodeBoundary(MeshNodeBoundary&&) = default; virtual ~MeshNodeBoundary() = default; }; template <size_t dimension> class MeshFlatNodeBoundary : public MeshNodeBoundary<dimension> { typedef TinyVector<dimension, double> Rd; private: const Rd m_outgoing_normal; template <typename MeshType> inline Rd _getNormal(const MeshType& mesh); template <typename MeshType> inline void _checkBoundaryIsFlat(const TinyVector<2,double>& normal, const TinyVector<2,double>& xmin, const TinyVector<2,double>& xmax, const MeshType& mesh) const; template <typename MeshType> inline Rd _getOutgoingNormal(const MeshType& mesh); public: MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default; MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default; template <typename MeshType> MeshFlatNodeBoundary(const MeshType& mesh, const Kokkos::View<const unsigned int*>& face_list) : MeshNodeBoundary<dimension>(mesh, face_list), m_outgoing_normal(_getOutgoingNormal(mesh)) { ; } MeshFlatNodeBoundary() = default; MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default; MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default; virtual ~MeshFlatNodeBoundary() = default; }; template<> template <typename MeshType> void MeshFlatNodeBoundary<2>:: _checkBoundaryIsFlat(const TinyVector<2,double>& normal, const TinyVector<2,double>& xmin, const TinyVector<2,double>& xmax, const MeshType& mesh) const { static_assert(MeshType::dimension == 2); typedef TinyVector<2,double> R2; const R2 origin = 0.5*(xmin+xmax); const double length = l2Norm(xmax-xmin); const Kokkos::View<const R2*> xr = mesh.xr(); Kokkos::parallel_for(m_node_list.extent(0), KOKKOS_LAMBDA(const size_t& r) { const R2& x = xr[m_node_list[r]]; if ((x-origin,normal)>1E-13*length) { std::cerr << "this FlatBoundary is not flat!\n"; std::exit(1); } }); } template <> template <typename MeshType> inline TinyVector<2,double> MeshFlatNodeBoundary<2>:: _getNormal(const MeshType& mesh) { static_assert(MeshType::dimension == 2); typedef TinyVector<2,double> R2; const Kokkos::View<const R2*> xr = mesh.xr(); R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max()); R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()); 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; } } if (xmin == xmax) { std::cerr << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal"; std::exit(1); } R2 dx = xmax-xmin; dx *= 1./l2Norm(dx); R2 normal(-dx[1], dx[0]); this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh); return normal; } template <> template <typename MeshType> inline TinyVector<2,double> MeshFlatNodeBoundary<2>:: _getOutgoingNormal(const MeshType& mesh) { static_assert(MeshType::dimension == 2); typedef TinyVector<2,double> R2; const R2 normal = this->_getNormal(mesh); const Kokkos::View<const R2*>& xr = mesh.xr(); const Kokkos::View<const unsigned int**>& node_cells = mesh.connectivity().nodeCells(); const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes(); const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes(); const size_t r0 = m_node_list[0]; const size_t j0 = node_cells(r0,0); double max_height = 0; for (int r=0; r<cell_nb_nodes(j0); ++r) { const double height = (xr(cell_nodes(j0, r))-xr(r0), normal); if (std::abs(height) > std::abs(max_height)) { max_height = height; } } if (max_height > 0) { return -normal; } else { return normal; } } #endif // MESH_NODE_BOUNDARY_HPP