*/
void clear_particles();
+ /**
+ * Updates all internally cached numbers. Note that all functions that
+ * modify internal data structures and act on multiple particles will
+ * call this function automatically (e.g. insert_particles), while
+ * functions that act on single particles will not call this function
+ * (e.g. insert_particle). This is done because the update is
+ * expensive compared to single operations.
+ */
+ void update_cache();
+
/**
* Return an iterator to the first particle.
*/
*/
unsigned int data_offset;
- /**
- * Calculates the number of particles in the global model domain.
- */
- void
- update_n_global_particles();
-
- /**
- * Calculates and stores the number of particles in the cell that
- * contains the most particles in the global model (stored in the
- * member variable global_max_particles_per_cell). This variable is a
- * state variable, because it is needed to serialize and deserialize
- * the particle data correctly in parallel (it determines the size of
- * the data chunks per cell that are stored and read). Before accessing
- * the variable this function has to be called, unless the state was
- * read from another source (e.g. after resuming from a checkpoint).
- */
- void
- update_global_max_particles_per_cell();
-
- /**
- * Calculates the next free particle index in the global model domain.
- * This equals one plus the highest particle index currently active.
- */
- void
- update_next_free_particle_index();
-
/**
* Transfer particles that have crossed subdomain boundaries to other
* processors.
+ template <int dim,int spacedim>
+ void
+ ParticleHandler<dim,spacedim>::update_cache()
+ {
+
+ types::particle_index locally_highest_index = 0;
+ unsigned int local_max_particles_per_cell = 0;
+ unsigned int current_particles_per_cell = 0;
+ typename Triangulation<dim,spacedim>::active_cell_iterator current_cell;
+
+ for (particle_iterator particle = begin(); particle != end(); ++particle)
+ {
+ locally_highest_index = std::max(locally_highest_index,particle->get_id());
+
+ if (particle->get_surrounding_cell(*triangulation) != current_cell)
+ current_particles_per_cell = 0;
+
+ ++current_particles_per_cell;
+ local_max_particles_per_cell = 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());
+ }
+
+
+
+
template <int dim,int spacedim>
typename ParticleHandler<dim,spacedim>::particle_iterator
ParticleHandler<dim,spacedim>::begin() const
ParticleHandler<dim,spacedim>::insert_particles(const std::multimap<types::LevelInd, Particle<dim,spacedim> > &new_particles)
{
particles.insert(new_particles.begin(),new_particles.end());
+ update_cache();
}
- template <int dim,int spacedim>
- void
- ParticleHandler<dim,spacedim>::update_n_global_particles()
- {
- global_number_of_particles = dealii::Utilities::MPI::sum (particles.size(), triangulation->get_communicator());
- }
-
-
-
template <int dim,int spacedim>
unsigned int
ParticleHandler<dim,spacedim>::n_particles_in_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell) const
- template <int dim,int spacedim>
- void
- ParticleHandler<dim,spacedim>::update_next_free_particle_index()
- {
- types::particle_index locally_highest_index = 0;
- for (particle_iterator particle = begin(); particle != end(); ++particle)
- locally_highest_index = std::max(locally_highest_index,particle->get_id());
-
- next_free_particle_index = dealii::Utilities::MPI::max (locally_highest_index, triangulation->get_communicator()) + 1;
- }
-
-
-
- template <int dim,int spacedim>
- void
- ParticleHandler<dim,spacedim>::update_global_max_particles_per_cell()
- {
- unsigned int local_max_particles_per_cell(0);
- typename Triangulation<dim,spacedim>::active_cell_iterator cell = triangulation->begin_active();
- for (; cell!=triangulation->end(); ++cell)
- if (cell->is_locally_owned())
- {
- local_max_particles_per_cell = std::max(local_max_particles_per_cell,
- n_particles_in_cell(cell));
- }
-
- global_max_particles_per_cell = dealii::Utilities::MPI::max(local_max_particles_per_cell,triangulation->get_communicator());
- }
-
-
-
template <int dim, int spacedim>
PropertyPool &
ParticleHandler<dim,spacedim>::get_property_pool () const
remove_particle(particles_out_of_cell[i]);
particles.insert(sorted_particles_map.begin(),sorted_particles_map.end());
+ update_cache();
}
// 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_global_max_particles_per_cell();
+ update_cache();
if (global_max_particles_per_cell > 0)
{
// Reset offset and update global number of particles. The number
// can change because of discarded or newly generated particles
data_offset = numbers::invalid_unsigned_int;
- update_n_global_particles();
+ update_cache();
}
}
reference_position(0) = 0.2;
reference_position(1) = 0.4;
- const Particles::types::particle_index index(7);
+ const types::particle_index index(7);
Particles::Particle<2> particle(position,reference_position,index);
reference_position(0) = 0.2;
reference_position(1) = 0.4;
- const Particles::types::particle_index index(7);
+ const types::particle_index index(7);
std::vector<double> properties = {0.15,0.45,0.75};
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the creation and destruction of particle within the particle handler class.
+
+#include "../tests.h"
+#include <deal.II/particles/particle_handler.h>
+
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/fe/mapping_q.h>
+
+template <int dim, int spacedim>
+void test ()
+{
+ {
+ parallel::distributed::Triangulation<dim,spacedim> tr(MPI_COMM_WORLD);
+
+ GridGenerator::hyper_cube(tr);
+ MappingQ<dim,spacedim> mapping(1);
+
+ Particles::ParticleHandler<dim,spacedim> particle_handler(tr,mapping);
+
+ Particles::Particle<dim,spacedim> particle;
+
+ 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(1) = 0.6;
+
+ particle.set_location(position);
+ particle.set_reference_location(reference_position);
+ 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_cache();
+
+ deallog << "Particle number: " << particle_handler.n_global_particles() << std::endl;
+
+ for (auto particle = particle_handler.begin(); particle != particle_handler.end(); ++particle)
+ {
+ 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[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+
+ initlog();
+
+ test<2,2>();
+ test<2,3>();
+
+ test<3,3>();
+}
--- /dev/null
+
+DEAL::Particle location: 0.300000 0.500000
+DEAL::Particle number: 1
+DEAL::Particle location: 0.300000 0.500000
+DEAL::Particle reference location: 0.200000 0.400000
+DEAL::OK
+DEAL::Particle location: 0.300000 0.500000 0.700000
+DEAL::Particle number: 1
+DEAL::Particle location: 0.300000 0.500000 0.700000
+DEAL::Particle reference location: 0.200000 0.400000
+DEAL::OK
+DEAL::Particle location: 0.300000 0.500000 0.700000
+DEAL::Particle number: 1
+DEAL::Particle location: 0.300000 0.500000 0.700000
+DEAL::Particle reference location: 0.200000 0.600000 0.00000
+DEAL::OK
if (dim>1)
reference_position(1) = 0.4;
- const Particles::types::particle_index index(7);
+ const types::particle_index index(7);
std::vector<double> properties = {0.15,0.45,0.75};
particle.set_property_pool(pool);
particle.set_properties(ArrayView<double>(&properties[0],properties.size()));
- std::multimap<Particles::types::LevelInd,Particles::Particle<dim> > map;
+ std::multimap<types::LevelInd,Particles::Particle<dim> > map;
- Particles::types::LevelInd level_index = std::make_pair(0,0);
+ types::LevelInd level_index = std::make_pair(0,0);
map.insert(std::make_pair(level_index,particle));
particle.get_properties()[0] = 0.05;
projects / dealii.git / commitdiff