--- /dev/null
+Improved: Allow using ParticleHandler with serial triangulations.
+<br>
+(Daniel Arndt, 2019/11/29)
DEAL_II_NAMESPACE_OPEN
-#ifdef DEAL_II_WITH_P4EST
-
namespace Particles
{
/**
} // namespace Generators
} // namespace Particles
-#endif // DEAL_II_WITH_P4EST
-
DEAL_II_NAMESPACE_CLOSE
#endif
DEAL_II_NAMESPACE_OPEN
-#ifdef DEAL_II_WITH_P4EST
-
namespace Particles
{
/**
* This constructor is equivalent to calling the default constructor and
* the initialize function.
*/
- ParticleHandler(
- const parallel::distributed::Triangulation<dim, spacedim> &tria,
- const Mapping<dim, spacedim> & mapping,
- const unsigned int n_properties = 0);
+ ParticleHandler(const Triangulation<dim, spacedim> &tria,
+ const Mapping<dim, spacedim> & mapping,
+ const unsigned int n_properties = 0);
/**
* Destructor.
/**
* 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<dim, spacedim> &tria,
- const Mapping<dim, spacedim> &mapping,
- const unsigned int n_properties = 0);
+ initialize(const Triangulation<dim, spacedim> &tria,
+ const Mapping<dim, spacedim> & mapping,
+ const unsigned int n_properties = 0);
/**
* Clear all particle related data.
/**
* Address of the triangulation to work on.
*/
- SmartPointer<const parallel::distributed::Triangulation<dim, spacedim>,
+ SmartPointer<const Triangulation<dim, spacedim>,
ParticleHandler<dim, spacedim>>
triangulation;
*/
unsigned int handle;
-# ifdef DEAL_II_WITH_MPI
+#ifdef DEAL_II_WITH_MPI
/**
* Transfer particles that have crossed subdomain boundaries to other
* processors.
types::subdomain_id,
std::vector<
typename Triangulation<dim, spacedim>::active_cell_iterator>>());
-# endif
+#endif
/**
* Called by listener functions from Triangulation for every cell
}
} // namespace Particles
-#endif // DEAL_II_WITH_P4EST
-
DEAL_II_NAMESPACE_CLOSE
#endif
DEAL_II_NAMESPACE_OPEN
-#ifdef DEAL_II_WITH_P4EST
-
namespace Particles
{
namespace Generators
{
types::particle_index particle_index = 0;
+#ifdef DEAL_II_WITH_MPI
if (const auto tria =
dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
&triangulation))
MPI_SUM,
tria->get_communicator());
}
+#endif
for (const auto &cell : triangulation.active_cell_iterators())
{
// 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<const parallel::Triangulation<dim, spacedim> *>(
&triangulation))
MPI_SUM,
tria->get_communicator());
}
+#endif
// Calculate start id and number of local particles
start_particle_id =
} // 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
DEAL_II_NAMESPACE_OPEN
-#ifdef DEAL_II_WITH_P4EST
-
namespace Particles
{
namespace
template <int dim, int spacedim>
ParticleHandler<dim, spacedim>::ParticleHandler(
- const parallel::distributed::Triangulation<dim, spacedim> &triangulation,
- const Mapping<dim, spacedim> & mapping,
- const unsigned int n_properties)
+ const Triangulation<dim, spacedim> &triangulation,
+ const Mapping<dim, spacedim> & mapping,
+ const unsigned int n_properties)
: triangulation(&triangulation, typeid(*this).name())
, mapping(&mapping, typeid(*this).name())
, particles()
template <int dim, int spacedim>
void
ParticleHandler<dim, spacedim>::initialize(
- const parallel::distributed::Triangulation<dim, spacedim>
- & new_triangulation,
- const Mapping<dim, spacedim> &new_mapping,
- const unsigned int n_properties)
+ const Triangulation<dim, spacedim> &new_triangulation,
+ const Mapping<dim, spacedim> & new_mapping,
+ const unsigned int n_properties)
{
triangulation = &new_triangulation;
mapping = &new_mapping;
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<const parallel::Triangulation<dim, spacedim> *>(
+ &*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;
+ }
}
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<const parallel::Triangulation<dim, spacedim> *>(
+ &*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;
sorted_particles.reserve(
static_cast<vector_size>(particles_out_of_cell.size() * 1.25));
- const std::set<types::subdomain_id> ghost_owners =
- triangulation->ghost_owners();
+ std::set<types::subdomain_id> ghost_owners;
+ if (const auto parallel_triangulation =
+ dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+ &*triangulation))
+ ghost_owners = parallel_triangulation->ghost_owners();
for (const auto ghost_owner : ghost_owners)
moved_particles[ghost_owner].reserve(
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<const parallel::Triangulation<dim, spacedim> *>(
+ &*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());
ParticleHandler<dim, spacedim>::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<const parallel::Triangulation<dim, spacedim> *>(
+ &*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();
ghost_particles_by_domain;
const std::set<types::subdomain_id> 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<typename std::vector<particle_iterator>::size_type>(
}
send_recv_particles(ghost_particles_by_domain, ghost_particles);
-# endif
+#endif
}
-# ifdef DEAL_II_WITH_MPI
+#ifdef DEAL_II_WITH_MPI
template <int dim, int spacedim>
void
ParticleHandler<dim, spacedim>::send_recv_particles(
std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
&send_cells)
{
+ const auto parallel_triangulation =
+ dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+ &*triangulation);
+ Assert(
+ parallel_triangulation,
+ ExcMessage(
+ "This function is only implemented for parallel::Triangulation objects."));
+
// Determine the communication pattern
const std::set<types::subdomain_id> ghost_owners =
- triangulation->ghost_owners();
+ parallel_triangulation->ghost_owners();
const std::vector<types::subdomain_id> neighbors(ghost_owners.begin(),
ghost_owners.end());
const unsigned int n_neighbors = neighbors.size();
MPI_UNSIGNED,
neighbors[i],
mpi_tag,
- triangulation->get_communicator(),
+ parallel_triangulation->get_communicator(),
&(n_requests[2 * i]));
AssertThrowMPI(ierr);
}
MPI_UNSIGNED,
neighbors[i],
mpi_tag,
- triangulation->get_communicator(),
+ parallel_triangulation->get_communicator(),
&(n_requests[2 * i + 1]));
AssertThrowMPI(ierr);
}
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++;
}
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++;
}
"The amount of data that was read into new particles "
"does not match the amount of data sent around."));
}
-# endif
+#endif
parallel::distributed::Triangulation<dim, spacedim>
*non_const_triangulation =
const_cast<parallel::distributed::Triangulation<dim, spacedim> *>(
- &(*triangulation));
+ dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim>
+ *>(&(*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();
handle = non_const_triangulation->register_data_attach(
callback_function, /*returns_variable_size_data=*/true);
}
+#endif
}
parallel::distributed::Triangulation<dim, spacedim>
*non_const_triangulation =
const_cast<parallel::distributed::Triangulation<dim, spacedim> *>(
- &(*triangulation));
+ dynamic_cast<const parallel::distributed::Triangulation<dim, spacedim>
+ *>(&(*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
handle = numbers::invalid_unsigned_int;
update_cached_numbers();
}
+#else
+ (void)serialization;
+#endif
}
// 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.
// 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.
// 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.
}
} // 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.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()
+{
+ Triangulation<dim, spacedim> tr;
+
+ GridGenerator::hyper_cube(tr);
+ MappingQ<dim, spacedim> mapping(1);
+
+ Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
+
+
+ Point<spacedim> position;
+ position(0) = 0.3;
+ if (spacedim > 1)
+ position(1) = 0.5;
+ if (spacedim > 2)
+ position(2) = 0.7;
+
+ Point<dim> 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<dim, spacedim> particle(position, reference_position, 7);
+ deallog << "Particle location: " << particle.get_location() << std::endl;
+
+
+ std::pair<typename parallel::distributed::Triangulation<dim, spacedim>::
+ active_cell_iterator,
+ Point<dim>>
+ 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();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.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()
+{
+ {
+ Triangulation<dim, spacedim> tr;
+
+ GridGenerator::hyper_cube(tr);
+ tr.refine_global(1);
+ MappingQ<dim, spacedim> mapping(1);
+
+ Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
+
+
+ Point<spacedim> position;
+ position(0) = 0.3;
+ if (spacedim > 1)
+ position(1) = 0.5;
+ if (spacedim > 2)
+ position(2) = 0.7;
+
+ Point<dim> 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<dim, spacedim> 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<dim, spacedim> 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();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.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()
+{
+ {
+ Triangulation<dim, spacedim> tr;
+
+ GridGenerator::hyper_cube(tr);
+ tr.refine_global(1);
+ MappingQ<dim, spacedim> mapping(1);
+
+ Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
+
+ std::vector<Point<spacedim>> position(2);
+ std::vector<Point<dim>> reference_position(2);
+
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ position[0](i) = 0.25;
+ position[1](i) = 0.75;
+ }
+
+ Particles::Particle<dim, spacedim> particle1(position[0],
+ reference_position[0],
+ 0);
+ Particles::Particle<dim, spacedim> particle2(position[1],
+ reference_position[1],
+ 1);
+
+ typename Triangulation<dim, spacedim>::active_cell_iterator cell1(&tr,
+ 1,
+ 0);
+ typename Triangulation<dim, spacedim>::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<spacedim> 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();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.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()
+{
+ {
+ Triangulation<dim, spacedim> tr;
+
+ GridGenerator::hyper_cube(tr);
+ tr.refine_global(2);
+
+ MappingQ<dim, spacedim> mapping(1);
+
+ Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
+
+ std::vector<Point<spacedim>> points(10);
+ for (unsigned int i = 0; i < 10; ++i)
+ {
+ const double coordinate = static_cast<double>(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();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff