Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
generators.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2019 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
20
22
25
27
29
30#include <limits>
31
33
34namespace Particles
35{
36 namespace Generators
37 {
38 namespace
39 {
43 template <int dim>
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.");
50
51
52 // This function integrates the given probability density
53 // function over all cells of the triangulation. For each
54 // cell it stores the cumulative sum (of all previous
55 // cells including the current cell) in a vector that is then
56 // returned. Therefore the returned vector has as many entries
57 // as active cells, the first entry being the integral over the
58 // first cell, and the last entry the integral over the whole
59 // locally owned domain. Cells that are not locally owned
60 // simply store the same value as the cell before (equivalent
61 // to assuming a probability density function value of 0).
62 template <int dim, int spacedim = dim>
63 std::vector<double>
64 compute_local_cumulative_cell_weights(
66 const Mapping<dim, spacedim> & mapping,
67 const Function<spacedim> & probability_density_function)
68 {
69 std::vector<double> cumulative_cell_weights(
70 triangulation.n_active_cells());
71 double cumulative_weight = 0.0;
72
73 // Evaluate function at all cell midpoints
74 const QMidpoint<dim> quadrature_formula;
75
76 // In the simplest case we do not even need a FEValues object, because
77 // using cell->center() and cell->measure() would be equivalent. This
78 // fails however for higher-order mappings.
79 FE_Nothing<dim, spacedim> alibi_finite_element;
80 FEValues<dim, spacedim> fe_values(mapping,
81 alibi_finite_element,
82 quadrature_formula,
85
86 // compute the integral weight
87 for (const auto &cell : triangulation.active_cell_iterators())
88 {
89 if (cell->is_locally_owned())
90 {
91 fe_values.reinit(cell);
92 const double quadrature_point_weight =
93 probability_density_function.value(
94 fe_values.get_quadrature_points()[0]);
95
96 AssertThrow(quadrature_point_weight >= 0.0,
97 ProbabilityFunctionNegative<dim>(
98 quadrature_formula.point(0)));
99
100 // get_cell_weight makes sure to return positive values
101 cumulative_weight += quadrature_point_weight * fe_values.JxW(0);
102 }
103 cumulative_cell_weights[cell->active_cell_index()] =
104 cumulative_weight;
105 }
106
107 return cumulative_cell_weights;
108 }
109
110
111
112 // This function generates a random position in the given cell and
113 // returns the position and its coordinates in the reference cell. It
114 // first tries to generate a random and uniformly distributed point in
115 // real space, but if that fails (e.g. because the cell has a bad aspect
116 // ratio) it reverts to generating a random point in the unit cell.
117 template <int dim, int spacedim>
118 std::pair<Point<spacedim>, Point<dim>>
119 random_location_in_cell(
121 const Mapping<dim, spacedim> &mapping,
122 std::mt19937 & random_number_generator)
123 {
124 Assert(cell->reference_cell().is_hyper_cube() == true,
126
127 // Uniform distribution on the interval [0,1]. This
128 // will be used to generate random particle locations.
129 std::uniform_real_distribution<double> uniform_distribution_01(0, 1);
130
131 const BoundingBox<spacedim> cell_bounding_box(cell->bounding_box());
132 const std::pair<Point<spacedim>, Point<spacedim>> &cell_bounds(
133 cell_bounding_box.get_boundary_points());
134
135 // Generate random points in these bounds until one is within the cell
136 // or we exceed the maximum number of attempts.
137 const unsigned int n_attempts = 100;
138 for (unsigned int i = 0; i < n_attempts; ++i)
139 {
140 Point<spacedim> position;
141 for (unsigned int d = 0; d < spacedim; ++d)
142 {
143 position[d] = uniform_distribution_01(random_number_generator) *
144 (cell_bounds.second[d] - cell_bounds.first[d]) +
145 cell_bounds.first[d];
146 }
147
148 try
149 {
150 const Point<dim> reference_position =
151 mapping.transform_real_to_unit_cell(cell, position);
152
153 if (cell->reference_cell().contains_point(reference_position))
154 return {position, reference_position};
155 }
157 {
158 // The point is not in this cell. Do nothing, just try again.
159 }
160 }
161
162 // If the above algorithm has not worked (e.g. because of badly
163 // deformed cells), retry generating particles
164 // randomly within the reference cell and then mapping it to to
165 // real space. This is not generating a
166 // uniform distribution in real space, but will always succeed.
167 Point<dim> reference_position;
168 for (unsigned int d = 0; d < dim; ++d)
169 reference_position[d] =
170 uniform_distribution_01(random_number_generator);
171
172 const Point<spacedim> position =
173 mapping.transform_unit_to_real_cell(cell, reference_position);
174
175 return {position, reference_position};
176 }
177 } // namespace
178
179
180
181 template <int dim, int spacedim>
182 void
185 const std::vector<Point<dim>> & particle_reference_locations,
186 ParticleHandler<dim, spacedim> & particle_handler,
187 const Mapping<dim, spacedim> & mapping)
188 {
189 types::particle_index particle_index = 0;
190 types::particle_index n_particles_to_generate =
191 triangulation.n_active_cells() * particle_reference_locations.size();
192
193#ifdef DEAL_II_WITH_MPI
194 if (const auto tria =
195 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
197 {
198 n_particles_to_generate = tria->n_locally_owned_active_cells() *
199 particle_reference_locations.size();
200
201 // The local particle start index is the number of all particles
202 // generated on lower MPI ranks.
203 const int ierr = MPI_Exscan(&n_particles_to_generate,
204 &particle_index,
205 1,
207 MPI_SUM,
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 =
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 double global_weight_integral;
316
317 if (const auto tria =
318 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
320 {
321 global_weight_integral =
322 Utilities::MPI::sum(local_weight_integral,
324 }
325 else
326 {
327 global_weight_integral = local_weight_integral;
328 }
329
330 AssertThrow(global_weight_integral > std::numeric_limits<double>::min(),
332 "The integral of the user prescribed probability "
333 "density function over the domain equals zero, "
334 "deal.II has no way to determine the cell of "
335 "generated particles. Please ensure that the "
336 "provided function is positive in at least a "
337 "part of the domain; also check the syntax of "
338 "the function."));
339
340 // Determine the starting weight of this process, which is the sum of
341 // the weights of all processes with a lower rank
342 double local_start_weight = 0.0;
343
344#ifdef DEAL_II_WITH_MPI
345 if (const auto tria =
346 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
348 {
349 const int ierr = MPI_Exscan(&local_weight_integral,
350 &local_start_weight,
351 1,
352 MPI_DOUBLE,
353 MPI_SUM,
355 AssertThrowMPI(ierr);
356 }
357#endif
358
359 // Calculate start id
360 start_particle_id =
361 std::llround(static_cast<double>(n_particles_to_create) *
362 local_start_weight / global_weight_integral);
363
364 // Calculate number of local particles
365 const types::particle_index end_particle_id =
366 std::llround(static_cast<double>(n_particles_to_create) *
367 ((local_start_weight + local_weight_integral) /
368 global_weight_integral));
369 n_local_particles = end_particle_id - start_particle_id;
370
371 if (random_cell_selection)
372 {
373 // Uniform distribution on the interval [0,local_weight_integral).
374 // This will be used to randomly select cells for all local
375 // particles.
376 std::uniform_real_distribution<double> uniform_distribution(
377 0.0, local_weight_integral);
378
379 // Loop over all particles to create locally and pick their cells
380 for (types::particle_index current_particle_index = 0;
381 current_particle_index < n_local_particles;
382 ++current_particle_index)
383 {
384 // Draw the random number that determines the cell of the
385 // particle
386 const double random_weight =
387 uniform_distribution(random_number_generator);
388
389 const std::vector<double>::const_iterator selected_cell =
390 std::lower_bound(cumulative_cell_weights.begin(),
391 cumulative_cell_weights.end(),
392 random_weight);
393 const unsigned int cell_index =
394 std::distance(cumulative_cell_weights.begin(), selected_cell);
395
396 ++particles_per_cell[cell_index];
397 }
398 }
399 else
400 {
401 // Compute number of particles per cell according to the ratio
402 // between their weight and the local weight integral
403 types::particle_index particles_created = 0;
404
405 for (const auto &cell : triangulation.active_cell_iterators() |
407 {
408 const types::particle_index cumulative_particles_to_create =
409 std::llround(
410 static_cast<double>(n_local_particles) *
411 cumulative_cell_weights[cell->active_cell_index()] /
412 local_weight_integral);
413
414 // Compute particles for this cell as difference between
415 // number of particles that should be created including this
416 // cell minus the number of particles already created.
417 particles_per_cell[cell->active_cell_index()] =
418 cumulative_particles_to_create - particles_created;
419 particles_created +=
420 particles_per_cell[cell->active_cell_index()];
421 }
422 }
423 }
424
425 // Now generate as many particles per cell as determined above
426 {
427 particle_handler.reserve(particle_handler.n_locally_owned_particles() +
428 n_local_particles);
429 unsigned int current_particle_index = start_particle_id;
430
431 for (const auto &cell : triangulation.active_cell_iterators() |
433 {
434 for (unsigned int i = 0;
435 i < particles_per_cell[cell->active_cell_index()];
436 ++i)
437 {
439 current_particle_index,
440 random_number_generator,
441 particle_handler,
442 mapping);
443
444 ++current_particle_index;
445 }
446 }
447
448 particle_handler.update_cached_numbers();
449 }
450 }
451
452
453
454 template <int dim, int spacedim>
455 void
457 const std::vector<std::vector<BoundingBox<spacedim>>>
458 & global_bounding_boxes,
459 ParticleHandler<dim, spacedim> &particle_handler,
460 const Mapping<dim, spacedim> & mapping,
461 const ComponentMask & components,
462 const std::vector<std::vector<double>> &properties)
463 {
464 const auto &fe = dof_handler.get_fe();
465
466 // Take care of components
467 const ComponentMask mask =
468 (components.size() == 0 ? ComponentMask(fe.n_components(), true) :
469 components);
470
471 const std::map<types::global_dof_index, Point<spacedim>>
472 support_points_map =
473 DoFTools::map_dofs_to_support_points(mapping, dof_handler, mask);
474
475 // Generate the vector of points from the map
476 // Memory is reserved for efficiency reasons
477 std::vector<Point<spacedim>> support_points_vec;
478 support_points_vec.reserve(support_points_map.size());
479 for (auto const &element : support_points_map)
480 support_points_vec.push_back(element.second);
481
482 particle_handler.insert_global_particles(support_points_vec,
483 global_bounding_boxes,
484 properties);
485 }
486
487
488 template <int dim, int spacedim>
489 void
492 const Quadrature<dim> & quadrature,
493 // const std::vector<Point<dim>> & particle_reference_locations,
494 const std::vector<std::vector<BoundingBox<spacedim>>>
495 & global_bounding_boxes,
496 ParticleHandler<dim, spacedim> & particle_handler,
497 const Mapping<dim, spacedim> & mapping,
498 const std::vector<std::vector<double>> &properties)
499 {
500 const std::vector<Point<dim>> &particle_reference_locations =
501 quadrature.get_points();
502 std::vector<Point<spacedim>> points_to_generate;
503 points_to_generate.reserve(particle_reference_locations.size() *
504 triangulation.n_active_cells());
505
506 // Loop through cells and gather gauss points
507 for (const auto &cell : triangulation.active_cell_iterators() |
509 {
510 for (const auto &reference_location : particle_reference_locations)
511 {
512 const Point<spacedim> position_real =
513 mapping.transform_unit_to_real_cell(cell, reference_location);
514 points_to_generate.push_back(position_real);
515 }
516 }
517 particle_handler.insert_global_particles(points_to_generate,
518 global_bounding_boxes,
519 properties);
520 }
521 } // namespace Generators
522} // namespace Particles
523
524#include "generators.inst"
525
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:317
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:112
const Point< dim > & point(const unsigned int i) const
const std::vector< Point< dim > > & get_points() const
virtual MPI_Comm get_communicator() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertThrowMPI(error_code)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
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=ComponentMask())
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)
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=ComponentMask(), const std::vector< std::vector< double > > &properties={})
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 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)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
#define DEAL_II_PARTICLE_INDEX_MPI_TYPE
const ::Triangulation< dim, spacedim > & tria