#ifndef dealii_particles_particle_generator_h
#define dealii_particles_particle_generator_h
+#include <deal.II/base/function.h>
+
#include <deal.II/distributed/tria.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/particles/particle.h>
#include <deal.II/particles/particle_handler.h>
+#include <random>
+
DEAL_II_NAMESPACE_OPEN
#ifdef DEAL_II_WITH_P4EST
const Mapping<dim, spacedim> & mapping =
StaticMappingQ1<dim, spacedim>::mapping);
+ /**
+ * A function that generates one particle at a random location in cell @p cell and with
+ * index @p id. The function expects a random number generator to avoid the expensive generation
+ * and destruction of a generator for every particle and optionally takes
+ * into account a mapping for the cell. The algorithm implemented in the
+ * function is described in @cite GLHPW2018. In short, the algorithm
+ * generates
+ * random locations within the bounding box of the @p cell. It then inverts the mapping
+ * to check if the generated particle is within the cell itself. This makes
+ * sure the algorithm produces statistically random locations even for
+ * nonlinear mappings and distorted cells. However, if the ratio between
+ * bounding box and cell volume becomes very large -- i.e. the cells become
+ * strongly deformed, for example a pencil shaped cell that lies diagonally
+ * in the domain -- then the algorithm can become very inefficient.
+ * Therefore, it only tries to find a location ni the cell a fixed number of
+ * times before throwing an error message.
+ *
+ * @param[in] cell The cell in which a particle is generated.
+ *
+ * @param[in] id The particle index that will be assigned to the new
+ * particle.
+ *
+ * @param[in,out] random_number_generator A random number generator that
+ * will be used for the creation of th particle.
+ *
+ * @param[in] mapping An optional mapping object that is used to map
+ * reference location in the unit cell to the real cell. If no mapping is
+ * provided a MappingQ1 is assumed.
+ */
+ template <int dim, int spacedim = dim>
+ Particle<dim, spacedim>
+ random_particle_in_cell(
+ const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+ const types::particle_index id,
+ std::mt19937 & random_number_generator,
+ const Mapping<dim, spacedim> &mapping =
+ StaticMappingQ1<dim, spacedim>::mapping);
+
+ /**
+ * A function that generates particles randomly in the domain with a
+ * particle density
+ * according to a provided probability density function @p probability_density_function.
+ * The total number of particles that is added to the @p particle_handler object is
+ * @p n_particles_to_create. An optional @p mapping argument
+ * can be used to map from @p particle_reference_locations to the real particle locations.
+ * The function can compute the number of particles per cell either
+ * deterministically by computing the integral of the probability density
+ * function for each cell and creating
+ * particles accordingly (if option @p random_cell_selection set to false), or it can
+ * select cells randomly based on the probability density function and the
+ * cell size
+ * (if option @p random_cell_selection set to true). In either case the position of
+ * individual particles inside the cell is computed randomly.
+ *
+ * The algorithm implemented in the function is described in @cite
+ * GLHPW2018.
+ *
+ * @param[in] triangulation The triangulation associated with the @p particle_handler.
+ *
+ * @param[in] probability_density_function A function with non-negative
+ * values that determines the probability density of a particle to be
+ * generated in this location. The function does not need to be normalized.
+ *
+ * @param[in] random_cell_selection A bool that determines, how the number
+ * of particles per cell is computed (see the description above).
+ *
+ * @param[in] n_particles_to_create The number of particles that will be
+ * created by this function.
+ *
+ * @param[in,out] particle_handler The particle handler that will take
+ * ownership of the generated particles.
+ *
+ * @param[in] mapping An optional mapping object that is used to map
+ * reference location in the unit cell to the real cells of the
+ * triangulation. If no mapping is provided a MappingQ1 is assumed.
+ *
+ * @param[in] random_number_seed An optional seed that determines the
+ * initial state of the random number generator. Use the same number to get
+ * a reproducible particle distribution, or a changing number (e.g. based on
+ * system time) to generate different particle distributions for each call
+ * to this function.
+ */
+ template <int dim, int spacedim = dim>
+ void
+ probabilistic_locations(
+ const Triangulation<dim, spacedim> &triangulation,
+ const Function<spacedim> & probability_density_function,
+ const bool random_cell_selection,
+ const types::particle_index n_particles_to_create,
+ ParticleHandler<dim, spacedim> & particle_handler,
+ const Mapping<dim, spacedim> & mapping =
+ StaticMappingQ1<dim, spacedim>::mapping,
+ const unsigned int random_number_seed = 5432);
+
} // namespace Generators
} // namespace Particles
//
// ---------------------------------------------------------------------
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/signaling_nan.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+
#include <deal.II/particles/generators.h>
DEAL_II_NAMESPACE_OPEN
{
namespace Generators
{
+ namespace
+ {
+ /**
+ * Exception
+ */
+ template <int dim>
+ DeclException1(
+ ProbabilityFunctionNegative,
+ Point<dim>,
+ << "Your probability density function in the particle generator "
+ "returned a negative value for the following position: "
+ << arg1 << ". Please check your function expression.");
+
+
+ // This function integrates the given probability density
+ // function over all cells of the triangulation. For each
+ // cell it stores the cumulative sum (of all previous
+ // cells including the current cell) in a vector that is then
+ // returned. Therefore the returned vector has as many entries
+ // as active cells, the first entry being the integral over the
+ // first cell, and the last entry the integral over the whole
+ // locally owned domain. Cells that are not locally owned
+ // simply store the same value as the cell before (equivalent
+ // to assuming a probability density function value of 0).
+ template <int dim, int spacedim = dim>
+ std::vector<double>
+ compute_local_cumulative_cell_weights(
+ const Triangulation<dim, spacedim> &triangulation,
+ const Mapping<dim, spacedim> & mapping,
+ const Function<spacedim> & probability_density_function)
+ {
+ std::vector<double> cumulative_cell_weights(
+ triangulation.n_active_cells());
+ double cumulative_weight = 0.0;
+
+ // Evaluate function at all cell midpoints
+ const QMidpoint<dim> quadrature_formula;
+
+ // In the simplest case we do not even need a FEValues object, because
+ // using cell->center() and cell->measure() would be equivalent. This
+ // fails however for higher-order mappings.
+ FE_Nothing<dim, spacedim> alibi_finite_element;
+ FEValues<dim, spacedim> fe_values(mapping,
+ alibi_finite_element,
+ quadrature_formula,
+ update_quadrature_points |
+ update_JxW_values);
+
+ // compute the integral weight
+ for (const auto &cell : triangulation.active_cell_iterators())
+ {
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit(cell);
+ const double quadrature_point_weight =
+ probability_density_function.value(
+ fe_values.get_quadrature_points()[0]);
+
+ AssertThrow(quadrature_point_weight >= 0.0,
+ ProbabilityFunctionNegative<dim>(
+ quadrature_formula.point(0)));
+
+ // get_cell_weight makes sure to return positive values
+ cumulative_weight += quadrature_point_weight * fe_values.JxW(0);
+ }
+ cumulative_cell_weights[cell->active_cell_index()] =
+ cumulative_weight;
+ }
+
+ return cumulative_cell_weights;
+ }
+ } // namespace
+
template <int dim, int spacedim>
void
regular_reference_locations(
{
types::particle_index particle_index = 0;
- if (const auto tria = dynamic_cast<
- const parallel::distributed::Triangulation<dim, spacedim> *>(
- &triangulation))
+ if (const auto tria =
+ dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+ &triangulation))
{
const types::particle_index n_particles_to_generate =
tria->n_locally_owned_active_cells() *
particle_handler.update_cached_numbers();
}
+
+
+
+ template <int dim, int spacedim = dim>
+ Particle<dim, spacedim>
+ random_particle_in_cell(
+ const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+ const types::particle_index id,
+ std::mt19937 & random_number_generator,
+ const Mapping<dim, spacedim> &mapping)
+ {
+ // Uniform distribution on the interval [0,1]. This
+ // will be used to generate random particle locations.
+ std::uniform_real_distribution<double> uniform_distribution_01(0, 1);
+
+ const BoundingBox<spacedim> cell_bounding_box(cell->bounding_box());
+ const std::pair<Point<spacedim>, Point<spacedim>> cell_bounds(
+ cell_bounding_box.get_boundary_points());
+
+ // Generate random points in these bounds until one is within the cell
+ unsigned int iteration = 0;
+ const unsigned int maximum_iterations = 100;
+ Point<spacedim> particle_position;
+ while (iteration < maximum_iterations)
+ {
+ for (unsigned int d = 0; d < spacedim; ++d)
+ {
+ particle_position[d] =
+ uniform_distribution_01(random_number_generator) *
+ (cell_bounds.second[d] - cell_bounds.first[d]) +
+ cell_bounds.first[d];
+ }
+ try
+ {
+ const Point<dim> p_unit =
+ mapping.transform_real_to_unit_cell(cell, particle_position);
+ if (GeometryInfo<dim>::is_inside_unit_cell(p_unit))
+ {
+ // Generate the particle
+ return Particle<dim, spacedim>(particle_position, p_unit, id);
+ }
+ }
+ catch (typename Mapping<dim>::ExcTransformationFailed &)
+ {
+ // The point is not in this cell. Do nothing, just try again.
+ }
+ ++iteration;
+ }
+ AssertThrow(
+ iteration < maximum_iterations,
+ ExcMessage(
+ "Couldn't generate a particle position within the maximum number of tries. "
+ "The ratio between the bounding box volume in which the particle is "
+ "generated and the actual cell volume is approximately: " +
+ std::to_string(
+ cell->measure() /
+ (cell_bounds.second - cell_bounds.first).norm_square())));
+
+ return Particle<dim, spacedim>();
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ probabilistic_locations(
+ const Triangulation<dim, spacedim> &triangulation,
+ const Function<spacedim> & probability_density_function,
+ const bool random_cell_selection,
+ const types::particle_index n_particles_to_create,
+ ParticleHandler<dim, spacedim> & particle_handler,
+ const Mapping<dim, spacedim> & mapping,
+ const unsigned int random_number_seed)
+ {
+ unsigned int combined_seed = random_number_seed;
+ if (const auto tria =
+ dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+ &triangulation))
+ {
+ const unsigned int my_rank =
+ Utilities::MPI::this_mpi_process(tria->get_communicator());
+ combined_seed += my_rank;
+ }
+ std::mt19937 random_number_generator(combined_seed);
+
+ types::particle_index start_particle_id(0);
+ types::particle_index n_local_particles(0);
+
+ std::vector<types::particle_index> particles_per_cell(
+ triangulation.n_active_cells(), 0);
+
+ // First determine how many particles to generate in which cell
+ {
+ // Get the local accumulated probabilities for every cell
+ const std::vector<double> cumulative_cell_weights =
+ compute_local_cumulative_cell_weights(triangulation,
+ mapping,
+ probability_density_function);
+
+ // Sum the local integrals over all nodes
+ double local_weight_integral = cumulative_cell_weights.back();
+ double global_weight_integral;
+
+ if (const auto tria =
+ dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+ &triangulation))
+ {
+ global_weight_integral =
+ Utilities::MPI::sum(local_weight_integral,
+ tria->get_communicator());
+ }
+ else
+ {
+ global_weight_integral = local_weight_integral;
+ }
+
+ AssertThrow(global_weight_integral > std::numeric_limits<double>::min(),
+ ExcMessage(
+ "The integral of the user prescribed probability "
+ "density function over the domain equals zero, "
+ "deal.II has no way to determine the cell of "
+ "generated particles. Please ensure that the "
+ "provided function is positive in at least a "
+ "part of the domain; also check the syntax of "
+ "the function."));
+
+ // Determine the starting weight of this process, which is the sum of
+ // the weights of all processes with a lower rank
+ double local_start_weight = 0.0;
+
+ if (const auto tria =
+ dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+ &triangulation))
+ {
+ MPI_Exscan(&local_weight_integral,
+ &local_start_weight,
+ 1,
+ MPI_DOUBLE,
+ MPI_SUM,
+ tria->get_communicator());
+ }
+
+ // Calculate start id and number of local particles
+ start_particle_id =
+ std::llround(static_cast<double>(n_particles_to_create) *
+ local_start_weight / global_weight_integral);
+ n_local_particles =
+ std::llround(static_cast<double>(n_particles_to_create) *
+ local_weight_integral / global_weight_integral);
+
+ if (random_cell_selection)
+ {
+ // Uniform distribution on the interval [0,local_weight_integral).
+ // This will be used to randomly select cells for all local
+ // particles.
+ std::uniform_real_distribution<double> uniform_distribution(
+ 0.0, local_weight_integral);
+
+ // Loop over all particles to create locally and pick their cells
+ for (types::particle_index current_particle_index = 0;
+ current_particle_index < n_local_particles;
+ ++current_particle_index)
+ {
+ // Draw the random number that determines the cell of the
+ // particle
+ const double random_weight =
+ uniform_distribution(random_number_generator);
+
+ const std::vector<double>::const_iterator selected_cell =
+ std::lower_bound(cumulative_cell_weights.begin(),
+ cumulative_cell_weights.end(),
+ random_weight);
+ const unsigned int cell_index =
+ std::distance(cumulative_cell_weights.begin(), selected_cell);
+
+ ++particles_per_cell[cell_index];
+ }
+ }
+ else
+ {
+ // Compute number of particles per cell according to the ratio
+ // between their weight and the local weight integral
+ types::particle_index particles_created = 0;
+
+ for (const auto &cell : triangulation.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ const types::particle_index cumulative_particles_to_create =
+ std::llround(
+ static_cast<double>(n_local_particles) *
+ cumulative_cell_weights[cell->active_cell_index()] /
+ local_weight_integral);
+
+ // Compute particles for this cell as difference between
+ // number of particles that should be created including this
+ // cell minus the number of particles already created.
+ particles_per_cell[cell->active_cell_index()] =
+ cumulative_particles_to_create - particles_created;
+ particles_created +=
+ particles_per_cell[cell->active_cell_index()];
+ }
+ }
+ }
+
+ // Now generate as many particles per cell as determined above
+ {
+ unsigned int current_particle_index = start_particle_id;
+
+ std::multimap<
+ typename Triangulation<dim, spacedim>::active_cell_iterator,
+ Particle<dim, spacedim>>
+ particles;
+
+ for (const auto &cell : triangulation.active_cell_iterators())
+ if (cell->is_locally_owned())
+ {
+ for (unsigned int i = 0;
+ i < particles_per_cell[cell->active_cell_index()];
+ ++i)
+ {
+ Particle<dim, spacedim> particle =
+ random_particle_in_cell(cell,
+ current_particle_index,
+ random_number_generator,
+ mapping);
+ particles.emplace_hint(particles.end(),
+ cell,
+ std::move(particle));
+ ++current_particle_index;
+ }
+ }
+
+ particle_handler.insert_particles(particles);
+ }
+ }
} // namespace Generators
} // namespace Particles
projects / dealii.git / commitdiff