From: Daniel Arndt Date: Fri, 29 Nov 2019 20:00:23 +0000 (-0500) Subject: Allow using serial triangulations with ParticleHandler X-Git-Tag: v9.2.0-rc1~831^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fb7ea534638ddd1ace22023ddd3c17dcfd25c71f;p=dealii.git Allow using serial triangulations with ParticleHandler --- diff --git a/doc/news/changes/minor/20191129DanielArndt b/doc/news/changes/minor/20191129DanielArndt new file mode 100644 index 0000000000..8346118c15 --- /dev/null +++ b/doc/news/changes/minor/20191129DanielArndt @@ -0,0 +1,3 @@ +Improved: Allow using ParticleHandler with serial triangulations. +
+(Daniel Arndt, 2019/11/29) diff --git a/include/deal.II/particles/generators.h b/include/deal.II/particles/generators.h index 272967c781..c1d84c0d87 100644 --- a/include/deal.II/particles/generators.h +++ b/include/deal.II/particles/generators.h @@ -30,8 +30,6 @@ DEAL_II_NAMESPACE_OPEN -#ifdef DEAL_II_WITH_P4EST - namespace Particles { /** @@ -165,8 +163,6 @@ namespace Particles } // namespace Generators } // namespace Particles -#endif // DEAL_II_WITH_P4EST - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h index 9988841514..d5a5f10c1a 100644 --- a/include/deal.II/particles/particle_handler.h +++ b/include/deal.II/particles/particle_handler.h @@ -34,8 +34,6 @@ DEAL_II_NAMESPACE_OPEN -#ifdef DEAL_II_WITH_P4EST - namespace Particles { /** @@ -77,10 +75,9 @@ namespace Particles * This constructor is equivalent to calling the default constructor and * the initialize function. */ - ParticleHandler( - const parallel::distributed::Triangulation &tria, - const Mapping & mapping, - const unsigned int n_properties = 0); + ParticleHandler(const Triangulation &tria, + const Mapping & mapping, + const unsigned int n_properties = 0); /** * Destructor. @@ -89,13 +86,13 @@ namespace Particles /** * Initialize the particle handler. This function does not clear the - * internal data structures, it just sets the connections to the - * MPI communicator and the triangulation. + * internal data structures, it just sets the triangulation and the mapping + * to be used. */ void - initialize(const parallel::distributed::Triangulation &tria, - const Mapping &mapping, - const unsigned int n_properties = 0); + initialize(const Triangulation &tria, + const Mapping & mapping, + const unsigned int n_properties = 0); /** * Clear all particle related data. @@ -367,7 +364,7 @@ namespace Particles /** * Address of the triangulation to work on. */ - SmartPointer, + SmartPointer, ParticleHandler> triangulation; @@ -464,7 +461,7 @@ namespace Particles */ unsigned int handle; -# ifdef DEAL_II_WITH_MPI +#ifdef DEAL_II_WITH_MPI /** * Transfer particles that have crossed subdomain boundaries to other * processors. @@ -500,7 +497,7 @@ namespace Particles types::subdomain_id, std::vector< typename Triangulation::active_cell_iterator>>()); -# endif +#endif /** * Called by listener functions from Triangulation for every cell @@ -546,8 +543,6 @@ namespace Particles } } // namespace Particles -#endif // DEAL_II_WITH_P4EST - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/particles/generators.cc b/source/particles/generators.cc index 9b35f33ab7..9619fa0e25 100644 --- a/source/particles/generators.cc +++ b/source/particles/generators.cc @@ -25,8 +25,6 @@ DEAL_II_NAMESPACE_OPEN -#ifdef DEAL_II_WITH_P4EST - namespace Particles { namespace Generators @@ -114,6 +112,7 @@ namespace Particles { types::particle_index particle_index = 0; +#ifdef DEAL_II_WITH_MPI if (const auto tria = dynamic_cast *>( &triangulation)) @@ -131,6 +130,7 @@ namespace Particles MPI_SUM, tria->get_communicator()); } +#endif for (const auto &cell : triangulation.active_cell_iterators()) { @@ -284,6 +284,7 @@ namespace Particles // the weights of all processes with a lower rank double local_start_weight = 0.0; +#ifdef DEAL_II_WITH_MPI if (const auto tria = dynamic_cast *>( &triangulation)) @@ -295,6 +296,7 @@ namespace Particles MPI_SUM, tria->get_communicator()); } +#endif // Calculate start id and number of local particles start_particle_id = @@ -392,14 +394,6 @@ namespace Particles } // namespace Generators } // namespace Particles -#endif // DEAL_II_WITH_P4EST - -DEAL_II_NAMESPACE_CLOSE - -DEAL_II_NAMESPACE_OPEN - -#ifdef DEAL_II_WITH_P4EST -# include "generators.inst" -#endif +#include "generators.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index 65aa37fbff..f64f192fcb 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -24,8 +24,6 @@ DEAL_II_NAMESPACE_OPEN -#ifdef DEAL_II_WITH_P4EST - namespace Particles { namespace @@ -106,9 +104,9 @@ namespace Particles template ParticleHandler::ParticleHandler( - const parallel::distributed::Triangulation &triangulation, - const Mapping & mapping, - const unsigned int n_properties) + const Triangulation &triangulation, + const Mapping & mapping, + const unsigned int n_properties) : triangulation(&triangulation, typeid(*this).name()) , mapping(&mapping, typeid(*this).name()) , particles() @@ -128,10 +126,9 @@ namespace Particles template void ParticleHandler::initialize( - const parallel::distributed::Triangulation - & new_triangulation, - const Mapping &new_mapping, - const unsigned int n_properties) + const Triangulation &new_triangulation, + const Mapping & new_mapping, + const unsigned int n_properties) { triangulation = &new_triangulation; mapping = &new_mapping; @@ -189,16 +186,26 @@ namespace Particles std::max(local_max_particles_per_cell, current_particles_per_cell); } - global_number_of_particles = - dealii::Utilities::MPI::sum(particles.size(), - triangulation->get_communicator()); - next_free_particle_index = - dealii::Utilities::MPI::max(locally_highest_index, - triangulation->get_communicator()) + - 1; - global_max_particles_per_cell = - dealii::Utilities::MPI::max(local_max_particles_per_cell, - triangulation->get_communicator()); + if (const auto parallel_triangulation = + dynamic_cast *>( + &*triangulation)) + { + global_number_of_particles = dealii::Utilities::MPI::sum( + particles.size(), parallel_triangulation->get_communicator()); + next_free_particle_index = + dealii::Utilities::MPI::max( + locally_highest_index, parallel_triangulation->get_communicator()) + + 1; + global_max_particles_per_cell = dealii::Utilities::MPI::max( + local_max_particles_per_cell, + parallel_triangulation->get_communicator()); + } + else + { + global_number_of_particles = particles.size(); + next_free_particle_index = locally_highest_index + 1; + global_max_particles_per_cell = local_max_particles_per_cell; + } } @@ -380,17 +387,22 @@ namespace Particles get_next_free_particle_index(); types::particle_index local_start_index = 0; -# ifdef DEAL_II_WITH_MPI - types::particle_index particles_to_add_locally = positions.size(); - const int ierr = MPI_Scan(&particles_to_add_locally, - &local_start_index, - 1, - DEAL_II_PARTICLE_INDEX_MPI_TYPE, - MPI_SUM, - triangulation->get_communicator()); - AssertThrowMPI(ierr); - local_start_index -= particles_to_add_locally; -# endif +#ifdef DEAL_II_WITH_MPI + if (const auto parallel_triangulation = + dynamic_cast *>( + &*triangulation)) + { + types::particle_index particles_to_add_locally = positions.size(); + const int ierr = MPI_Scan(&particles_to_add_locally, + &local_start_index, + 1, + DEAL_II_PARTICLE_INDEX_MPI_TYPE, + MPI_SUM, + parallel_triangulation->get_communicator()); + AssertThrowMPI(ierr); + local_start_index -= particles_to_add_locally; + } +#endif local_start_index += local_next_particle_index; @@ -598,8 +610,11 @@ namespace Particles sorted_particles.reserve( static_cast(particles_out_of_cell.size() * 1.25)); - const std::set ghost_owners = - triangulation->ghost_owners(); + std::set ghost_owners; + if (const auto parallel_triangulation = + dynamic_cast *>( + &*triangulation)) + ghost_owners = parallel_triangulation->ghost_owners(); for (const auto ghost_owner : ghost_owners) moved_particles[ghost_owner].reserve( @@ -743,11 +758,18 @@ namespace Particles sorted_particles_map; // Exchange particles between processors if we have more than one process -# ifdef DEAL_II_WITH_MPI - if (dealii::Utilities::MPI::n_mpi_processes( - triangulation->get_communicator()) > 1) - send_recv_particles(moved_particles, sorted_particles_map, moved_cells); -# endif +#ifdef DEAL_II_WITH_MPI + if (const auto parallel_triangulation = + dynamic_cast *>( + &*triangulation)) + { + if (dealii::Utilities::MPI::n_mpi_processes( + parallel_triangulation->get_communicator()) > 1) + send_recv_particles(moved_particles, + sorted_particles_map, + moved_cells); + } +#endif sorted_particles_map.insert(sorted_particles.begin(), sorted_particles.end()); @@ -766,11 +788,19 @@ namespace Particles ParticleHandler::exchange_ghost_particles() { // Nothing to do in serial computations - if (dealii::Utilities::MPI::n_mpi_processes( - triangulation->get_communicator()) == 1) + const auto parallel_triangulation = + dynamic_cast *>( + &*triangulation); + if (parallel_triangulation != nullptr) + { + if (dealii::Utilities::MPI::n_mpi_processes( + parallel_triangulation->get_communicator()) == 1) + return; + } + else return; -# ifdef DEAL_II_WITH_MPI +#ifdef DEAL_II_WITH_MPI // First clear the current ghost_particle information ghost_particles.clear(); @@ -778,7 +808,7 @@ namespace Particles ghost_particles_by_domain; const std::set ghost_owners = - triangulation->ghost_owners(); + parallel_triangulation->ghost_owners(); for (const auto ghost_owner : ghost_owners) ghost_particles_by_domain[ghost_owner].reserve( static_cast::size_type>( @@ -827,12 +857,12 @@ namespace Particles } send_recv_particles(ghost_particles_by_domain, ghost_particles); -# endif +#endif } -# ifdef DEAL_II_WITH_MPI +#ifdef DEAL_II_WITH_MPI template void ParticleHandler::send_recv_particles( @@ -845,9 +875,17 @@ namespace Particles std::vector::active_cell_iterator>> &send_cells) { + const auto parallel_triangulation = + dynamic_cast *>( + &*triangulation); + Assert( + parallel_triangulation, + ExcMessage( + "This function is only implemented for parallel::Triangulation objects.")); + // Determine the communication pattern const std::set ghost_owners = - triangulation->ghost_owners(); + parallel_triangulation->ghost_owners(); const std::vector neighbors(ghost_owners.begin(), ghost_owners.end()); const unsigned int n_neighbors = neighbors.size(); @@ -948,7 +986,7 @@ namespace Particles MPI_UNSIGNED, neighbors[i], mpi_tag, - triangulation->get_communicator(), + parallel_triangulation->get_communicator(), &(n_requests[2 * i])); AssertThrowMPI(ierr); } @@ -959,7 +997,7 @@ namespace Particles MPI_UNSIGNED, neighbors[i], mpi_tag, - triangulation->get_communicator(), + parallel_triangulation->get_communicator(), &(n_requests[2 * i + 1])); AssertThrowMPI(ierr); } @@ -991,13 +1029,14 @@ namespace Particles for (unsigned int i = 0; i < n_neighbors; ++i) if (n_recv_data[i] > 0) { - const int ierr = MPI_Irecv(&(recv_data[recv_offsets[i]]), - n_recv_data[i], - MPI_CHAR, - neighbors[i], - mpi_tag, - triangulation->get_communicator(), - &(requests[send_ops])); + const int ierr = + MPI_Irecv(&(recv_data[recv_offsets[i]]), + n_recv_data[i], + MPI_CHAR, + neighbors[i], + mpi_tag, + parallel_triangulation->get_communicator(), + &(requests[send_ops])); AssertThrowMPI(ierr); send_ops++; } @@ -1005,13 +1044,14 @@ namespace Particles for (unsigned int i = 0; i < n_neighbors; ++i) if (n_send_data[i] > 0) { - const int ierr = MPI_Isend(&(send_data[send_offsets[i]]), - n_send_data[i], - MPI_CHAR, - neighbors[i], - mpi_tag, - triangulation->get_communicator(), - &(requests[send_ops + recv_ops])); + const int ierr = + MPI_Isend(&(send_data[send_offsets[i]]), + n_send_data[i], + MPI_CHAR, + neighbors[i], + mpi_tag, + parallel_triangulation->get_communicator(), + &(requests[send_ops + recv_ops])); AssertThrowMPI(ierr); recv_ops++; } @@ -1053,7 +1093,7 @@ namespace Particles "The amount of data that was read into new particles " "does not match the amount of data sent around.")); } -# endif +#endif @@ -1079,8 +1119,12 @@ namespace Particles parallel::distributed::Triangulation *non_const_triangulation = const_cast *>( - &(*triangulation)); + dynamic_cast + *>(&(*triangulation))); + Assert(non_const_triangulation != nullptr, dealii::ExcNotImplemented()); + +#ifdef DEAL_II_WITH_P4EST // Only save and load particles if there are any, we might get here for // example if somebody created a ParticleHandler but generated 0 particles. update_cached_numbers(); @@ -1098,6 +1142,7 @@ namespace Particles handle = non_const_triangulation->register_data_attach( callback_function, /*returns_variable_size_data=*/true); } +#endif } @@ -1114,8 +1159,12 @@ namespace Particles parallel::distributed::Triangulation *non_const_triangulation = const_cast *>( - &(*triangulation)); + dynamic_cast + *>(&(*triangulation))); + Assert(non_const_triangulation != nullptr, dealii::ExcNotImplemented()); + +#ifdef DEAL_II_WITH_P4EST // If we are resuming from a checkpoint, we first have to register the // store function again, to set the triangulation in the same state as // before the serialization. Only by this it knows how to deserialize the @@ -1155,6 +1204,9 @@ namespace Particles handle = numbers::invalid_unsigned_int; update_cached_numbers(); } +#else + (void)serialization; +#endif } @@ -1279,19 +1331,19 @@ namespace Particles // particles. This is a C++11 function, but not all compilers // that report a -std=c++11 (like gcc 4.6) implement it, so // require C++14 instead. -# ifdef DEAL_II_WITH_CXX14 +#ifdef DEAL_II_WITH_CXX14 position_hint = particles.emplace_hint(position_hint, std::make_pair(cell->level(), cell->index()), std::move(particle)); -# else +#else position_hint = particles.insert(position_hint, std::make_pair(std::make_pair(cell->level(), cell->index()), std::move(particle))); -# endif +#endif // Move the hint position forward by one, i.e., for the next // particle. The 'hint' position will thus be right after the // one just inserted. @@ -1315,19 +1367,19 @@ namespace Particles // particles. This is a C++11 function, but not all compilers // that report a -std=c++11 (like gcc 4.6) implement it, so // require C++14 instead. -# ifdef DEAL_II_WITH_CXX14 +#ifdef DEAL_II_WITH_CXX14 position_hint = particles.emplace_hint(position_hint, std::make_pair(cell->level(), cell->index()), std::move(particle)); -# else +#else position_hint = particles.insert(position_hint, std::make_pair(std::make_pair(cell->level(), cell->index()), std::move(particle))); -# endif +#endif // Move the hint position forward by one, i.e., for the next // particle. The 'hint' position will thus be right after the // one just inserted. @@ -1374,19 +1426,19 @@ namespace Particles // but not all compilers that report a -std=c++11 // (like gcc 4.6) implement it, so require C++14 // instead. -# ifdef DEAL_II_WITH_CXX14 +#ifdef DEAL_II_WITH_CXX14 position_hints[child_index] = particles.emplace_hint( position_hints[child_index], std::make_pair(child->level(), child->index()), std::move(particle)); -# else +#else position_hints[child_index] = particles.insert( position_hints[child_index], std::make_pair(std::make_pair(child->level(), child->index()), std::move(particle))); -# endif +#endif // Move the hint position forward by one, i.e., for // the next particle. The 'hint' position will thus // be right after the one just inserted. @@ -1408,14 +1460,6 @@ namespace Particles } } // namespace Particles -#endif - -DEAL_II_NAMESPACE_CLOSE - -DEAL_II_NAMESPACE_OPEN - -#ifdef DEAL_II_WITH_P4EST -# include "particle_handler.inst" -#endif +#include "particle_handler.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/tests/particles/particle_handler_serial_01.cc b/tests/particles/particle_handler_serial_01.cc new file mode 100644 index 0000000000..ff93d2e881 --- /dev/null +++ b/tests/particles/particle_handler_serial_01.cc @@ -0,0 +1,102 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2018 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. +// +// --------------------------------------------------------------------- + + + +// check the creation and destruction of particle within the particle handler +// class for a serial triangulation. + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + Triangulation tr; + + GridGenerator::hyper_cube(tr); + MappingQ mapping(1); + + Particles::ParticleHandler particle_handler(tr, mapping); + + + Point position; + position(0) = 0.3; + if (spacedim > 1) + position(1) = 0.5; + if (spacedim > 2) + position(2) = 0.7; + + Point reference_position; + reference_position(0) = 0.2; + if (dim > 1) + reference_position(1) = 0.4; + if (dim > 2) + reference_position(2) = 0.6; + + Particles::Particle particle(position, reference_position, 7); + deallog << "Particle location: " << particle.get_location() << std::endl; + + + std::pair:: + active_cell_iterator, + Point> + cell_position = + GridTools::find_active_cell_around_point(mapping, + tr, + particle.get_location()); + + particle_handler.insert_particle(particle, cell_position.first); + particle_handler.update_cached_numbers(); + + deallog << "Particle number: " << particle_handler.n_global_particles() + << std::endl; + + for (const auto &particle : particle_handler) + { + deallog << "Particle location: " << particle.get_location() << std::endl; + deallog << "Particle reference location: " + << particle.get_reference_location() << std::endl; + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("2d/2d"); + test<2, 2>(); + deallog.pop(); + deallog.push("2d/3d"); + test<2, 3>(); + deallog.pop(); + deallog.push("3d/3d"); + test<3, 3>(); + deallog.pop(); +} diff --git a/tests/particles/particle_handler_serial_01.output b/tests/particles/particle_handler_serial_01.output new file mode 100644 index 0000000000..c8dd406e46 --- /dev/null +++ b/tests/particles/particle_handler_serial_01.output @@ -0,0 +1,16 @@ + +DEAL:2d/2d::Particle location: 0.300000 0.500000 +DEAL:2d/2d::Particle number: 1 +DEAL:2d/2d::Particle location: 0.300000 0.500000 +DEAL:2d/2d::Particle reference location: 0.200000 0.400000 +DEAL:2d/2d::OK +DEAL:2d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:2d/3d::Particle number: 1 +DEAL:2d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:2d/3d::Particle reference location: 0.200000 0.400000 +DEAL:2d/3d::OK +DEAL:3d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:3d/3d::Particle number: 1 +DEAL:3d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:3d/3d::Particle reference location: 0.200000 0.400000 0.600000 +DEAL:3d/3d::OK diff --git a/tests/particles/particle_handler_serial_02.cc b/tests/particles/particle_handler_serial_02.cc new file mode 100644 index 0000000000..6ff66e7207 --- /dev/null +++ b/tests/particles/particle_handler_serial_02.cc @@ -0,0 +1,122 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2018 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. +// +// --------------------------------------------------------------------- + + + +// check the cached numbers inside of the particle handler class if +// particles are distributed in different cells for a serial triangulation. + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + { + Triangulation tr; + + GridGenerator::hyper_cube(tr); + tr.refine_global(1); + MappingQ mapping(1); + + Particles::ParticleHandler particle_handler(tr, mapping); + + + Point position; + position(0) = 0.3; + if (spacedim > 1) + position(1) = 0.5; + if (spacedim > 2) + position(2) = 0.7; + + Point reference_position; + reference_position(0) = 0.2; + if (dim > 1) + reference_position(1) = 0.4; + if (dim > 2) + reference_position(2) = 0.6; + + Particles::Particle particle(position, + reference_position, + 7); + deallog << "Particle location: " << particle.get_location() << std::endl; + + + auto cell_position = + GridTools::find_active_cell_around_point(mapping, + tr, + particle.get_location()); + + particle_handler.insert_particle(particle, cell_position.first); + particle_handler.insert_particle(particle, cell_position.first); + + position(0) = 0.7; + Particles::Particle particle2(position, + reference_position, + 9); + + cell_position = + GridTools::find_active_cell_around_point(mapping, + tr, + particle2.get_location()); + particle_handler.insert_particle(particle2, cell_position.first); + + particle_handler.update_cached_numbers(); + + deallog << "Particle number: " << particle_handler.n_global_particles() + << std::endl; + deallog << "Next free particle index: " + << particle_handler.get_next_free_particle_index() << std::endl; + deallog << "Max particles per cell: " + << particle_handler.n_global_max_particles_per_cell() << std::endl; + + for (const auto &particle : particle_handler) + { + deallog << "Particle location: " << particle.get_location() + << std::endl; + deallog << "Particle reference location: " + << particle.get_reference_location() << std::endl; + } + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("2d/2d"); + test<2, 2>(); + deallog.pop(); + deallog.push("2d/3d"); + test<2, 3>(); + deallog.pop(); + deallog.push("3d/3d"); + test<3, 3>(); + deallog.pop(); +} diff --git a/tests/particles/particle_handler_serial_02.output b/tests/particles/particle_handler_serial_02.output new file mode 100644 index 0000000000..6c118fe038 --- /dev/null +++ b/tests/particles/particle_handler_serial_02.output @@ -0,0 +1,34 @@ + +DEAL:2d/2d::Particle location: 0.300000 0.500000 +DEAL:2d/2d::Particle number: 3 +DEAL:2d/2d::Next free particle index: 10 +DEAL:2d/2d::Max particles per cell: 2 +DEAL:2d/2d::Particle location: 0.300000 0.500000 +DEAL:2d/2d::Particle reference location: 0.200000 0.400000 +DEAL:2d/2d::Particle location: 0.300000 0.500000 +DEAL:2d/2d::Particle reference location: 0.200000 0.400000 +DEAL:2d/2d::Particle location: 0.700000 0.500000 +DEAL:2d/2d::Particle reference location: 0.200000 0.400000 +DEAL:2d/2d::OK +DEAL:2d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:2d/3d::Particle number: 3 +DEAL:2d/3d::Next free particle index: 10 +DEAL:2d/3d::Max particles per cell: 2 +DEAL:2d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:2d/3d::Particle reference location: 0.200000 0.400000 +DEAL:2d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:2d/3d::Particle reference location: 0.200000 0.400000 +DEAL:2d/3d::Particle location: 0.700000 0.500000 0.700000 +DEAL:2d/3d::Particle reference location: 0.200000 0.400000 +DEAL:2d/3d::OK +DEAL:3d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:3d/3d::Particle number: 3 +DEAL:3d/3d::Next free particle index: 10 +DEAL:3d/3d::Max particles per cell: 2 +DEAL:3d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:3d/3d::Particle reference location: 0.200000 0.400000 0.600000 +DEAL:3d/3d::Particle location: 0.300000 0.500000 0.700000 +DEAL:3d/3d::Particle reference location: 0.200000 0.400000 0.600000 +DEAL:3d/3d::Particle location: 0.700000 0.500000 0.700000 +DEAL:3d/3d::Particle reference location: 0.200000 0.400000 0.600000 +DEAL:3d/3d::OK diff --git a/tests/particles/particle_handler_serial_03.cc b/tests/particles/particle_handler_serial_03.cc new file mode 100644 index 0000000000..c154384397 --- /dev/null +++ b/tests/particles/particle_handler_serial_03.cc @@ -0,0 +1,117 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2018 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. +// +// --------------------------------------------------------------------- + + + +// check the particle sort within the particle handler class in serial models. +// Same as particle_handler_03 for a serial triangulation. + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + { + Triangulation tr; + + GridGenerator::hyper_cube(tr); + tr.refine_global(1); + MappingQ mapping(1); + + Particles::ParticleHandler particle_handler(tr, mapping); + + std::vector> position(2); + std::vector> reference_position(2); + + for (unsigned int i = 0; i < dim; ++i) + { + position[0](i) = 0.25; + position[1](i) = 0.75; + } + + Particles::Particle particle1(position[0], + reference_position[0], + 0); + Particles::Particle particle2(position[1], + reference_position[1], + 1); + + typename Triangulation::active_cell_iterator cell1(&tr, + 1, + 0); + typename Triangulation::active_cell_iterator cell2(&tr, + 1, + 0); + + particle_handler.insert_particle(particle1, cell1); + particle_handler.insert_particle(particle2, cell2); + + for (const auto &particle : particle_handler) + deallog << "Before sort particle id " << particle.get_id() + << " is in cell " << particle.get_surrounding_cell(tr) + << std::endl; + + particle_handler.sort_particles_into_subdomains_and_cells(); + + for (const auto &particle : particle_handler) + deallog << "After sort particle id " << particle.get_id() + << " is in cell " << particle.get_surrounding_cell(tr) + << std::endl; + + // Move all points up by 0.5. This will change cell for particle 1, and will + // move particle 2 out of the domain. Note that we need to change the + // coordinate dim-1 despite having a spacedim point. + Point shift; + shift(dim - 1) = 0.5; + for (auto &particle : particle_handler) + particle.set_location(particle.get_location() + shift); + + particle_handler.sort_particles_into_subdomains_and_cells(); + for (const auto &particle : particle_handler) + deallog << "After shift particle id " << particle.get_id() + << " is in cell " << particle.get_surrounding_cell(tr) + << std::endl; + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("2d/2d"); + test<2, 2>(); + deallog.pop(); + deallog.push("2d/3d"); + test<2, 3>(); + deallog.pop(); + deallog.push("3d/3d"); + test<3, 3>(); + deallog.pop(); +} diff --git a/tests/particles/particle_handler_serial_03.output b/tests/particles/particle_handler_serial_03.output new file mode 100644 index 0000000000..c63087437d --- /dev/null +++ b/tests/particles/particle_handler_serial_03.output @@ -0,0 +1,19 @@ + +DEAL:2d/2d::Before sort particle id 0 is in cell 1.0 +DEAL:2d/2d::Before sort particle id 1 is in cell 1.0 +DEAL:2d/2d::After sort particle id 0 is in cell 1.0 +DEAL:2d/2d::After sort particle id 1 is in cell 1.3 +DEAL:2d/2d::After shift particle id 0 is in cell 1.2 +DEAL:2d/2d::OK +DEAL:2d/3d::Before sort particle id 0 is in cell 1.0 +DEAL:2d/3d::Before sort particle id 1 is in cell 1.0 +DEAL:2d/3d::After sort particle id 0 is in cell 1.0 +DEAL:2d/3d::After sort particle id 1 is in cell 1.3 +DEAL:2d/3d::After shift particle id 0 is in cell 1.2 +DEAL:2d/3d::OK +DEAL:3d/3d::Before sort particle id 0 is in cell 1.0 +DEAL:3d/3d::Before sort particle id 1 is in cell 1.0 +DEAL:3d/3d::After sort particle id 0 is in cell 1.0 +DEAL:3d/3d::After sort particle id 1 is in cell 1.7 +DEAL:3d/3d::After shift particle id 0 is in cell 1.4 +DEAL:3d/3d::OK diff --git a/tests/particles/particle_handler_serial_07.cc b/tests/particles/particle_handler_serial_07.cc new file mode 100644 index 0000000000..b555d69a33 --- /dev/null +++ b/tests/particles/particle_handler_serial_07.cc @@ -0,0 +1,89 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 2018 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. +// +// --------------------------------------------------------------------- + + + +// check the creation and destruction of particle within the particle handler +// class for a serial triangulation. + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + { + Triangulation tr; + + GridGenerator::hyper_cube(tr); + tr.refine_global(2); + + MappingQ mapping(1); + + Particles::ParticleHandler particle_handler(tr, mapping); + + std::vector> points(10); + for (unsigned int i = 0; i < 10; ++i) + { + const double coordinate = static_cast(i) / 10.0; + for (unsigned int j = 0; j < spacedim; ++j) + points[i][j] = 0.05 + coordinate; + } + + particle_handler.insert_particles(points); + particle_handler.update_cached_numbers(); + + deallog << "Particle number: " << particle_handler.n_global_particles() + << std::endl; + + for (auto particle = particle_handler.begin(); + particle != particle_handler.end(); + ++particle) + deallog << "Particle id " << particle->get_id() << " is in cell " + << particle->get_surrounding_cell(tr) << std::endl + << " at location " << particle->get_location() << std::endl + << " at reference location " + << particle->get_reference_location() << std::endl; + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + initlog(); + + deallog.push("2d/2d"); + test<2, 2>(); + deallog.pop(); + deallog.push("2d/3d"); + test<2, 3>(); + deallog.pop(); + deallog.push("3d/3d"); + test<3, 3>(); + deallog.pop(); +} diff --git a/tests/particles/particle_handler_serial_07.output b/tests/particles/particle_handler_serial_07.output new file mode 100644 index 0000000000..066bf60dcc --- /dev/null +++ b/tests/particles/particle_handler_serial_07.output @@ -0,0 +1,97 @@ + +DEAL:2d/2d::Particle number: 10 +DEAL:2d/2d::Particle id 3 is in cell 2.0 +DEAL:2d/2d:: at location 0.250000 0.250000 +DEAL:2d/2d:: at reference location 1.00000 1.00000 +DEAL:2d/2d::Particle id 2 is in cell 2.0 +DEAL:2d/2d:: at location 0.150000 0.150000 +DEAL:2d/2d:: at reference location 0.600000 0.600000 +DEAL:2d/2d::Particle id 1 is in cell 2.0 +DEAL:2d/2d:: at location 0.0500000 0.0500000 +DEAL:2d/2d:: at reference location 0.200000 0.200000 +DEAL:2d/2d::Particle id 5 is in cell 2.3 +DEAL:2d/2d:: at location 0.450000 0.450000 +DEAL:2d/2d:: at reference location 0.800000 0.800000 +DEAL:2d/2d::Particle id 4 is in cell 2.3 +DEAL:2d/2d:: at location 0.350000 0.350000 +DEAL:2d/2d:: at reference location 0.400000 0.400000 +DEAL:2d/2d::Particle id 8 is in cell 2.12 +DEAL:2d/2d:: at location 0.750000 0.750000 +DEAL:2d/2d:: at reference location 1.00000 1.00000 +DEAL:2d/2d::Particle id 7 is in cell 2.12 +DEAL:2d/2d:: at location 0.650000 0.650000 +DEAL:2d/2d:: at reference location 0.600000 0.600000 +DEAL:2d/2d::Particle id 6 is in cell 2.12 +DEAL:2d/2d:: at location 0.550000 0.550000 +DEAL:2d/2d:: at reference location 0.200000 0.200000 +DEAL:2d/2d::Particle id 10 is in cell 2.15 +DEAL:2d/2d:: at location 0.950000 0.950000 +DEAL:2d/2d:: at reference location 0.800000 0.800000 +DEAL:2d/2d::Particle id 9 is in cell 2.15 +DEAL:2d/2d:: at location 0.850000 0.850000 +DEAL:2d/2d:: at reference location 0.400000 0.400000 +DEAL:2d/2d::OK +DEAL:2d/3d::Particle number: 10 +DEAL:2d/3d::Particle id 3 is in cell 2.0 +DEAL:2d/3d:: at location 0.250000 0.250000 0.250000 +DEAL:2d/3d:: at reference location 1.00000 1.00000 +DEAL:2d/3d::Particle id 2 is in cell 2.0 +DEAL:2d/3d:: at location 0.150000 0.150000 0.150000 +DEAL:2d/3d:: at reference location 0.600000 0.600000 +DEAL:2d/3d::Particle id 1 is in cell 2.0 +DEAL:2d/3d:: at location 0.0500000 0.0500000 0.0500000 +DEAL:2d/3d:: at reference location 0.200000 0.200000 +DEAL:2d/3d::Particle id 5 is in cell 2.3 +DEAL:2d/3d:: at location 0.450000 0.450000 0.450000 +DEAL:2d/3d:: at reference location 0.800000 0.800000 +DEAL:2d/3d::Particle id 4 is in cell 2.3 +DEAL:2d/3d:: at location 0.350000 0.350000 0.350000 +DEAL:2d/3d:: at reference location 0.400000 0.400000 +DEAL:2d/3d::Particle id 8 is in cell 2.12 +DEAL:2d/3d:: at location 0.750000 0.750000 0.750000 +DEAL:2d/3d:: at reference location 1.00000 1.00000 +DEAL:2d/3d::Particle id 7 is in cell 2.12 +DEAL:2d/3d:: at location 0.650000 0.650000 0.650000 +DEAL:2d/3d:: at reference location 0.600000 0.600000 +DEAL:2d/3d::Particle id 6 is in cell 2.12 +DEAL:2d/3d:: at location 0.550000 0.550000 0.550000 +DEAL:2d/3d:: at reference location 0.200000 0.200000 +DEAL:2d/3d::Particle id 10 is in cell 2.15 +DEAL:2d/3d:: at location 0.950000 0.950000 0.950000 +DEAL:2d/3d:: at reference location 0.800000 0.800000 +DEAL:2d/3d::Particle id 9 is in cell 2.15 +DEAL:2d/3d:: at location 0.850000 0.850000 0.850000 +DEAL:2d/3d:: at reference location 0.400000 0.400000 +DEAL:2d/3d::OK +DEAL:3d/3d::Particle number: 10 +DEAL:3d/3d::Particle id 3 is in cell 2.0 +DEAL:3d/3d:: at location 0.250000 0.250000 0.250000 +DEAL:3d/3d:: at reference location 1.00000 1.00000 1.00000 +DEAL:3d/3d::Particle id 2 is in cell 2.0 +DEAL:3d/3d:: at location 0.150000 0.150000 0.150000 +DEAL:3d/3d:: at reference location 0.600000 0.600000 0.600000 +DEAL:3d/3d::Particle id 1 is in cell 2.0 +DEAL:3d/3d:: at location 0.0500000 0.0500000 0.0500000 +DEAL:3d/3d:: at reference location 0.200000 0.200000 0.200000 +DEAL:3d/3d::Particle id 5 is in cell 2.7 +DEAL:3d/3d:: at location 0.450000 0.450000 0.450000 +DEAL:3d/3d:: at reference location 0.800000 0.800000 0.800000 +DEAL:3d/3d::Particle id 4 is in cell 2.7 +DEAL:3d/3d:: at location 0.350000 0.350000 0.350000 +DEAL:3d/3d:: at reference location 0.400000 0.400000 0.400000 +DEAL:3d/3d::Particle id 8 is in cell 2.56 +DEAL:3d/3d:: at location 0.750000 0.750000 0.750000 +DEAL:3d/3d:: at reference location 1.00000 1.00000 1.00000 +DEAL:3d/3d::Particle id 7 is in cell 2.56 +DEAL:3d/3d:: at location 0.650000 0.650000 0.650000 +DEAL:3d/3d:: at reference location 0.600000 0.600000 0.600000 +DEAL:3d/3d::Particle id 6 is in cell 2.56 +DEAL:3d/3d:: at location 0.550000 0.550000 0.550000 +DEAL:3d/3d:: at reference location 0.200000 0.200000 0.200000 +DEAL:3d/3d::Particle id 10 is in cell 2.63 +DEAL:3d/3d:: at location 0.950000 0.950000 0.950000 +DEAL:3d/3d:: at reference location 0.800000 0.800000 0.800000 +DEAL:3d/3d::Particle id 9 is in cell 2.63 +DEAL:3d/3d:: at location 0.850000 0.850000 0.850000 +DEAL:3d/3d:: at reference location 0.400000 0.400000 0.400000 +DEAL:3d/3d::OK