#ifndef MESH_DATA_HPP #define MESH_DATA_HPP #include <Kokkos_Core.hpp> #include <TinyVector.hpp> template <typename M> class MeshData { public: typedef M MeshType; typedef TinyVector<dimension> Rd; static constexpr size_t dimension = MeshType::dimension; static_assert(dimension>0, "dimension must be strictly positive"); static constexpr double inv_dimension = 1./dimension; private: const MeshType& m_mesh; Kokkos::View<Rd**> m_Cjr; Kokkos::View<Rd*> m_xj; Kokkos::View<double*> m_Vj; KOKKOS_INLINE_FUNCTION void _updateCenter() { // Computes vertices isobarycenter if(dimension == 1) { const Kokkos::View<const Rd*> xr = m_mesh.xr(); const Kokkos::View<const unsigned int**>& cell_nodes = m_mesh.connectivity().cellNodes(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ m_xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]); }); } else { const Kokkos::View<const Rd*> xr = m_mesh.xr(); const Kokkos::View<const unsigned int**>& cell_nodes = m_mesh.connectivity().cellNodes(); const Kokkos::View<const unsigned short*>& cell_nb_nodes = m_mesh.connectivity().cellNbNodes(); const Kokkos::View<const double*>& inv_cell_nb_nodes = m_mesh.connectivity().invCellNbNodes(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Rd X = zero; for (int R=0; R<cell_nb_nodes[j]; ++R) { X += xr[cell_nodes(j,R)]; } m_xj[j] = inv_cell_nb_nodes[j]*X; }); } } KOKKOS_INLINE_FUNCTION void _updateVolume() { const Kokkos::View<const unsigned int**>& cell_nodes = m_mesh.connectivity().cellNodes(); const Kokkos::View<const unsigned short*> cell_nb_nodes = m_mesh.connectivity().cellNbNodes(); const Kokkos::View<const Rd*> xr = m_mesh.xr(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ double sum_cjr_xr = 0; for (int R=0; R<cell_nb_nodes[j]; ++R) { sum_cjr_xr += (xr[cell_nodes(j,R)], m_Cjr(j,R)); } m_Vj[j] = inv_dimension * sum_cjr_xr; }); } KOKKOS_INLINE_FUNCTION void _updateCjr() { if(dimension == 1) { // Cjr are constant overtime } static_assert(dimension==1, "only 1d is implemented"); } public: const MeshType& mesh() const { return m_mesh; } const Kokkos::View<const Rd**> Cjr() const { return m_Cjr; } const Kokkos::View<const Rd*> xj() const { return m_xj; } const Kokkos::View<const double*> Vj() const { return m_Vj; } void updateAllData() { this->_updateCenter(); this->_updateVolume(); this->_updateCjr(); } MeshData(const MeshType& mesh) : m_mesh(mesh), m_Cjr("Cjr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()), m_xj("xj", mesh.numberOfCells()), m_Vj("Vj", mesh.numberOfCells()) { if (dimension==1) { // in 1d Cjr are computed once for all Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { m_Cjr(j,0)=-1; m_Cjr(j,1)= 1; }); } this->updateAllData(); } ~MeshData() { ; } }; #endif // MESH_DATA_HPP