Reference documentation for deal.II version GIT e22cb6a53e 2023-09-21 14:00:02+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\}}\)
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 
21 #include <deal.II/dofs/dof_tools.h>
22 
23 #include <deal.II/fe/fe_nothing.h>
24 #include <deal.II/fe/fe_values.h>
25 
27 
29 
30 #include <limits>
31 
33 
34 namespace Particles
35 {
36  namespace Generators
37  {
38  namespace
39  {
43  template <int dim>
45  ProbabilityFunctionNegative,
46  Point<dim>,
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  }
156  catch (typename Mapping<dim>::ExcTransformationFailed &)
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  {
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> *>(
196  &triangulation))
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,
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,
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> *>(
288  &triangulation))
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> *>(
319  &triangulation))
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(),
331  ExcMessage(
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> *>(
347  &triangulation))
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 (const auto &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
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
void reserve(const std::size_t n_particles)
types::particle_index n_locally_owned_particles() const
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={})
particle_iterator insert_particle(const Particle< dim, spacedim > &particle, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
Definition: point.h:112
const std::vector< Point< dim > > & get_points() const
const Point< dim > & point(const unsigned int i) const
virtual MPI_Comm get_communicator() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
unsigned int cell_index
Definition: grid_tools.cc:1192
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1916
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:512
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
@ 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={})
Definition: dof_tools.cc:2293
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)
Definition: generators.cc:276
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={})
Definition: generators.cc:456
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 >()))
Definition: generators.cc:183
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 >()))
Definition: generators.cc:256
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 >()))
Definition: generators.cc:239
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={})
Definition: generators.cc:490
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:164
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1008
unsigned int particle_index
Definition: property_pool.h:65
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
#define DEAL_II_PARTICLE_INDEX_MPI_TYPE
Definition: property_pool.h:72
const ::Triangulation< dim, spacedim > & tria