diff --git a/experimental/Connectivity1D.hpp b/experimental/Connectivity1D.hpp index 6431faae7e538366ef3063058b1c6da849ec9fbc..8af879957dc3cd74cddb65dcb2faac1af6733631 100644 --- a/experimental/Connectivity1D.hpp +++ b/experimental/Connectivity1D.hpp @@ -19,6 +19,7 @@ private: Kokkos::View<unsigned short*> m_cell_nb_nodes; const Kokkos::View<unsigned short*> m_cell_nb_faces; + Kokkos::View<double*> m_inv_cell_nb_nodes; Kokkos::View<unsigned short*> m_node_nb_cells; const Kokkos::View<unsigned short*>& m_face_nb_cells; @@ -72,6 +73,11 @@ public: 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; @@ -111,6 +117,7 @@ public: m_cell_nodes ("cell_nodes", m_number_of_cells), m_cell_faces (m_cell_nodes), m_cell_nb_nodes ("cell_nb_nodes", m_number_of_cells), + m_inv_cell_nb_nodes ("inv_cell_nb_nodes", m_number_of_cells), m_cell_nb_faces (m_cell_nb_nodes), m_node_nb_cells ("node_nb_cells",m_number_of_nodes), m_face_nb_cells (m_node_nb_cells), @@ -126,6 +133,10 @@ public: m_cell_nb_nodes[j] = 2; }); + Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) { + m_inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j]; + }); + Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) { for (int R=0; R<2; ++R) { m_cell_nodes(j,R)=j+R; diff --git a/experimental/MeshData.hpp b/experimental/MeshData.hpp index c2125f5f0b04a47b7b982102949f81cf94b27d90..d4b4d28020e4cfbc5a078a9c68eeeed865f7d053 100644 --- a/experimental/MeshData.hpp +++ b/experimental/MeshData.hpp @@ -11,37 +11,46 @@ public: typedef M MeshType; typedef TinyVector<dimension> Rd; - typedef double R; static constexpr size_t dimension = MeshType::dimension; static_assert(dimension>0, "dimension must be strictly positive"); - static constexpr R inv_dimension = 1./dimension; + static constexpr double inv_dimension = 1./dimension; private: const MeshType& m_mesh; Kokkos::View<Rd**> m_Cjr; Kokkos::View<Rd*> m_xj; - Kokkos::View<R*> m_Vj; + Kokkos::View<double*> m_Vj; - template<size_t Dim> KOKKOS_INLINE_FUNCTION - void _updateCenter(); - - template <> - KOKKOS_INLINE_FUNCTION - void _updateCenter<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)]); - }); + 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; + }); + } } - template<size_t Dim> KOKKOS_INLINE_FUNCTION void _updateVolume() { @@ -54,7 +63,7 @@ private: const Kokkos::View<const Rd*> xr = m_mesh.xr(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - R sum_cjr_xr = 0; + 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)); } @@ -62,38 +71,20 @@ private: }); } - template<> KOKKOS_INLINE_FUNCTION - void _updateVolume<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_Vj[j] - = (xr[cell_nodes(j,1)], m_Cjr(j,1)) - + (xr[cell_nodes(j,0)], m_Cjr(j,0)); - }); + void _updateCjr() { + if(dimension == 1) { + // Cjr are constant overtime + } + static_assert(dimension==1, "only 1d is implemented"); } - template<size_t Dim> - KOKKOS_INLINE_FUNCTION - void _updateCjr(); - - template<> - KOKKOS_INLINE_FUNCTION - void _updateCjr<1>() {} - public: - const MeshType& mesh() const { return m_mesh; } - const Kokkos::View<const Rd**> Cjr() const { return m_Cjr; @@ -104,16 +95,16 @@ public: return m_xj; } - const Kokkos::View<const R*> Vj() const + const Kokkos::View<const double*> Vj() const { return m_Vj; } void updateAllData() { - this->_updateCenter<dimension>(); - this->_updateVolume<dimension>(); - this->_updateCjr<dimension>(); + this->_updateCenter(); + this->_updateVolume(); + this->_updateCjr(); } MeshData(const MeshType& mesh) @@ -122,10 +113,13 @@ public: m_xj("xj", mesh.numberOfCells()), m_Vj("Vj", mesh.numberOfCells()) { - Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - m_Cjr(j,0)=-1; - m_Cjr(j,1)= 1; - }); + 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(); }