46 ProbabilityFunctionNegative,
48 <<
"Your probability density function in the particle generator "
49 "returned a negative value for the following position: "
50 << arg1 <<
". Please check your function expression.");
63 template <
int dim,
int spacedim = dim>
65 compute_local_cumulative_cell_weights(
70 std::vector<double> cumulative_cell_weights(
72 double cumulative_weight = 0.0;
88 for (
const auto &cell :
triangulation.active_cell_iterators())
90 if (cell->is_locally_owned())
92 fe_values.reinit(cell);
93 const double quadrature_point_weight =
94 probability_density_function.
value(
95 fe_values.get_quadrature_points()[0]);
98 ProbabilityFunctionNegative<dim>(
99 quadrature_formula.point(0)));
102 cumulative_weight += quadrature_point_weight * fe_values.JxW(0);
104 cumulative_cell_weights[cell->active_cell_index()] =
108 return cumulative_cell_weights;
118 template <
int dim,
int spacedim>
120 random_location_in_cell(
123 std::mt19937 &random_number_generator)
125 Assert(cell->reference_cell().is_hyper_cube() ==
true,
130 std::uniform_real_distribution<double> uniform_distribution_01(0, 1);
134 cell_bounding_box.get_boundary_points());
138 const unsigned int n_attempts = 100;
139 for (
unsigned int i = 0; i < n_attempts; ++i)
142 for (
unsigned int d = 0;
d < spacedim; ++
d)
144 position[
d] = uniform_distribution_01(random_number_generator) *
145 (cell_bounds.second[
d] - cell_bounds.first[
d]) +
146 cell_bounds.first[d];
154 if (cell->reference_cell().contains_point(reference_position))
155 return {position, reference_position};
169 for (
unsigned int d = 0;
d < dim; ++
d)
170 reference_position[d] =
171 uniform_distribution_01(random_number_generator);
176 return {position, reference_position};
182 template <
int dim,
int spacedim>
186 const std::vector<
Point<dim>> &particle_reference_locations,
192 triangulation.n_active_cells() * particle_reference_locations.size();
194#ifdef DEAL_II_WITH_MPI
195 if (
const auto tria =
199 n_particles_to_generate = tria->n_locally_owned_active_cells() *
200 particle_reference_locations.size();
204 const int ierr = MPI_Exscan(
205 &n_particles_to_generate,
208 Utilities::MPI::mpi_type_id_for_type<types::particle_index>,
210 tria->get_mpi_communicator());
216 n_particles_to_generate);
218 for (
const auto &cell :
triangulation.active_cell_iterators() |
221 for (
const auto &reference_location : particle_reference_locations)
239 template <
int dim,
int spacedim>
244 std::mt19937 &random_number_generator,
247 const auto position_and_reference_position =
248 random_location_in_cell(cell, mapping, random_number_generator);
250 position_and_reference_position.second,
256 template <
int dim,
int spacedim>
261 std::mt19937 &random_number_generator,
265 const auto position_and_reference_position =
266 random_location_in_cell(cell, mapping, random_number_generator);
268 position_and_reference_position.first,
269 position_and_reference_position.second,
276 template <
int dim,
int spacedim>
281 const bool random_cell_selection,
285 const unsigned int random_number_seed)
287 unsigned int combined_seed = random_number_seed;
288 if (
const auto tria =
296 std::mt19937 random_number_generator(combined_seed);
301 std::vector<types::particle_index> particles_per_cell(
307 const std::vector<double> cumulative_cell_weights =
310 probability_density_function);
313 double local_weight_integral = (cumulative_cell_weights.size() > 0) ?
314 cumulative_cell_weights.back() :
318 double local_start_weight = numbers::signaling_nan<double>();
319 double global_weight_integral = numbers::signaling_nan<double>();
321 if (
const auto tria =
325 std::tie(local_start_weight, global_weight_integral) =
327 local_weight_integral, tria->get_mpi_communicator());
331 local_start_weight = 0;
332 global_weight_integral = local_weight_integral;
335 AssertThrow(global_weight_integral > std::numeric_limits<double>::min(),
337 "The integral of the user prescribed probability "
338 "density function over the domain equals zero, "
339 "deal.II has no way to determine the cell of "
340 "generated particles. Please ensure that the "
341 "provided function is positive in at least a "
342 "part of the domain; also check the syntax of "
347 std::llround(
static_cast<double>(n_particles_to_create) *
348 local_start_weight / global_weight_integral);
352 std::llround(
static_cast<double>(n_particles_to_create) *
353 ((local_start_weight + local_weight_integral) /
354 global_weight_integral));
355 n_local_particles = end_particle_id - start_particle_id;
357 if (random_cell_selection)
362 std::uniform_real_distribution<double> uniform_distribution(
363 0.0, local_weight_integral);
367 current_particle_index < n_local_particles;
368 ++current_particle_index)
372 const double random_weight =
373 uniform_distribution(random_number_generator);
375 const std::vector<double>::const_iterator selected_cell =
376 std::lower_bound(cumulative_cell_weights.begin(),
377 cumulative_cell_weights.end(),
380 std::distance(cumulative_cell_weights.begin(), selected_cell);
387 if (local_weight_integral > 0)
393 for (
const auto &cell :
triangulation.active_cell_iterators() |
398 static_cast<double>(n_local_particles) *
399 cumulative_cell_weights[cell->active_cell_index()] /
400 local_weight_integral);
405 particles_per_cell[cell->active_cell_index()] =
406 cumulative_particles_to_create - particles_created;
408 particles_per_cell[cell->active_cell_index()];
418 unsigned int current_particle_index = start_particle_id;
420 for (
const auto &cell :
triangulation.active_cell_iterators() |
423 for (
unsigned int i = 0;
424 i < particles_per_cell[cell->active_cell_index()];
428 current_particle_index,
429 random_number_generator,
433 ++current_particle_index;
443 template <
int dim,
int spacedim>
447 &global_bounding_boxes,
451 const std::vector<std::vector<double>> &properties)
453 const auto &fe = dof_handler.
get_fe();
460 const std::map<types::global_dof_index, Point<spacedim>>
466 std::vector<Point<spacedim>> support_points_vec;
467 support_points_vec.reserve(support_points_map.size());
468 for (
const auto &element : support_points_map)
469 support_points_vec.push_back(element.second);
472 global_bounding_boxes,
477 template <
int dim,
int spacedim>
484 &global_bounding_boxes,
487 const std::vector<std::vector<double>> &properties)
489 const std::vector<Point<dim>> &particle_reference_locations =
491 std::vector<Point<spacedim>> points_to_generate;
492 points_to_generate.reserve(particle_reference_locations.size() *
496 for (
const auto &cell :
triangulation.active_cell_iterators() |
499 for (
const auto &reference_location : particle_reference_locations)
503 points_to_generate.push_back(position_real);
507 global_bounding_boxes,
513#include "particles/generators.inst"
unsigned int size() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Abstract base class for mapping classes.
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
void update_cached_numbers()
void reserve(const std::size_t n_particles)
std::map< unsigned int, IndexSet > insert_global_particles(const std::vector< Point< spacedim > > &positions, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, const std::vector< std::vector< double > > &properties={}, const std::vector< types::particle_index > &ids={})
types::particle_index n_locally_owned_particles() const
particle_iterator insert_particle(const Particle< dim, spacedim > &particle, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
const std::vector< Point< dim > > & get_points() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertThrowMPI(error_code)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
const unsigned int my_rank
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=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const unsigned int random_number_seed=5432)
ParticleIterator< dim, spacedim > random_particle_in_cell_insert(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell, const types::particle_index id, std::mt19937 &random_number_generator, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
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=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
void dof_support_points(const DoFHandler< dim, spacedim > &dof_handler, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const ComponentMask &components={}, const std::vector< std::vector< double > > &properties={})
void regular_reference_locations(const Triangulation< dim, spacedim > &triangulation, const std::vector< Point< dim > > &particle_reference_locations, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::pair< T, T > partial_and_total_sum(const T &value, const MPI_Comm comm)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation