Reference documentation for deal.II version Git 7439f2b0b4 2021-09-15 12:22:58 +0200
\(\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 - 2021 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::TriangulationBase<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  const int ierr = MPI_Exscan(&n_particles_to_generate,
133  &particle_index,
134  1,
136  MPI_SUM,
137  tria->get_communicator());
138  AssertThrowMPI(ierr);
139  }
140 #endif
141 
142  for (const auto &cell : triangulation.active_cell_iterators())
143  {
144  if (cell->is_locally_owned())
145  {
146  for (const auto &reference_location :
147  particle_reference_locations)
148  {
149  const Point<spacedim> position_real =
150  mapping.transform_unit_to_real_cell(cell,
151  reference_location);
152 
153  const Particle<dim, spacedim> particle(position_real,
154  reference_location,
155  particle_index);
156  particle_handler.insert_particle(particle, cell);
157  ++particle_index;
158  }
159  }
160  }
161 
162  particle_handler.update_cached_numbers();
163  }
164 
165 
166 
167  template <int dim, int spacedim>
171  const types::particle_index id,
172  std::mt19937 & random_number_generator,
173  const Mapping<dim, spacedim> &mapping)
174  {
175  // Uniform distribution on the interval [0,1]. This
176  // will be used to generate random particle locations.
177  std::uniform_real_distribution<double> uniform_distribution_01(0, 1);
178 
179  const BoundingBox<spacedim> cell_bounding_box(cell->bounding_box());
180  const std::pair<Point<spacedim>, Point<spacedim>> &cell_bounds(
181  cell_bounding_box.get_boundary_points());
182 
183  // Generate random points in these bounds until one is within the cell
184  unsigned int iteration = 0;
185  const unsigned int maximum_iterations = 100;
186  Point<spacedim> particle_position;
187  while (iteration < maximum_iterations)
188  {
189  for (unsigned int d = 0; d < spacedim; ++d)
190  {
191  particle_position[d] =
192  uniform_distribution_01(random_number_generator) *
193  (cell_bounds.second[d] - cell_bounds.first[d]) +
194  cell_bounds.first[d];
195  }
196  try
197  {
198  const Point<dim> p_unit =
199  mapping.transform_real_to_unit_cell(cell, particle_position);
201  {
202  // Generate the particle
203  return Particle<dim, spacedim>(particle_position, p_unit, id);
204  }
205  }
206  catch (typename Mapping<dim>::ExcTransformationFailed &)
207  {
208  // The point is not in this cell. Do nothing, just try again.
209  }
210  ++iteration;
211  }
212  AssertThrow(
213  iteration < maximum_iterations,
214  ExcMessage(
215  "Couldn't generate a particle position within the maximum number of tries. "
216  "The ratio between the bounding box volume in which the particle is "
217  "generated and the actual cell volume is approximately: " +
219  cell->measure() /
220  (cell_bounds.second - cell_bounds.first).norm_square())));
221 
222  return Particle<dim, spacedim>();
223  }
224 
225 
226 
227  template <int dim, int spacedim>
228  void
230  const Triangulation<dim, spacedim> &triangulation,
231  const Function<spacedim> & probability_density_function,
232  const bool random_cell_selection,
233  const types::particle_index n_particles_to_create,
234  ParticleHandler<dim, spacedim> & particle_handler,
235  const Mapping<dim, spacedim> & mapping,
236  const unsigned int random_number_seed)
237  {
238  unsigned int combined_seed = random_number_seed;
239  if (const auto tria =
240  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
241  &triangulation))
242  {
243  const unsigned int my_rank =
244  Utilities::MPI::this_mpi_process(tria->get_communicator());
245  combined_seed += my_rank;
246  }
247  std::mt19937 random_number_generator(combined_seed);
248 
249  types::particle_index start_particle_id(0);
250  types::particle_index n_local_particles(0);
251 
252  std::vector<types::particle_index> particles_per_cell(
253  triangulation.n_active_cells(), 0);
254 
255  // First determine how many particles to generate in which cell
256  {
257  // Get the local accumulated probabilities for every cell
258  const std::vector<double> cumulative_cell_weights =
259  compute_local_cumulative_cell_weights(triangulation,
260  mapping,
261  probability_density_function);
262 
263  // Sum the local integrals over all nodes
264  double local_weight_integral = (cumulative_cell_weights.size() > 0) ?
265  cumulative_cell_weights.back() :
266  0.0;
267 
268  double global_weight_integral;
269 
270  if (const auto tria =
271  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
272  &triangulation))
273  {
274  global_weight_integral =
275  Utilities::MPI::sum(local_weight_integral,
276  tria->get_communicator());
277  }
278  else
279  {
280  global_weight_integral = local_weight_integral;
281  }
282 
283  AssertThrow(global_weight_integral > std::numeric_limits<double>::min(),
284  ExcMessage(
285  "The integral of the user prescribed probability "
286  "density function over the domain equals zero, "
287  "deal.II has no way to determine the cell of "
288  "generated particles. Please ensure that the "
289  "provided function is positive in at least a "
290  "part of the domain; also check the syntax of "
291  "the function."));
292 
293  // Determine the starting weight of this process, which is the sum of
294  // the weights of all processes with a lower rank
295  double local_start_weight = 0.0;
296 
297 #ifdef DEAL_II_WITH_MPI
298  if (const auto tria =
299  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
300  &triangulation))
301  {
302  const int ierr = MPI_Exscan(&local_weight_integral,
303  &local_start_weight,
304  1,
305  MPI_DOUBLE,
306  MPI_SUM,
307  tria->get_communicator());
308  AssertThrowMPI(ierr);
309  }
310 #endif
311 
312  // Calculate start id
313  start_particle_id =
314  std::llround(static_cast<double>(n_particles_to_create) *
315  local_start_weight / global_weight_integral);
316 
317  // Calcualate number of local particles
318  const types::particle_index end_particle_id =
319  std::llround(static_cast<double>(n_particles_to_create) *
320  ((local_start_weight + local_weight_integral) /
321  global_weight_integral));
322  n_local_particles = end_particle_id - start_particle_id;
323 
324  if (random_cell_selection)
325  {
326  // Uniform distribution on the interval [0,local_weight_integral).
327  // This will be used to randomly select cells for all local
328  // particles.
329  std::uniform_real_distribution<double> uniform_distribution(
330  0.0, local_weight_integral);
331 
332  // Loop over all particles to create locally and pick their cells
333  for (types::particle_index current_particle_index = 0;
334  current_particle_index < n_local_particles;
335  ++current_particle_index)
336  {
337  // Draw the random number that determines the cell of the
338  // particle
339  const double random_weight =
340  uniform_distribution(random_number_generator);
341 
342  const std::vector<double>::const_iterator selected_cell =
343  std::lower_bound(cumulative_cell_weights.begin(),
344  cumulative_cell_weights.end(),
345  random_weight);
346  const unsigned int cell_index =
347  std::distance(cumulative_cell_weights.begin(), selected_cell);
348 
349  ++particles_per_cell[cell_index];
350  }
351  }
352  else
353  {
354  // Compute number of particles per cell according to the ratio
355  // between their weight and the local weight integral
356  types::particle_index particles_created = 0;
357 
358  for (const auto &cell : triangulation.active_cell_iterators())
359  if (cell->is_locally_owned())
360  {
361  const types::particle_index cumulative_particles_to_create =
362  std::llround(
363  static_cast<double>(n_local_particles) *
364  cumulative_cell_weights[cell->active_cell_index()] /
365  local_weight_integral);
366 
367  // Compute particles for this cell as difference between
368  // number of particles that should be created including this
369  // cell minus the number of particles already created.
370  particles_per_cell[cell->active_cell_index()] =
371  cumulative_particles_to_create - particles_created;
372  particles_created +=
373  particles_per_cell[cell->active_cell_index()];
374  }
375  }
376  }
377 
378  // Now generate as many particles per cell as determined above
379  {
380  unsigned int current_particle_index = start_particle_id;
381 
382  std::multimap<
385  particles;
386 
387  for (const auto &cell : triangulation.active_cell_iterators())
388  if (cell->is_locally_owned())
389  {
390  for (unsigned int i = 0;
391  i < particles_per_cell[cell->active_cell_index()];
392  ++i)
393  {
394  Particle<dim, spacedim> particle =
396  current_particle_index,
397  random_number_generator,
398  mapping);
399  particles.emplace_hint(particles.end(),
400  cell,
401  std::move(particle));
402  ++current_particle_index;
403  }
404  }
405 
406  particle_handler.insert_particles(particles);
407  }
408  }
409 
410 
411 
412  template <int dim, int spacedim>
413  void
415  const std::vector<std::vector<BoundingBox<spacedim>>>
416  & global_bounding_boxes,
417  ParticleHandler<dim, spacedim> &particle_handler,
418  const Mapping<dim, spacedim> & mapping,
419  const ComponentMask & components,
420  const std::vector<std::vector<double>> &properties)
421  {
422  const auto &fe = dof_handler.get_fe();
423 
424  // Take care of components
425  const ComponentMask mask =
426  (components.size() == 0 ? ComponentMask(fe.n_components(), true) :
427  components);
428 
429  std::map<types::global_dof_index, Point<spacedim>> support_points_map;
430 
432  dof_handler,
433  support_points_map,
434  mask);
435 
436  // Generate the vector of points from the map
437  // Memory is reserved for efficiency reasons
438  std::vector<Point<spacedim>> support_points_vec;
439  support_points_vec.reserve(support_points_map.size());
440  for (auto const &element : support_points_map)
441  support_points_vec.push_back(element.second);
442 
443  particle_handler.insert_global_particles(support_points_vec,
444  global_bounding_boxes,
445  properties);
446  }
447 
448 
449  template <int dim, int spacedim>
450  void
452  const Triangulation<dim, spacedim> &triangulation,
453  const Quadrature<dim> & quadrature,
454  // const std::vector<Point<dim>> & particle_reference_locations,
455  const std::vector<std::vector<BoundingBox<spacedim>>>
456  & global_bounding_boxes,
457  ParticleHandler<dim, spacedim> & particle_handler,
458  const Mapping<dim, spacedim> & mapping,
459  const std::vector<std::vector<double>> &properties)
460  {
461  const std::vector<Point<dim>> &particle_reference_locations =
462  quadrature.get_points();
463  std::vector<Point<spacedim>> points_to_generate;
464 
465  // Loop through cells and gather gauss points
466  for (const auto &cell : triangulation.active_cell_iterators())
467  {
468  if (cell->is_locally_owned())
469  {
470  for (const auto &reference_location :
471  particle_reference_locations)
472  {
473  const Point<spacedim> position_real =
474  mapping.transform_unit_to_real_cell(cell,
475  reference_location);
476  points_to_generate.push_back(position_real);
477  }
478  }
479  }
480  particle_handler.insert_global_particles(points_to_generate,
481  global_bounding_boxes,
482  properties);
483  }
484  } // namespace Generators
485 } // namespace Particles
486 
487 #include "generators.inst"
488 
Transformed quadrature weights.
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1111
unsigned int n_active_cells() const
Definition: tria.cc:13754
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:2222
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:451
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:13264
const ::Triangulation< dim, spacedim > & tria
const Point< dim > & point(const unsigned int i) const
Transformed quadrature points.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
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
#define DEAL_II_PARTICLE_INDEX_MPI_TYPE
Definition: property_pool.h:72
unsigned int particle_index
Definition: property_pool.h:65
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
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:303
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={})
Definition: generators.cc:414
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
std::string to_string(const T &t)
Definition: patterns.h:2329
unsigned int size() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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:113
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1776
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
Definition: fe.h:38
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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:169
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) 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=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const unsigned int random_number_seed=5432)
Definition: generators.cc:229
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
unsigned int cell_index
Definition: grid_tools.cc:1090
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={})