Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
generators.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2019 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
19
21
24
26
28
29#include <limits>
30
32
33namespace Particles
34{
35 namespace Generators
36 {
37 namespace
38 {
42 template <int dim>
44 ProbabilityFunctionNegative,
46 << "Your probability density function in the particle generator "
47 "returned a negative value for the following position: "
48 << arg1 << ". Please check your function expression.");
49
50
51 // This function integrates the given probability density
52 // function over all cells of the triangulation. For each
53 // cell it stores the cumulative sum (of all previous
54 // cells including the current cell) in a vector that is then
55 // returned. Therefore the returned vector has as many entries
56 // as active cells, the first entry being the integral over the
57 // first cell, and the last entry the integral over the whole
58 // locally owned domain. Cells that are not locally owned
59 // simply store the same value as the cell before (equivalent
60 // to assuming a probability density function value of 0).
61 template <int dim, int spacedim = dim>
62 std::vector<double>
63 compute_local_cumulative_cell_weights(
65 const Mapping<dim, spacedim> &mapping,
66 const Function<spacedim> &probability_density_function)
67 {
68 std::vector<double> cumulative_cell_weights(
69 triangulation.n_active_cells());
70 double cumulative_weight = 0.0;
71
72 // Evaluate function at all cell midpoints
73 const QMidpoint<dim> quadrature_formula;
74
75 // In the simplest case we do not even need a FEValues object, because
76 // using cell->center() and cell->measure() would be equivalent. This
77 // fails however for higher-order mappings.
78 FE_Nothing<dim, spacedim> alibi_finite_element;
79 FEValues<dim, spacedim> fe_values(mapping,
80 alibi_finite_element,
81 quadrature_formula,
84
85 // compute the integral weight
86 for (const auto &cell : triangulation.active_cell_iterators())
87 {
88 if (cell->is_locally_owned())
89 {
90 fe_values.reinit(cell);
91 const double quadrature_point_weight =
92 probability_density_function.value(
93 fe_values.get_quadrature_points()[0]);
94
95 AssertThrow(quadrature_point_weight >= 0.0,
96 ProbabilityFunctionNegative<dim>(
97 quadrature_formula.point(0)));
98
99 // get_cell_weight makes sure to return positive values
100 cumulative_weight += quadrature_point_weight * fe_values.JxW(0);
101 }
102 cumulative_cell_weights[cell->active_cell_index()] =
103 cumulative_weight;
104 }
105
106 return cumulative_cell_weights;
107 }
108
109
110
111 // This function generates a random position in the given cell and
112 // returns the position and its coordinates in the reference cell. It
113 // first tries to generate a random and uniformly distributed point in
114 // real space, but if that fails (e.g. because the cell has a bad aspect
115 // ratio) it reverts to generating a random point in the unit cell.
116 template <int dim, int spacedim>
117 std::pair<Point<spacedim>, Point<dim>>
118 random_location_in_cell(
120 const Mapping<dim, spacedim> &mapping,
121 std::mt19937 &random_number_generator)
122 {
123 Assert(cell->reference_cell().is_hyper_cube() == true,
125
126 // Uniform distribution on the interval [0,1]. This
127 // will be used to generate random particle locations.
128 std::uniform_real_distribution<double> uniform_distribution_01(0, 1);
129
130 const BoundingBox<spacedim> cell_bounding_box(cell->bounding_box());
131 const std::pair<Point<spacedim>, Point<spacedim>> &cell_bounds(
132 cell_bounding_box.get_boundary_points());
133
134 // Generate random points in these bounds until one is within the cell
135 // or we exceed the maximum number of attempts.
136 const unsigned int n_attempts = 100;
137 for (unsigned int i = 0; i < n_attempts; ++i)
138 {
139 Point<spacedim> position;
140 for (unsigned int d = 0; d < spacedim; ++d)
141 {
142 position[d] = uniform_distribution_01(random_number_generator) *
143 (cell_bounds.second[d] - cell_bounds.first[d]) +
144 cell_bounds.first[d];
145 }
146
147 try
148 {
149 const Point<dim> reference_position =
150 mapping.transform_real_to_unit_cell(cell, position);
151
152 if (cell->reference_cell().contains_point(reference_position))
153 return {position, reference_position};
154 }
156 {
157 // The point is not in this cell. Do nothing, just try again.
158 }
159 }
160
161 // If the above algorithm has not worked (e.g. because of badly
162 // deformed cells), retry generating particles
163 // randomly within the reference cell and then mapping it to to
164 // real space. This is not generating a
165 // uniform distribution in real space, but will always succeed.
166 Point<dim> reference_position;
167 for (unsigned int d = 0; d < dim; ++d)
168 reference_position[d] =
169 uniform_distribution_01(random_number_generator);
170
171 const Point<spacedim> position =
172 mapping.transform_unit_to_real_cell(cell, reference_position);
173
174 return {position, reference_position};
175 }
176 } // namespace
177
178
179
180 template <int dim, int spacedim>
181 void
184 const std::vector<Point<dim>> &particle_reference_locations,
185 ParticleHandler<dim, spacedim> &particle_handler,
186 const Mapping<dim, spacedim> &mapping)
187 {
188 types::particle_index particle_index = 0;
189 types::particle_index n_particles_to_generate =
190 triangulation.n_active_cells() * particle_reference_locations.size();
191
192#ifdef DEAL_II_WITH_MPI
193 if (const auto tria =
194 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
196 {
197 n_particles_to_generate = tria->n_locally_owned_active_cells() *
198 particle_reference_locations.size();
199
200 // The local particle start index is the number of all particles
201 // generated on lower MPI ranks.
202 const int ierr = MPI_Exscan(
203 &n_particles_to_generate,
204 &particle_index,
205 1,
206 Utilities::MPI::mpi_type_id_for_type<types::particle_index>,
207 MPI_SUM,
208 tria->get_mpi_communicator());
209 AssertThrowMPI(ierr);
210 }
211#endif
212
213 particle_handler.reserve(particle_handler.n_locally_owned_particles() +
214 n_particles_to_generate);
215
216 for (const auto &cell : triangulation.active_cell_iterators() |
218 {
219 for (const auto &reference_location : particle_reference_locations)
220 {
221 const Point<spacedim> position_real =
222 mapping.transform_unit_to_real_cell(cell, reference_location);
223
224 particle_handler.insert_particle(position_real,
225 reference_location,
226 particle_index,
227 cell);
228 ++particle_index;
229 }
230 }
231
232 particle_handler.update_cached_numbers();
233 }
234
235
236
237 template <int dim, int spacedim>
241 const types::particle_index id,
242 std::mt19937 &random_number_generator,
243 const Mapping<dim, spacedim> &mapping)
244 {
245 const auto position_and_reference_position =
246 random_location_in_cell(cell, mapping, random_number_generator);
247 return Particle<dim, spacedim>(position_and_reference_position.first,
248 position_and_reference_position.second,
249 id);
250 }
251
252
253
254 template <int dim, int spacedim>
258 const types::particle_index id,
259 std::mt19937 &random_number_generator,
260 ParticleHandler<dim, spacedim> &particle_handler,
261 const Mapping<dim, spacedim> &mapping)
262 {
263 const auto position_and_reference_position =
264 random_location_in_cell(cell, mapping, random_number_generator);
265 return particle_handler.insert_particle(
266 position_and_reference_position.first,
267 position_and_reference_position.second,
268 id,
269 cell);
270 }
271
272
273
274 template <int dim, int spacedim>
275 void
278 const Function<spacedim> &probability_density_function,
279 const bool random_cell_selection,
280 const types::particle_index n_particles_to_create,
281 ParticleHandler<dim, spacedim> &particle_handler,
282 const Mapping<dim, spacedim> &mapping,
283 const unsigned int random_number_seed)
284 {
285 unsigned int combined_seed = random_number_seed;
286 if (const auto tria =
287 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
289 {
290 const unsigned int my_rank =
291 Utilities::MPI::this_mpi_process(tria->get_mpi_communicator());
292 combined_seed += my_rank;
293 }
294 std::mt19937 random_number_generator(combined_seed);
295
296 types::particle_index start_particle_id(0);
297 types::particle_index n_local_particles(0);
298
299 std::vector<types::particle_index> particles_per_cell(
300 triangulation.n_active_cells(), 0);
301
302 // First determine how many particles to generate in which cell
303 {
304 // Get the local accumulated probabilities for every cell
305 const std::vector<double> cumulative_cell_weights =
306 compute_local_cumulative_cell_weights(triangulation,
307 mapping,
308 probability_density_function);
309
310 // Sum the local integrals over all nodes
311 double local_weight_integral = (cumulative_cell_weights.size() > 0) ?
312 cumulative_cell_weights.back() :
313 0.0;
314
315
316 double local_start_weight = numbers::signaling_nan<double>();
317 double global_weight_integral = numbers::signaling_nan<double>();
318
319 if (const auto tria =
320 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
322 {
323 std::tie(local_start_weight, global_weight_integral) =
325 local_weight_integral, tria->get_mpi_communicator());
326 }
327 else
328 {
329 local_start_weight = 0;
330 global_weight_integral = local_weight_integral;
331 }
332
333 AssertThrow(global_weight_integral > std::numeric_limits<double>::min(),
335 "The integral of the user prescribed probability "
336 "density function over the domain equals zero, "
337 "deal.II has no way to determine the cell of "
338 "generated particles. Please ensure that the "
339 "provided function is positive in at least a "
340 "part of the domain; also check the syntax of "
341 "the function."));
342
343 // Calculate start id
344 start_particle_id =
345 std::llround(static_cast<double>(n_particles_to_create) *
346 local_start_weight / global_weight_integral);
347
348 // Calculate number of local particles
349 const types::particle_index end_particle_id =
350 std::llround(static_cast<double>(n_particles_to_create) *
351 ((local_start_weight + local_weight_integral) /
352 global_weight_integral));
353 n_local_particles = end_particle_id - start_particle_id;
354
355 if (random_cell_selection)
356 {
357 // Uniform distribution on the interval [0,local_weight_integral).
358 // This will be used to randomly select cells for all local
359 // particles.
360 std::uniform_real_distribution<double> uniform_distribution(
361 0.0, local_weight_integral);
362
363 // Loop over all particles to create locally and pick their cells
364 for (types::particle_index current_particle_index = 0;
365 current_particle_index < n_local_particles;
366 ++current_particle_index)
367 {
368 // Draw the random number that determines the cell of the
369 // particle
370 const double random_weight =
371 uniform_distribution(random_number_generator);
372
373 const std::vector<double>::const_iterator selected_cell =
374 std::lower_bound(cumulative_cell_weights.begin(),
375 cumulative_cell_weights.end(),
376 random_weight);
377 const unsigned int cell_index =
378 std::distance(cumulative_cell_weights.begin(), selected_cell);
379
380 ++particles_per_cell[cell_index];
381 }
382 }
383 else
384 {
385 if (local_weight_integral > 0)
386 {
387 types::particle_index particles_created = 0;
388
389 // Compute number of particles per cell according to the ratio
390 // between their weight and the local weight integral
391 for (const auto &cell : triangulation.active_cell_iterators() |
393 {
394 const types::particle_index cumulative_particles_to_create =
395 std::llround(
396 static_cast<double>(n_local_particles) *
397 cumulative_cell_weights[cell->active_cell_index()] /
398 local_weight_integral);
399
400 // Compute particles for this cell as difference between
401 // number of particles that should be created including this
402 // cell minus the number of particles already created.
403 particles_per_cell[cell->active_cell_index()] =
404 cumulative_particles_to_create - particles_created;
405 particles_created +=
406 particles_per_cell[cell->active_cell_index()];
407 }
408 }
409 }
410 }
411
412 // Now generate as many particles per cell as determined above
413 {
414 particle_handler.reserve(particle_handler.n_locally_owned_particles() +
415 n_local_particles);
416 unsigned int current_particle_index = start_particle_id;
417
418 for (const auto &cell : triangulation.active_cell_iterators() |
420 {
421 for (unsigned int i = 0;
422 i < particles_per_cell[cell->active_cell_index()];
423 ++i)
424 {
426 current_particle_index,
427 random_number_generator,
428 particle_handler,
429 mapping);
430
431 ++current_particle_index;
432 }
433 }
434
435 particle_handler.update_cached_numbers();
436 }
437 }
438
439
440
441 template <int dim, int spacedim>
442 void
444 const std::vector<std::vector<BoundingBox<spacedim>>>
445 &global_bounding_boxes,
446 ParticleHandler<dim, spacedim> &particle_handler,
447 const Mapping<dim, spacedim> &mapping,
448 const ComponentMask &components,
449 const std::vector<std::vector<double>> &properties)
450 {
451 const auto &fe = dof_handler.get_fe();
452
453 // Take care of components
454 const ComponentMask mask =
455 (components.size() == 0 ? ComponentMask(fe.n_components(), true) :
456 components);
457
458 const std::map<types::global_dof_index, Point<spacedim>>
459 support_points_map =
460 DoFTools::map_dofs_to_support_points(mapping, dof_handler, mask);
461
462 // Generate the vector of points from the map
463 // Memory is reserved for efficiency reasons
464 std::vector<Point<spacedim>> support_points_vec;
465 support_points_vec.reserve(support_points_map.size());
466 for (const auto &element : support_points_map)
467 support_points_vec.push_back(element.second);
468
469 particle_handler.insert_global_particles(support_points_vec,
470 global_bounding_boxes,
471 properties);
472 }
473
474
475 template <int dim, int spacedim>
476 void
479 const Quadrature<dim> &quadrature,
480 // const std::vector<Point<dim>> & particle_reference_locations,
481 const std::vector<std::vector<BoundingBox<spacedim>>>
482 &global_bounding_boxes,
483 ParticleHandler<dim, spacedim> &particle_handler,
484 const Mapping<dim, spacedim> &mapping,
485 const std::vector<std::vector<double>> &properties)
486 {
487 const std::vector<Point<dim>> &particle_reference_locations =
488 quadrature.get_points();
489 std::vector<Point<spacedim>> points_to_generate;
490 points_to_generate.reserve(particle_reference_locations.size() *
491 triangulation.n_active_cells());
492
493 // Loop through cells and gather gauss points
494 for (const auto &cell : triangulation.active_cell_iterators() |
496 {
497 for (const auto &reference_location : particle_reference_locations)
498 {
499 const Point<spacedim> position_real =
500 mapping.transform_unit_to_real_cell(cell, reference_location);
501 points_to_generate.push_back(position_real);
502 }
503 }
504 particle_handler.insert_global_particles(points_to_generate,
505 global_bounding_boxes,
506 properties);
507 }
508 } // namespace Generators
509} // namespace Particles
510
511#include "particles/generators.inst"
512
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.
Definition mapping.h:320
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 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)
Definition point.h:113
const std::vector< Point< dim > > & get_points() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
unsigned int cell_index
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
Definition mpi.cc:929
void map_dofs_to_support_points(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, std::vector< Point< spacedim > > &support_points, const ComponentMask &mask={})
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)
Definition mpi.cc:114
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation