Reference documentation for deal.II version 9.2.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\}}\)
generators.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2019 - 2020 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 
29 
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  } // namespace
110 
111  template <int dim, int spacedim>
112  void
114  const Triangulation<dim, spacedim> &triangulation,
115  const std::vector<Point<dim>> & particle_reference_locations,
116  ParticleHandler<dim, spacedim> & particle_handler,
117  const Mapping<dim, spacedim> & mapping)
118  {
120 
121 #ifdef DEAL_II_WITH_MPI
122  if (const auto tria =
123  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
124  &triangulation))
125  {
126  const types::particle_index n_particles_to_generate =
127  tria->n_locally_owned_active_cells() *
128  particle_reference_locations.size();
129 
130  // The local particle start index is the number of all particles
131  // generated on lower MPI ranks.
132  MPI_Exscan(&n_particles_to_generate,
133  &particle_index,
134  1,
136  MPI_SUM,
137  tria->get_communicator());
138  }
139 #endif
140 
141  for (const auto &cell : triangulation.active_cell_iterators())
142  {
143  if (cell->is_locally_owned())
144  {
145  for (const auto &reference_location :
146  particle_reference_locations)
147  {
148  const Point<spacedim> position_real =
149  mapping.transform_unit_to_real_cell(cell,
150  reference_location);
151 
152  const Particle<dim, spacedim> particle(position_real,
153  reference_location,
154  particle_index);
155  particle_handler.insert_particle(particle, cell);
156  ++particle_index;
157  }
158  }
159  }
160 
161  particle_handler.update_cached_numbers();
162  }
163 
164 
165 
166  template <int dim, int spacedim>
170  const types::particle_index id,
171  std::mt19937 & random_number_generator,
172  const Mapping<dim, spacedim> &mapping)
173  {
174  // Uniform distribution on the interval [0,1]. This
175  // will be used to generate random particle locations.
176  std::uniform_real_distribution<double> uniform_distribution_01(0, 1);
177 
178  const BoundingBox<spacedim> cell_bounding_box(cell->bounding_box());
179  const std::pair<Point<spacedim>, Point<spacedim>> cell_bounds(
180  cell_bounding_box.get_boundary_points());
181 
182  // Generate random points in these bounds until one is within the cell
183  unsigned int iteration = 0;
184  const unsigned int maximum_iterations = 100;
185  Point<spacedim> particle_position;
186  while (iteration < maximum_iterations)
187  {
188  for (unsigned int d = 0; d < spacedim; ++d)
189  {
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];
194  }
195  try
196  {
197  const Point<dim> p_unit =
198  mapping.transform_real_to_unit_cell(cell, particle_position);
200  {
201  // Generate the particle
202  return Particle<dim, spacedim>(particle_position, p_unit, id);
203  }
204  }
205  catch (typename Mapping<dim>::ExcTransformationFailed &)
206  {
207  // The point is not in this cell. Do nothing, just try again.
208  }
209  ++iteration;
210  }
211  AssertThrow(
212  iteration < maximum_iterations,
213  ExcMessage(
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: " +
218  cell->measure() /
219  (cell_bounds.second - cell_bounds.first).norm_square())));
220 
221  return Particle<dim, spacedim>();
222  }
223 
224 
225 
226  template <int dim, int spacedim>
227  void
229  const Triangulation<dim, spacedim> &triangulation,
230  const Function<spacedim> & probability_density_function,
231  const bool random_cell_selection,
232  const types::particle_index n_particles_to_create,
233  ParticleHandler<dim, spacedim> & particle_handler,
234  const Mapping<dim, spacedim> & mapping,
235  const unsigned int random_number_seed)
236  {
237  unsigned int combined_seed = random_number_seed;
238  if (const auto tria =
239  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
240  &triangulation))
241  {
242  const unsigned int my_rank =
243  Utilities::MPI::this_mpi_process(tria->get_communicator());
244  combined_seed += my_rank;
245  }
246  std::mt19937 random_number_generator(combined_seed);
247 
248  types::particle_index start_particle_id(0);
249  types::particle_index n_local_particles(0);
250 
251  std::vector<types::particle_index> particles_per_cell(
252  triangulation.n_active_cells(), 0);
253 
254  // First determine how many particles to generate in which cell
255  {
256  // Get the local accumulated probabilities for every cell
257  const std::vector<double> cumulative_cell_weights =
258  compute_local_cumulative_cell_weights(triangulation,
259  mapping,
260  probability_density_function);
261 
262  // Sum the local integrals over all nodes
263  double local_weight_integral = cumulative_cell_weights.back();
264  double global_weight_integral;
265 
266  if (const auto tria =
267  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
268  &triangulation))
269  {
270  global_weight_integral =
271  Utilities::MPI::sum(local_weight_integral,
272  tria->get_communicator());
273  }
274  else
275  {
276  global_weight_integral = local_weight_integral;
277  }
278 
279  AssertThrow(global_weight_integral > std::numeric_limits<double>::min(),
280  ExcMessage(
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 "
287  "the function."));
288 
289  // Determine the starting weight of this process, which is the sum of
290  // the weights of all processes with a lower rank
291  double local_start_weight = 0.0;
292 
293 #ifdef DEAL_II_WITH_MPI
294  if (const auto tria =
295  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
296  &triangulation))
297  {
298  MPI_Exscan(&local_weight_integral,
299  &local_start_weight,
300  1,
301  MPI_DOUBLE,
302  MPI_SUM,
303  tria->get_communicator());
304  }
305 #endif
306 
307  // Calculate start id
308  start_particle_id =
309  std::llround(static_cast<double>(n_particles_to_create) *
310  local_start_weight / global_weight_integral);
311 
312  // Calcualate number of local particles
313  const types::particle_index end_particle_id =
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;
318 
319  if (random_cell_selection)
320  {
321  // Uniform distribution on the interval [0,local_weight_integral).
322  // This will be used to randomly select cells for all local
323  // particles.
324  std::uniform_real_distribution<double> uniform_distribution(
325  0.0, local_weight_integral);
326 
327  // Loop over all particles to create locally and pick their cells
328  for (types::particle_index current_particle_index = 0;
329  current_particle_index < n_local_particles;
330  ++current_particle_index)
331  {
332  // Draw the random number that determines the cell of the
333  // particle
334  const double random_weight =
335  uniform_distribution(random_number_generator);
336 
337  const std::vector<double>::const_iterator selected_cell =
338  std::lower_bound(cumulative_cell_weights.begin(),
339  cumulative_cell_weights.end(),
340  random_weight);
341  const unsigned int cell_index =
342  std::distance(cumulative_cell_weights.begin(), selected_cell);
343 
344  ++particles_per_cell[cell_index];
345  }
346  }
347  else
348  {
349  // Compute number of particles per cell according to the ratio
350  // between their weight and the local weight integral
351  types::particle_index particles_created = 0;
352 
353  for (const auto &cell : triangulation.active_cell_iterators())
354  if (cell->is_locally_owned())
355  {
356  const types::particle_index cumulative_particles_to_create =
357  std::llround(
358  static_cast<double>(n_local_particles) *
359  cumulative_cell_weights[cell->active_cell_index()] /
360  local_weight_integral);
361 
362  // Compute particles for this cell as difference between
363  // number of particles that should be created including this
364  // cell minus the number of particles already created.
365  particles_per_cell[cell->active_cell_index()] =
366  cumulative_particles_to_create - particles_created;
367  particles_created +=
368  particles_per_cell[cell->active_cell_index()];
369  }
370  }
371  }
372 
373  // Now generate as many particles per cell as determined above
374  {
375  unsigned int current_particle_index = start_particle_id;
376 
377  std::multimap<
380  particles;
381 
382  for (const auto &cell : triangulation.active_cell_iterators())
383  if (cell->is_locally_owned())
384  {
385  for (unsigned int i = 0;
386  i < particles_per_cell[cell->active_cell_index()];
387  ++i)
388  {
389  Particle<dim, spacedim> particle =
391  current_particle_index,
392  random_number_generator,
393  mapping);
394  particles.emplace_hint(particles.end(),
395  cell,
396  std::move(particle));
397  ++current_particle_index;
398  }
399  }
400 
401  particle_handler.insert_particles(particles);
402  }
403  }
404 
405  template <int dim, int spacedim>
406  void
408  const std::vector<std::vector<BoundingBox<spacedim>>>
409  & global_bounding_boxes,
410  ParticleHandler<dim, spacedim> &particle_handler,
411  const Mapping<dim, spacedim> & mapping,
412  const ComponentMask & components)
413  {
414  const auto &fe = dof_handler.get_fe();
415 
416  // Take care of components
417  const ComponentMask mask =
418  (components.size() == 0 ? ComponentMask(fe.n_components(), true) :
419  components);
420 
421  std::map<types::global_dof_index, Point<spacedim>> support_points_map;
422 
424  dof_handler,
425  support_points_map,
426  mask);
427 
428  // Generate the vector of points from the map
429  // Memory is reserved for efficiency reasons
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);
434 
435  particle_handler.insert_global_particles(support_points_vec,
436  global_bounding_boxes);
437  }
438 
439 
440  template <int dim, int spacedim>
441  void
443  const Triangulation<dim, spacedim> &triangulation,
444  const Quadrature<dim> & quadrature,
445  // const std::vector<Point<dim>> & particle_reference_locations,
446  const std::vector<std::vector<BoundingBox<spacedim>>>
447  & global_bounding_boxes,
448  ParticleHandler<dim, spacedim> &particle_handler,
449  const Mapping<dim, spacedim> & mapping)
450  {
451  const std::vector<Point<dim>> &particle_reference_locations =
452  quadrature.get_points();
453  std::vector<Point<spacedim>> points_to_generate;
454 
455  // Loop through cells and gather gauss points
456  for (const auto &cell : triangulation.active_cell_iterators())
457  {
458  if (cell->is_locally_owned())
459  {
460  for (const auto &reference_location :
461  particle_reference_locations)
462  {
463  const Point<spacedim> position_real =
464  mapping.transform_unit_to_real_cell(cell,
465  reference_location);
466  points_to_generate.push_back(position_real);
467  }
468  }
469  }
470  particle_handler.insert_global_particles(points_to_generate,
471  global_bounding_boxes);
472  }
473  } // namespace Generators
474 } // namespace Particles
475 
476 #include "generators.inst"
477 
Transformed quadrature weights.
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1102
unsigned int n_active_cells() const
Definition: tria.cc:12675
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=StaticMappingQ1< dim, spacedim >::mapping)
Definition: generators.cc:113
const std::vector< Point< dim > > & get_points() const
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())
Definition: dof_tools.cc:2304
particle_iterator insert_particle(const Particle< dim, spacedim > &particle, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12185
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=StaticMappingQ1< dim, spacedim >::mapping)
Definition: generators.cc:442
const Point< dim > & point(const unsigned int i) const
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=StaticMappingQ1< dim, spacedim >::mapping, const unsigned int random_number_seed=5432)
Definition: generators.cc:228
Transformed quadrature points.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
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=StaticMappingQ1< dim, spacedim >::mapping)
Definition: generators.cc:168
unsigned int particle_index
Definition: particle.h:69
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
void insert_particles(const std::multimap< typename Triangulation< dim, spacedim >::active_cell_iterator, Particle< dim, spacedim >> &particles)
T sum(const T &t, const MPI_Comm &mpi_communicator)
Abstract base class for mapping classes.
Definition: mapping.h:302
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
std::string to_string(const T &t)
Definition: patterns.h:2360
unsigned int size() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
#define DEAL_II_PARTICLE_INDEX_MPI_TYPE
Definition: particle.h:76
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
T min(const T &t, const MPI_Comm &mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
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=StaticMappingQ1< dim, spacedim >::mapping, const ComponentMask &components=ComponentMask())
Definition: generators.cc:407
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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={})
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0