From 9e3690e1bb831e7996e94c07bbddf9a6a0566646 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 1 Dec 2021 10:48:00 -0500 Subject: [PATCH] Faster and simpler sort of lost particles --- source/particles/particle_handler.cc | 83 +++++++++++++++++----------- 1 file changed, 52 insertions(+), 31 deletions(-) diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index bdc3b69301..caa254fece 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -1310,17 +1310,23 @@ namespace Particles std::vector neighbor_permutation; + // Reuse these vectors below, but only with a single element. + // Avoid resizing for every particle. + reference_locations.resize(1); + real_locations.resize(1); + // Find the cells that the particles moved to. for (auto &out_particle : particles_out_of_cell) { // make a copy of the current cell, since we will modify the - // variable current_cell in the following but we need the original in + // variable current_cell below, but we need the original in // the case the particle is not found auto current_cell = out_particle->get_surrounding_cell(); - // The cell the particle is in - Point current_reference_position; - bool found_cell = false; + real_locations[0] = out_particle->get_location(); + + // Record if the new cell was found + bool found_cell = false; // Check if the particle is in one of the old cell's neighbors // that are adjacent to the closest vertex @@ -1356,26 +1362,22 @@ namespace Particles // previous cell Most likely we will find the particle in them. for (unsigned int i = 0; i < n_neighbor_cells; ++i) { - try + typename std::set:: + active_cell_iterator>::const_iterator cell = + vertex_to_cells[closest_vertex_index].begin(); + + std::advance(cell, neighbor_permutation[i]); + mapping->transform_points_real_to_unit_cell(*cell, + real_locations, + reference_locations); + + if (GeometryInfo::is_inside_unit_cell( + reference_locations[0])) { - typename std::set:: - active_cell_iterator>::const_iterator - cell = vertex_to_cells[closest_vertex_index].begin(); - - std::advance(cell, neighbor_permutation[i]); - const Point p_unit = - mapping->transform_real_to_unit_cell( - *cell, out_particle->get_location()); - if (GeometryInfo::is_inside_unit_cell(p_unit)) - { - current_cell = *cell; - current_reference_position = p_unit; - found_cell = true; - break; - } + current_cell = *cell; + found_cell = true; + break; } - catch (typename Mapping::ExcTransformationFailed &) - {} } if (!found_cell) @@ -1383,14 +1385,33 @@ namespace Particles // The particle is not in a neighbor of the old cell. // Look for the new cell in the whole local domain. // This case is rare. - const std::pair:: - active_cell_iterator, - Point> - current_cell_and_position = - GridTools::find_active_cell_around_point<>( - *triangulation_cache, out_particle->get_location()); - current_cell = current_cell_and_position.first; - current_reference_position = current_cell_and_position.second; + std::vector, unsigned int>> + closest_vertex_in_domain; + triangulation_cache->get_used_vertices_rtree().query( + boost::geometry::index::nearest(out_particle->get_location(), + 1), + std::back_inserter(closest_vertex_in_domain)); + + // We should have one and only one result + AssertDimension(closest_vertex_in_domain.size(), 1); + const unsigned int closest_vertex_index_in_domain = + closest_vertex_in_domain[0].second; + + // Search all of the cells adjacent to the closest vertex of the + // domain. Most likely we will find the particle in them. + for (const auto &cell : + vertex_to_cells[closest_vertex_index_in_domain]) + { + mapping->transform_points_real_to_unit_cell( + cell, real_locations, reference_locations); + + if (GeometryInfo::is_inside_unit_cell( + reference_locations[0])) + { + current_cell = cell; + break; + } + } if (current_cell.state() != IteratorState::valid) { @@ -1405,7 +1426,7 @@ namespace Particles // If we are here, we found a cell and reference position for this // particle - out_particle->set_reference_location(current_reference_position); + out_particle->set_reference_location(reference_locations[0]); // Reinsert the particle into our domain if we own its cell. // Mark it for MPI transfer otherwise -- 2.39.5