Merge pull request #16691 from mschreter/particle_handler_tolerance_i…

…nside_cell

`ParticleHandler`: use tolerance for is_inside_unit_cell() check

doc/news/changes/minor/20240222SchreterBrotz

@@ -0,0 +1,4 @@
+Fixed: In a special case where a particle lies on a vertex, the ParticleHandler detected it as lost.
+We introduce a tolerance for GeometryInfo<dim>::is_inside_cell() inside the ParticleHandler to fix this issue.
+<br>
+(Magdalena Schreter-Fleischhacker, Julian Brotz, 2024/02/22)

include/deal.II/particles/particle_handler.h

@@ -1018,6 +1018,11 @@ namespace Particles
      */
     unsigned int handle;
 
+    /**
+     * Tolerance to be used for GeometryInfo<dim>::is_inside_cell().
+     */
+    const double tolerance_inside_cell = 1e-12;
+
     /**
      * The GridTools::Cache is used to store the information about the
      * vertex_to_cells set and the vertex_to_cell_centers vectors to prevent

source/particles/particle_handler.cc

@@ -1262,7 +1262,8 @@ namespace Particles
         for (const auto &p_unit : reference_locations)
           {
             if (numbers::is_finite(p_unit[0]) &&
-                GeometryInfo<dim>::is_inside_unit_cell(p_unit))
+                GeometryInfo<dim>::is_inside_unit_cell(p_unit,
+                                                       tolerance_inside_cell))
               particle->set_reference_location(p_unit);
             else
               particles_out_of_cell.push_back(particle);
@@ -1384,8 +1385,8 @@ namespace Particles
                                                           real_locations,
                                                           reference_locations);
 
-              if (GeometryInfo<dim>::is_inside_unit_cell(
-                    reference_locations[0]))
+              if (GeometryInfo<dim>::is_inside_unit_cell(reference_locations[0],
+                                                         tolerance_inside_cell))
                 {
                   current_cell = *cell;
                   found_cell   = true;
@@ -1433,7 +1434,7 @@ namespace Particles
                     cell, real_locations, reference_locations);
 
                   if (GeometryInfo<dim>::is_inside_unit_cell(
-                        reference_locations[0]))
+                        reference_locations[0], tolerance_inside_cell))
                     {
                       current_cell = cell;
                       found_cell   = true;
@@ -2431,8 +2432,8 @@ namespace Particles
                         const Point<dim> p_unit =
                           mapping->transform_real_to_unit_cell(
                             child, particle->get_location());
-                        if (GeometryInfo<dim>::is_inside_unit_cell(p_unit,
-                                                                   1e-12))
+                        if (GeometryInfo<dim>::is_inside_unit_cell(
+                              p_unit, tolerance_inside_cell))
                           {
                             found_new_cell = true;
                             particle->set_reference_location(p_unit);

tests/particles/particle_handler_sort_02.cc

@@ -0,0 +1,93 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test if sort_particles_into_subdomains_and_cells() works for a particle
+// located at a vertex.
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/particles/particle_handler.h>
+
+#include "../tests.h"
+
+template <int dim, int spacedim>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim, spacedim> tr(MPI_COMM_WORLD);
+
+  GridGenerator::subdivided_hyper_cube(tr, 20, 0, 4);
+
+  MappingQ<dim, spacedim> mapping(3);
+
+  // both processes create a particle handler with a particle in a
+  // position that is in a ghost cell to the other process
+  Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
+
+  particle_handler.signals.particle_lost.connect(
+    [](const typename Particles::ParticleIterator<dim>         &particle,
+       const typename Triangulation<dim>::active_cell_iterator &cell) {
+      AssertThrow(false, ExcMessage("Particle lost."));
+    });
+
+  std::vector<Point<spacedim>> position(1);
+
+  if (Utilities::MPI::this_mpi_process(tr.get_communicator()) == 0)
+    for (unsigned int i = 0; i < dim; ++i)
+      position[0][i] = 2;
+
+  auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
+    tr, IteratorFilters::LocallyOwnedCell());
+
+  auto global_bounding_boxes =
+    Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box);
+
+  particle_handler.insert_global_particles(position, global_bounding_boxes);
+
+  // Reverse the refinement and check again
+  for (auto cell = tr.begin_active(); cell != tr.end(); ++cell)
+    {
+      if (cell->is_locally_owned())
+        if (cell->center().distance(position[0]) <= 1)
+          cell->set_refine_flag();
+    }
+
+  particle_handler.prepare_for_coarsening_and_refinement();
+  tr.prepare_coarsening_and_refinement();
+  tr.execute_coarsening_and_refinement();
+  particle_handler.unpack_after_coarsening_and_refinement();
+  particle_handler.sort_particles_into_subdomains_and_cells();
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  MPILogInitAll all;
+
+  test<2, 2>();
+}

tests/particles/particle_handler_sort_02.with_mpi=true.with_p4est=true.output

@@ -0,0 +1,2 @@
+
+DEAL:0::OK