Reference documentation for deal.II version GIT relicensing-1399-g2b3a69e058 2024-08-09 08:50:02+00:00
Searching...
No Matches
generators.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
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>
46 << "Your probability density function in the particle generator "
47 "returned a negative value for the following position: "
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>
65 const Mapping<dim, spacedim> &mapping,
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
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.
79 FEValues<dim, spacedim> fe_values(mapping,
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);
94
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()] =
104 }
105
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>>
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 {
143 (cell_bounds.second[d] - cell_bounds.first[d]) +
144 cell_bounds.first[d];
145 }
146
147 try
148 {
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.
167 for (unsigned int d = 0; d < dim; ++d)
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,
186 const Mapping<dim, spacedim> &mapping)
187 {
188 types::particle_index particle_index = 0;
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() *
199
200 // The local particle start index is the number of all particles
201 // generated on lower MPI ranks.
203 &particle_index,
204 1,
206 MPI_SUM,
207 tria->get_communicator());
209 }
210#endif
211
212 particle_handler.reserve(particle_handler.n_locally_owned_particles() +
214
215 for (const auto &cell : triangulation.active_cell_iterators() |
217 {
219 {
221 mapping.transform_unit_to_real_cell(cell, reference_location);
222
223 particle_handler.insert_particle(position_real,
225 particle_index,
226 cell);
227 ++particle_index;
228 }
229 }
230
231 particle_handler.update_cached_numbers();
232 }
233
234
235
236 template <int dim, int spacedim>
250
251
252
253 template <int dim, int spacedim>
257 const types::particle_index id,
258 std::mt19937 &random_number_generator,
260 const Mapping<dim, spacedim> &mapping)
261 {
264 return particle_handler.insert_particle(
267 id,
268 cell);
269 }
270
271
272
273 template <int dim, int spacedim>
274 void
278 const bool random_cell_selection,
281 const Mapping<dim, spacedim> &mapping,
282 const unsigned int random_number_seed)
283 {
284 unsigned int combined_seed = random_number_seed;
285 if (const auto tria =
286 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
288 {
289 const unsigned int my_rank =
290 Utilities::MPI::this_mpi_process(tria->get_communicator());
291 combined_seed += my_rank;
292 }
294
297
298 std::vector<types::particle_index> particles_per_cell(
299 triangulation.n_active_cells(), 0);
300
301 // First determine how many particles to generate in which cell
302 {
303 // Get the local accumulated probabilities for every cell
304 const std::vector<double> cumulative_cell_weights =
306 mapping,
308
309 // Sum the local integrals over all nodes
310 double local_weight_integral = (cumulative_cell_weights.size() > 0) ?
312 0.0;
313
314
315 double local_start_weight = numbers::signaling_nan<double>();
316 double global_weight_integral = numbers::signaling_nan<double>();
317
318 if (const auto tria =
319 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
321 {
324 tria->get_communicator());
325 }
326 else
327 {
330 }
331
332 AssertThrow(global_weight_integral > std::numeric_limits<double>::min(),
334 "The integral of the user prescribed probability "
335 "density function over the domain equals zero, "
336 "deal.II has no way to determine the cell of "
337 "generated particles. Please ensure that the "
338 "provided function is positive in at least a "
339 "part of the domain; also check the syntax of "
340 "the function."));
341
342 // Calculate start id
344 std::llround(static_cast<double>(n_particles_to_create) *
346
347 // Calculate number of local particles
349 std::llround(static_cast<double>(n_particles_to_create) *
353
355 {
356 // Uniform distribution on the interval [0,local_weight_integral).
357 // This will be used to randomly select cells for all local
358 // particles.
359 std::uniform_real_distribution<double> uniform_distribution(
361
362 // Loop over all particles to create locally and pick their cells
366 {
367 // Draw the random number that determines the cell of the
368 // particle
369 const double random_weight =
371
372 const std::vector<double>::const_iterator selected_cell =
373 std::lower_bound(cumulative_cell_weights.begin(),
376 const unsigned int cell_index =
377 std::distance(cumulative_cell_weights.begin(), selected_cell);
378
380 }
381 }
382 else
383 {
384 // Compute number of particles per cell according to the ratio
385 // between their weight and the local weight integral
387
388 for (const auto &cell : triangulation.active_cell_iterators() |
390 {
392 std::llround(
393 static_cast<double>(n_local_particles) *
394 cumulative_cell_weights[cell->active_cell_index()] /
396
397 // Compute particles for this cell as difference between
398 // number of particles that should be created including this
399 // cell minus the number of particles already created.
400 particles_per_cell[cell->active_cell_index()] =
403 particles_per_cell[cell->active_cell_index()];
404 }
405 }
406 }
407
408 // Now generate as many particles per cell as determined above
409 {
410 particle_handler.reserve(particle_handler.n_locally_owned_particles() +
413
414 for (const auto &cell : triangulation.active_cell_iterators() |
416 {
417 for (unsigned int i = 0;
418 i < particles_per_cell[cell->active_cell_index()];
419 ++i)
420 {
425 mapping);
426
428 }
429 }
430
431 particle_handler.update_cached_numbers();
432 }
433 }
434
435
436
437 template <int dim, int spacedim>
438 void
440 const std::vector<std::vector<BoundingBox<spacedim>>>
443 const Mapping<dim, spacedim> &mapping,
445 const std::vector<std::vector<double>> &properties)
446 {
447 const auto &fe = dof_handler.get_fe();
448
449 // Take care of components
451 (components.size() == 0 ? ComponentMask(fe.n_components(), true) :
452 components);
453
454 const std::map<types::global_dof_index, Point<spacedim>>
457
458 // Generate the vector of points from the map
459 // Memory is reserved for efficiency reasons
460 std::vector<Point<spacedim>> support_points_vec;
462 for (const auto &element : support_points_map)
463 support_points_vec.push_back(element.second);
464
465 particle_handler.insert_global_particles(support_points_vec,
467 properties);
468 }
469
470
471 template <int dim, int spacedim>
472 void
476 // const std::vector<Point<dim>> & particle_reference_locations,
477 const std::vector<std::vector<BoundingBox<spacedim>>>
480 const Mapping<dim, spacedim> &mapping,
481 const std::vector<std::vector<double>> &properties)
482 {
483 const std::vector<Point<dim>> &particle_reference_locations =
485 std::vector<Point<spacedim>> points_to_generate;
487 triangulation.n_active_cells());
488
489 // Loop through cells and gather gauss points
490 for (const auto &cell : triangulation.active_cell_iterators() |
492 {
494 {
496 mapping.transform_unit_to_real_cell(cell, reference_location);
498 }
499 }
500 particle_handler.insert_global_particles(points_to_generate,
502 properties);
503 }
504 } // namespace Generators
505} // namespace Particles
506
507#include "generators.inst"
508
unsigned int size() const
Abstract base class for mapping classes.
Definition mapping.h:318
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertThrowMPI(error_code)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_JxW_values