From: Magdalena Schreter Date: Thu, 22 Feb 2024 22:07:01 +0000 (+0100) Subject: add tolerance for inside cell check X-Git-Tag: relicensing~1^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=449c986928dc059105e0278353adaa386785cc51;p=dealii.git add tolerance for inside cell check --- diff --git a/doc/news/changes/minor/20240222SchreterBrotz b/doc/news/changes/minor/20240222SchreterBrotz new file mode 100644 index 0000000000..8293f645bf --- /dev/null +++ b/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::is_inside_cell() inside the ParticleHandler to fix this issue. +
+(Magdalena Schreter-Fleischhacker, Julian Brotz, 2024/02/22) diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h index 473f6a13b8..4443f26223 100644 --- a/include/deal.II/particles/particle_handler.h +++ b/include/deal.II/particles/particle_handler.h @@ -1018,6 +1018,11 @@ namespace Particles */ unsigned int handle; + /** + * Tolerance to be used for GeometryInfo::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 diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index 0d47946dc6..5824031769 100644 --- a/source/particles/particle_handler.cc +++ b/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::is_inside_unit_cell(p_unit)) + GeometryInfo::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::is_inside_unit_cell( - reference_locations[0])) + if (GeometryInfo::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::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 p_unit = mapping->transform_real_to_unit_cell( child, particle->get_location()); - if (GeometryInfo::is_inside_unit_cell(p_unit, - 1e-12)) + if (GeometryInfo::is_inside_unit_cell( + p_unit, tolerance_inside_cell)) { found_new_cell = true; particle->set_reference_location(p_unit); diff --git a/tests/particles/particle_handler_sort_02.cc b/tests/particles/particle_handler_sort_02.cc new file mode 100644 index 0000000000..72440e1ecd --- /dev/null +++ b/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 + +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::subdivided_hyper_cube(tr, 20, 0, 4); + + MappingQ 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 particle_handler(tr, mapping); + + particle_handler.signals.particle_lost.connect( + [](const typename Particles::ParticleIterator &particle, + const typename Triangulation::active_cell_iterator &cell) { + AssertThrow(false, ExcMessage("Particle lost.")); + }); + + std::vector> 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>(); +} diff --git a/tests/particles/particle_handler_sort_02.with_mpi=true.with_p4est=true.output b/tests/particles/particle_handler_sort_02.with_mpi=true.with_p4est=true.output new file mode 100644 index 0000000000..be8d055f86 --- /dev/null +++ b/tests/particles/particle_handler_sort_02.with_mpi=true.with_p4est=true.output @@ -0,0 +1,2 @@ + +DEAL:0::OK