From c93635588ee1d33ba0753bbfcdc9667bb3aa43db Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 2 Oct 2020 12:15:32 +0200 Subject: [PATCH] Use multiple point evaluation in ParticleHandler::sort_particles_... --- source/particles/particle_handler.cc | 49 ++++++++++++++++------------ 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index 956a3c08fb..2376c51dc0 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -925,29 +925,38 @@ namespace Particles particles_out_of_cell.reserve(n_locally_owned_particles()); // Now update the reference locations of the moved particles - for (particle_iterator it = begin(); it != end(); ++it) + std::vector> real_locations; + std::vector> reference_locations; + for (auto particle = begin(); particle != end();) { - const typename Triangulation::cell_iterator cell = - it->get_surrounding_cell(*triangulation); - - try + const auto cell = particle->get_surrounding_cell(*triangulation); + real_locations.clear(); + + // Since we might also work on artificial cells when we initialize the + // particles on a remote processor, we cannot use the + // particles_in_cell method. Thus, We instead simply go through the + // particles and check if the next one belongs to the same cell as the + // current one. + for (auto it = particle; + it != end() && it->get_surrounding_cell(*triangulation) == cell; + ++it) + real_locations.push_back(it->get_location()); + + reference_locations.resize(real_locations.size()); + ArrayView> reference(reference_locations.data(), + reference_locations.size()); + mapping->transform_points_real_to_unit_cell(cell, + real_locations, + reference); + + for (const auto &p_unit : reference_locations) { - const Point p_unit = - mapping->transform_real_to_unit_cell(cell, it->get_location()); - if (GeometryInfo::is_inside_unit_cell(p_unit)) - { - it->set_reference_location(p_unit); - } + if (p_unit[0] == std::numeric_limits::infinity() || + !GeometryInfo::is_inside_unit_cell(p_unit)) + particles_out_of_cell.push_back(particle); else - { - // The particle has left the cell - particles_out_of_cell.push_back(it); - } - } - catch (typename Mapping::ExcTransformationFailed &) - { - // The particle has left the cell - particles_out_of_cell.push_back(it); + particle->set_reference_location(p_unit); + ++particle; } } -- 2.39.5