45 ProbabilityFunctionNegative,
47 <<
"Your probability density function in the particle generator "
48 "returned a negative value for the following position: "
49 << arg1 <<
". Please check your function expression.");
62 template <
int dim,
int spacedim = dim>
64 compute_local_cumulative_cell_weights(
69 std::vector<double> cumulative_cell_weights(
71 double cumulative_weight = 0.0;
87 for (
const auto &cell :
triangulation.active_cell_iterators())
89 if (cell->is_locally_owned())
91 fe_values.reinit(cell);
92 const double quadrature_point_weight =
93 probability_density_function.
value(
94 fe_values.get_quadrature_points()[0]);
97 ProbabilityFunctionNegative<dim>(
98 quadrature_formula.
point(0)));
101 cumulative_weight += quadrature_point_weight * fe_values.JxW(0);
103 cumulative_cell_weights[cell->active_cell_index()] =
107 return cumulative_cell_weights;
111 template <
int dim,
int spacedim>
115 const std::vector<
Point<dim>> & particle_reference_locations,
121 #ifdef DEAL_II_WITH_MPI
122 if (
const auto tria =
127 tria->n_locally_owned_active_cells() *
128 particle_reference_locations.size();
132 MPI_Exscan(&n_particles_to_generate,
137 tria->get_communicator());
141 for (
const auto &cell :
triangulation.active_cell_iterators())
143 if (cell->is_locally_owned())
145 for (
const auto &reference_location :
146 particle_reference_locations)
166 template <
int dim,
int spacedim>
171 std::mt19937 & random_number_generator,
176 std::uniform_real_distribution<double> uniform_distribution_01(0, 1);
183 unsigned int iteration = 0;
184 const unsigned int maximum_iterations = 100;
186 while (iteration < maximum_iterations)
188 for (
unsigned int d = 0;
d < spacedim; ++
d)
190 particle_position[
d] =
191 uniform_distribution_01(random_number_generator) *
192 (cell_bounds.second[
d] - cell_bounds.first[
d]) +
193 cell_bounds.first[
d];
212 iteration < maximum_iterations,
214 "Couldn't generate a particle position within the maximum number of tries. "
215 "The ratio between the bounding box volume in which the particle is "
216 "generated and the actual cell volume is approximately: " +
219 (cell_bounds.second - cell_bounds.first).norm_square())));
226 template <
int dim,
int spacedim>
231 const bool random_cell_selection,
235 const unsigned int random_number_seed)
237 unsigned int combined_seed = random_number_seed;
238 if (
const auto tria =
242 const unsigned int my_rank =
244 combined_seed += my_rank;
246 std::mt19937 random_number_generator(combined_seed);
251 std::vector<types::particle_index> particles_per_cell(
257 const std::vector<double> cumulative_cell_weights =
260 probability_density_function);
263 double local_weight_integral = cumulative_cell_weights.back();
264 double global_weight_integral;
266 if (
const auto tria =
270 global_weight_integral =
272 tria->get_communicator());
276 global_weight_integral = local_weight_integral;
281 "The integral of the user prescribed probability "
282 "density function over the domain equals zero, "
283 "deal.II has no way to determine the cell of "
284 "generated particles. Please ensure that the "
285 "provided function is positive in at least a "
286 "part of the domain; also check the syntax of "
291 double local_start_weight = 0.0;
293 #ifdef DEAL_II_WITH_MPI
294 if (
const auto tria =
298 MPI_Exscan(&local_weight_integral,
303 tria->get_communicator());
309 std::llround(
static_cast<double>(n_particles_to_create) *
310 local_start_weight / global_weight_integral);
314 std::llround(
static_cast<double>(n_particles_to_create) *
315 ((local_start_weight + local_weight_integral) /
316 global_weight_integral));
317 n_local_particles = end_particle_id - start_particle_id;
319 if (random_cell_selection)
324 std::uniform_real_distribution<double> uniform_distribution(
325 0.0, local_weight_integral);
329 current_particle_index < n_local_particles;
330 ++current_particle_index)
334 const double random_weight =
335 uniform_distribution(random_number_generator);
337 const std::vector<double>::const_iterator selected_cell =
339 cumulative_cell_weights.end(),
341 const unsigned int cell_index =
342 std::distance(cumulative_cell_weights.begin(), selected_cell);
344 ++particles_per_cell[cell_index];
353 for (
const auto &cell :
triangulation.active_cell_iterators())
354 if (cell->is_locally_owned())
358 static_cast<double>(n_local_particles) *
359 cumulative_cell_weights[cell->active_cell_index()] /
360 local_weight_integral);
365 particles_per_cell[cell->active_cell_index()] =
366 cumulative_particles_to_create - particles_created;
368 particles_per_cell[cell->active_cell_index()];
375 unsigned int current_particle_index = start_particle_id;
382 for (
const auto &cell :
triangulation.active_cell_iterators())
383 if (cell->is_locally_owned())
385 for (
unsigned int i = 0;
386 i < particles_per_cell[cell->active_cell_index()];
391 current_particle_index,
392 random_number_generator,
394 particles.emplace_hint(particles.end(),
396 std::move(particle));
397 ++current_particle_index;
405 template <
int dim,
int spacedim>
409 & global_bounding_boxes,
414 const auto &fe = dof_handler.
get_fe();
421 std::map<types::global_dof_index, Point<spacedim>> support_points_map;
430 std::vector<Point<spacedim>> support_points_vec;
431 support_points_vec.reserve(support_points_map.size());
432 for (
auto const &element : support_points_map)
433 support_points_vec.push_back(element.second);
436 global_bounding_boxes);
440 template <
int dim,
int spacedim>
447 & global_bounding_boxes,
451 const std::vector<Point<dim>> &particle_reference_locations =
453 std::vector<Point<spacedim>> points_to_generate;
456 for (
const auto &cell :
triangulation.active_cell_iterators())
458 if (cell->is_locally_owned())
460 for (
const auto &reference_location :
461 particle_reference_locations)
466 points_to_generate.push_back(position_real);
471 global_bounding_boxes);
476 #include "generators.inst"