diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp index 4351e61b082dd3d892a0eed23671b388b8857d9c..49513c0549421e29b4307140bd9159fce312f7ba 100644 --- a/src/mesh/MeshData.hpp +++ b/src/mesh/MeshData.hpp @@ -124,40 +124,54 @@ class MeshData } }); + Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){ + m_Cjr[jr] = zero; + }); + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j); - for (size_t R=0; R<cell_nodes.length; ++R) { - m_Cjr(j,R) = zero; - } - std::map<unsigned int, unsigned short> node_id_to_local; - for (size_t R=0; R<cell_nodes.length; ++R) { - node_id_to_local[cell_nodes(R)] = R; - } + const auto& cell_faces = m_mesh.connectivity().m_cell_to_face_matrix.rowConst(j); const auto& cell_face_is_reversed = m_mesh.connectivity().m_cell_face_is_reversed.itemValues(j); + for (size_t L=0; L<cell_faces.length; ++L) { const size_t l = cell_faces(L); const auto& face_nodes = m_mesh.connectivity().m_face_to_node_matrix.rowConst(l); + +#warning should this lambda be replaced by a precomputed correspondance? + std::function local_node_number_in_cell + = [&](const size_t& node_number) { + for (size_t i_node=0; i_node<cell_nodes.length; ++i_node) { + if (node_number == cell_nodes(i_node)) { + return i_node; + break; + } + } + return std::numeric_limits<size_t>::max(); + }; + if (cell_face_is_reversed[L]) { for (size_t rl = 0; rl<face_nodes.length; ++rl) { - m_Cjr(j, node_id_to_local[face_nodes(rl)]) -= Nlr(l,rl); + const size_t R = local_node_number_in_cell(face_nodes(rl)); + m_Cjr(j, R) -= Nlr(l,rl); } } else { for (size_t rl = 0; rl<face_nodes.length; ++rl) { - m_Cjr(j, node_id_to_local[face_nodes(rl)]) += Nlr(l, rl); + const size_t R = local_node_number_in_cell(face_nodes(rl)); + m_Cjr(j, R) += Nlr(l,rl); } } } }); const NodeValuePerCell<Rd>& Cjr = m_Cjr; - Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){ - m_ljr[j] = l2Norm(Cjr[j]); + Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){ + m_ljr[jr] = l2Norm(Cjr[jr]); }); const NodeValuePerCell<double>& ljr = m_ljr; - Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){ - m_njr[j] = (1./ljr[j])*Cjr[j]; + Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){ + m_njr[jr] = (1./ljr[jr])*Cjr[jr]; }); } static_assert((dimension<=3), "only 1d, 2d and 3d are implemented");