Reference documentation for deal.II version GIT c14369f203 2023-10-01 07:40: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\}}\)
grid_refinement.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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 
16 
17 #include <deal.II/base/utilities.h>
18 
22 #include <deal.II/lac/vector.h>
23 
24 #ifdef DEAL_II_WITH_P4EST
25 
27 
30 # include <deal.II/grid/tria.h>
33 
34 # include <algorithm>
35 # include <functional>
36 # include <limits>
37 # include <numeric>
38 
39 
41 
42 
43 namespace
44 {
45  template <int dim, int spacedim>
46  unsigned int
47  n_locally_owned_active_cells(const Triangulation<dim, spacedim> &tria)
48  {
49  if (const auto parallel_tria =
50  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
51  &tria))
52  return parallel_tria->n_locally_owned_active_cells();
53  else
54  return tria.n_active_cells();
55  }
56 
57  template <typename number>
58  inline number
59  max_element(const ::Vector<number> &criteria)
60  {
61  return (criteria.size() > 0) ?
62  (*std::max_element(criteria.begin(), criteria.end())) :
64  }
65 
66 
67 
68  template <typename number>
69  inline number
70  min_element(const ::Vector<number> &criteria)
71  {
72  return (criteria.size() > 0) ?
73  (*std::min_element(criteria.begin(), criteria.end())) :
75  }
76 
77 
78 
84  template <typename number>
85  double
86  compute_global_sum(const ::Vector<number> &criteria,
87  const MPI_Comm mpi_communicator)
88  {
89  double my_sum =
90  std::accumulate(criteria.begin(),
91  criteria.end(),
92  /* do accumulation in the correct data type: */
93  number());
94 
95  double result = 0;
96  // compute the minimum on processor zero
97  const int ierr =
98  MPI_Reduce(&my_sum, &result, 1, MPI_DOUBLE, MPI_SUM, 0, mpi_communicator);
99  AssertThrowMPI(ierr);
100 
101  // make sure only processor zero got something
102  if (Utilities::MPI::this_mpi_process(mpi_communicator) != 0)
103  Assert(result == 0, ExcInternalError());
104 
105  return result;
106  }
107 
108 
109 
114  template <int dim, int spacedim, typename Number>
115  void
116  get_locally_owned_indicators(const ::Triangulation<dim, spacedim> &tria,
117  const ::Vector<Number> &criteria,
118  Vector<Number> &locally_owned_indicators)
119  {
120  Assert(locally_owned_indicators.size() ==
121  n_locally_owned_active_cells(tria),
122  ExcInternalError());
123 
124  unsigned int owned_index = 0;
125  for (const auto &cell :
127  {
128  locally_owned_indicators(owned_index) =
129  criteria(cell->active_cell_index());
130  ++owned_index;
131  }
132  Assert(owned_index == n_locally_owned_active_cells(tria),
133  ExcInternalError());
134  }
135 
136 
137  // we compute refinement thresholds by bisection of the interval spanned by
138  // the smallest and largest error indicator. this leads to a small problem:
139  // if, for example, we want to coarsen zero per cent of the cells, then we
140  // need to pick a threshold equal to the smallest indicator, but of course
141  // the bisection algorithm can never find a threshold equal to one of the
142  // end points of the interval. So we slightly increase the interval before
143  // we even start
144  void
145  adjust_interesting_range(double (&interesting_range)[2])
146  {
147  Assert(interesting_range[0] <= interesting_range[1], ExcInternalError());
148 
149  if (interesting_range[0] > 0)
150  {
151  // In this case, we calculate the first interval split point `m` in the
152  // `compute_threshold` functions in the optimized way: We exploit that
153  // the logarithms of all criteria are more uniformly distributed than
154  // their actual values, i.e. m=sqrt(b*e).
155  //
156  // Both factors will modify the split point only slightly by a factor of
157  // sqrt(1.01*0.99) = sqrt(0.9999) ~ 0.9950.
158  interesting_range[0] *= 0.99;
159  interesting_range[1] *= 1.01;
160  }
161  else
162  {
163  // In all other cases, we begin with an the arithmetic mean as the
164  // standard interval split point, i.e. m=(b+e)/2.
165  //
166  // Both increments will add up to zero when calculating the initial
167  // split point in the `compute_threshold` functions.
168  const double difference =
169  std::abs(interesting_range[1] - interesting_range[0]);
170  interesting_range[0] -= 0.01 * difference;
171  interesting_range[1] += 0.01 * difference;
172  }
173  }
174 
175 
176 
182  template <int dim, int spacedim, typename Number>
183  void
184  mark_cells(::Triangulation<dim, spacedim> &tria,
185  const ::Vector<Number> &criteria,
186  const double top_threshold,
187  const double bottom_threshold)
188  {
189  ::GridRefinement::refine(tria, criteria, top_threshold);
190  ::GridRefinement::coarsen(tria, criteria, bottom_threshold);
191 
192  // as a final good measure, delete all flags again from cells that we don't
193  // locally own
194  for (const auto &cell : tria.active_cell_iterators())
195  if (cell->subdomain_id() != tria.locally_owned_subdomain())
196  {
197  cell->clear_refine_flag();
198  cell->clear_coarsen_flag();
199  }
200  }
201 
202 
203 
211  template <int dim, int spacedim, typename Number>
212  void
213  refine_and_coarsen_fixed_fraction_via_l1_norm(
215  const ::Vector<Number> &criteria,
216  const double top_fraction_of_error,
217  const double bottom_fraction_of_error)
218  {
219  // first extract from the vector of indicators the ones that correspond
220  // to cells that we locally own
221  Vector<Number> locally_owned_indicators(n_locally_owned_active_cells(tria));
222  get_locally_owned_indicators(tria, criteria, locally_owned_indicators);
223 
224  MPI_Comm mpi_communicator = tria.get_communicator();
225 
226  // figure out the global max and min of the indicators. we don't need it
227  // here, but it's a collective communication call
228  const std::pair<double, double> global_min_and_max =
230  compute_global_min_and_max_at_root(locally_owned_indicators,
231  mpi_communicator);
232 
233  const double total_error =
234  compute_global_sum(locally_owned_indicators, mpi_communicator);
235 
236  double top_target_error = top_fraction_of_error * total_error,
237  bottom_target_error = (1. - bottom_fraction_of_error) * total_error;
238 
239  double top_threshold, bottom_threshold;
242  global_min_and_max,
243  top_target_error,
244  mpi_communicator);
245 
246  // compute bottom threshold only if necessary. otherwise use the lowest
247  // threshold possible
248  if (bottom_fraction_of_error > 0)
249  bottom_threshold = ::internal::parallel::distributed::
251  locally_owned_indicators,
252  global_min_and_max,
253  bottom_target_error,
254  mpi_communicator);
255  else
256  bottom_threshold = std::numeric_limits<Number>::lowest();
257 
258  // now refine the mesh
259  mark_cells(tria, criteria, top_threshold, bottom_threshold);
260  }
261 } // namespace
262 
263 
264 
265 namespace internal
266 {
267  namespace parallel
268  {
269  namespace distributed
270  {
271  namespace GridRefinement
272  {
273  template <typename number>
274  std::pair<number, number>
276  const ::Vector<number> &criteria,
277  const MPI_Comm mpi_communicator)
278  {
279  // we'd like to compute the global max and min from the local ones in
280  // one MPI communication. we can do that by taking the elementwise
281  // minimum of the local min and the negative maximum over all
282  // processors
283 
284  const double local_min = min_element(criteria),
285  local_max = max_element(criteria);
286  double comp[2] = {local_min, -local_max};
287  double result[2] = {0, 0};
288 
289  // compute the minimum on processor zero
290  const int ierr = MPI_Reduce(
291  comp, result, 2, MPI_DOUBLE, MPI_MIN, 0, mpi_communicator);
292  AssertThrowMPI(ierr);
293 
294  // make sure only processor zero got something
295  if (Utilities::MPI::this_mpi_process(mpi_communicator) != 0)
296  Assert((result[0] == 0) && (result[1] == 0), ExcInternalError());
297 
298  return std::make_pair(result[0], -result[1]);
299  }
300 
301 
302 
303  namespace RefineAndCoarsenFixedNumber
304  {
305  template <typename number>
306  number
307  compute_threshold(const ::Vector<number> &criteria,
308  const std::pair<double, double> &global_min_and_max,
309  const types::global_cell_index n_target_cells,
310  const MPI_Comm mpi_communicator)
311  {
312  double interesting_range[2] = {global_min_and_max.first,
313  global_min_and_max.second};
314  adjust_interesting_range(interesting_range);
315 
316  const unsigned int root_mpi_rank = 0;
317  unsigned int iteration = 0;
318 
319  do
320  {
321  int ierr = MPI_Bcast(interesting_range,
322  2,
323  MPI_DOUBLE,
324  root_mpi_rank,
325  mpi_communicator);
326  AssertThrowMPI(ierr);
327 
328  if (interesting_range[0] == interesting_range[1])
329  return interesting_range[0];
330 
331  const double test_threshold =
332  (interesting_range[0] > 0 ?
333  std::sqrt(interesting_range[0] * interesting_range[1]) :
334  (interesting_range[0] + interesting_range[1]) / 2);
335 
336  // Count how many of our own elements would be above this
337  // threshold. Use a 64bit result type if we are compiling with
338  // 64bit indices to avoid an overflow when computing the sum
339  // below.
340  const types::global_cell_index my_count =
341  std::count_if(criteria.begin(),
342  criteria.end(),
343  [test_threshold](const double c) {
344  return c > test_threshold;
345  });
346  const types::global_cell_index total_count =
347  Utilities::MPI::sum(my_count, mpi_communicator);
348 
349  // now adjust the range. if we have too many cells, we take the
350  // upper half of the previous range, otherwise the lower half.
351  // if we have hit the right number, then set the range to the
352  // exact value. non-root nodes also update their own
353  // interesting_range, however their results are not significant
354  // since the values will be overwritten by MPI_Bcast from the
355  // root node in next loop.
356  if (total_count > n_target_cells)
357  interesting_range[0] = test_threshold;
358  else if (total_count < n_target_cells)
359  interesting_range[1] = test_threshold;
360  else
361  interesting_range[0] = interesting_range[1] = test_threshold;
362 
363  // terminate the iteration after 25 go-arounds. this is
364  // necessary because oftentimes error indicators on cells have
365  // exactly the same value, and so there may not be a particular
366  // value that cuts the indicators in such a way that we can
367  // achieve the desired number of cells. using a maximal number
368  // of iterations means that we terminate the iteration after a
369  // fixed number N of steps if the indicators were perfectly
370  // badly distributed, and we make at most a mistake of 1/2^N in
371  // the number of cells flagged if indicators are perfectly
372  // equidistributed
373  ++iteration;
374  if (iteration == 25)
375  interesting_range[0] = interesting_range[1] = test_threshold;
376  }
377  while (true);
378 
379  Assert(false, ExcInternalError());
380  return -1;
381  }
382  } // namespace RefineAndCoarsenFixedNumber
383 
384 
385 
386  namespace RefineAndCoarsenFixedFraction
387  {
388  template <typename number>
389  number
390  compute_threshold(const ::Vector<number> &criteria,
391  const std::pair<double, double> &global_min_and_max,
392  const double target_error,
393  const MPI_Comm mpi_communicator)
394  {
395  double interesting_range[2] = {global_min_and_max.first,
396  global_min_and_max.second};
397  adjust_interesting_range(interesting_range);
398 
399  const unsigned int root_mpi_rank = 0;
400  unsigned int iteration = 0;
401 
402  do
403  {
404  int ierr = MPI_Bcast(interesting_range,
405  2,
406  MPI_DOUBLE,
407  root_mpi_rank,
408  mpi_communicator);
409  AssertThrowMPI(ierr);
410 
411  if (interesting_range[0] == interesting_range[1])
412  {
413  // so we have found our threshold. since we adjust the range
414  // at the top of the function to be slightly larger than the
415  // actual extremes of the refinement criteria values, we can
416  // end up in a situation where the threshold is in fact
417  // larger than the maximal refinement indicator. in such
418  // cases, we get no refinement at all. thus, cap the
419  // threshold by the actual largest value
420  double final_threshold =
421  std::min(interesting_range[0], global_min_and_max.second);
422  ierr = MPI_Bcast(&final_threshold,
423  1,
424  MPI_DOUBLE,
425  root_mpi_rank,
426  mpi_communicator);
427  AssertThrowMPI(ierr);
428 
429  return final_threshold;
430  }
431 
432  const double test_threshold =
433  (interesting_range[0] > 0 ?
434  std::sqrt(interesting_range[0] * interesting_range[1]) :
435  (interesting_range[0] + interesting_range[1]) / 2);
436 
437  // accumulate the error of those our own elements above this
438  // threshold and then add to it the number for all the others
439  double my_error = 0;
440  for (unsigned int i = 0; i < criteria.size(); ++i)
441  if (criteria(i) > test_threshold)
442  my_error += criteria(i);
443 
444  double total_error = 0.;
445 
446  ierr = MPI_Reduce(&my_error,
447  &total_error,
448  1,
449  MPI_DOUBLE,
450  MPI_SUM,
451  root_mpi_rank,
452  mpi_communicator);
453  AssertThrowMPI(ierr);
454 
455  // now adjust the range. if we have too many cells, we take the
456  // upper half of the previous range, otherwise the lower half.
457  // if we have hit the right number, then set the range to the
458  // exact value. non-root nodes also update their own
459  // interesting_range, however their results are not significant
460  // since the values will be overwritten by MPI_Bcast from the
461  // root node in next loop.
462  if (total_error > target_error)
463  interesting_range[0] = test_threshold;
464  else if (total_error < target_error)
465  interesting_range[1] = test_threshold;
466  else
467  interesting_range[0] = interesting_range[1] = test_threshold;
468 
469  // terminate the iteration after 25 go-arounds. this is
470  // necessary because oftentimes error indicators on cells
471  // have exactly the same value, and so there may not be a
472  // particular value that cuts the indicators in such a way
473  // that we can achieve the desired number of cells. using a
474  // max of 25 iterations means that we terminate the
475  // iteration after 25 steps if the indicators were perfectly
476  // badly distributed, and we make at most a mistake of
477  // 1/2^25 in the number of cells flagged if indicators are
478  // perfectly equidistributed
479  ++iteration;
480  if (iteration == 25)
481  interesting_range[0] = interesting_range[1] = test_threshold;
482  }
483  while (true);
484 
485  Assert(false, ExcInternalError());
486  return -1;
487  }
488  } // namespace RefineAndCoarsenFixedFraction
489  } // namespace GridRefinement
490  } // namespace distributed
491  } // namespace parallel
492 } // namespace internal
493 
494 
495 
496 namespace parallel
497 {
498  namespace distributed
499  {
500  namespace GridRefinement
501  {
502  template <int dim, typename Number, int spacedim>
503  void
506  const ::Vector<Number> &criteria,
507  const double top_fraction_of_cells,
508  const double bottom_fraction_of_cells,
509  const types::global_cell_index max_n_cells)
510  {
511  Assert(criteria.size() == tria.n_active_cells(),
512  ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
513  Assert((top_fraction_of_cells >= 0) && (top_fraction_of_cells <= 1),
515  Assert((bottom_fraction_of_cells >= 0) &&
516  (bottom_fraction_of_cells <= 1),
518  Assert(top_fraction_of_cells + bottom_fraction_of_cells <= 1,
520  Assert(criteria.is_non_negative(),
522 
523  const std::pair<double, double> adjusted_fractions =
526  max_n_cells,
527  top_fraction_of_cells,
528  bottom_fraction_of_cells);
529 
530  // first extract from the vector of indicators the ones that correspond
531  // to cells that we locally own
532  Vector<Number> locally_owned_indicators(
533  n_locally_owned_active_cells(tria));
534  get_locally_owned_indicators(tria, criteria, locally_owned_indicators);
535 
536  MPI_Comm mpi_communicator = tria.get_communicator();
537 
538  // figure out the global max and min of the indicators. we don't need it
539  // here, but it's a collective communication call
540  const std::pair<Number, Number> global_min_and_max =
542  compute_global_min_and_max_at_root(locally_owned_indicators,
543  mpi_communicator);
544 
545 
546  double top_threshold, bottom_threshold;
549  locally_owned_indicators,
550  global_min_and_max,
551  static_cast<types::global_cell_index>(adjusted_fractions.first *
553  mpi_communicator);
554 
555  // compute bottom threshold only if necessary. otherwise use the lowest
556  // threshold possible
557  if (adjusted_fractions.second > 0)
558  bottom_threshold = ::internal::parallel::distributed::
560  locally_owned_indicators,
561  global_min_and_max,
562  static_cast<types::global_cell_index>(
563  std::ceil((1. - adjusted_fractions.second) *
565  mpi_communicator);
566  else
567  bottom_threshold = std::numeric_limits<Number>::lowest();
568 
569  // now refine the mesh
570  mark_cells(tria, criteria, top_threshold, bottom_threshold);
571  }
572 
573 
574 
575  template <int dim, typename Number, int spacedim>
576  void
579  const ::Vector<Number> &criteria,
580  const double top_fraction_of_error,
581  const double bottom_fraction_of_error,
582  const VectorTools::NormType norm_type)
583  {
584  Assert(criteria.size() == tria.n_active_cells(),
585  ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
586  Assert((top_fraction_of_error >= 0) && (top_fraction_of_error <= 1),
588  Assert((bottom_fraction_of_error >= 0) &&
589  (bottom_fraction_of_error <= 1),
591  Assert(top_fraction_of_error + bottom_fraction_of_error <= 1,
593  Assert(criteria.is_non_negative(),
595 
596  switch (norm_type)
597  {
599  // evaluate norms on subsets and compare them as
600  // c_0 + c_1 + ... < fraction * l1-norm(c)
601  refine_and_coarsen_fixed_fraction_via_l1_norm(
602  tria,
603  criteria,
604  top_fraction_of_error,
605  bottom_fraction_of_error);
606  break;
607 
609  {
610  // we do not want to evaluate norms on subsets as:
611  // sqrt(c_0^2 + c_1^2 + ...) < fraction * l2-norm(c)
612  // instead take the square of both sides of the equation
613  // and evaluate:
614  // c_0^2 + c_1^2 + ... < fraction^2 * l1-norm(c.c)
615  // we adjust all parameters accordingly
616  Vector<Number> criteria_squared(criteria.size());
617  std::transform(criteria.begin(),
618  criteria.end(),
619  criteria_squared.begin(),
620  [](Number c) { return c * c; });
621 
622  refine_and_coarsen_fixed_fraction_via_l1_norm(
623  tria,
624  criteria_squared,
625  top_fraction_of_error * top_fraction_of_error,
626  bottom_fraction_of_error * bottom_fraction_of_error);
627  }
628  break;
629 
630  default:
631  Assert(false, ExcNotImplemented());
632  break;
633  }
634  }
635  } // namespace GridRefinement
636  } // namespace distributed
637 } // namespace parallel
638 
639 
640 // explicit instantiations
641 # include "grid_refinement.inst"
642 
644 
645 #endif
virtual types::global_cell_index n_global_active_cells() const
virtual MPI_Comm get_communicator() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_cells() const
Definition: vector.h:110
virtual size_type size() const override
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1916
static ::ExceptionBase & ExcNegativeCriteria()
static ::ExceptionBase & ExcInvalidParameterValue()
Expression ceil(const Expression &x)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
std::pair< double, double > adjust_refine_and_coarsen_number_fraction(const types::global_cell_index current_n_cells, const types::global_cell_index max_n_cells, const double top_fraction_of_cells, const double bottom_fraction_of_cells)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition: mpi.cc:164
number compute_threshold(const ::Vector< number > &criteria, const std::pair< double, double > &global_min_and_max, const double target_error, const MPI_Comm mpi_communicator)
number compute_threshold(const ::Vector< number > &criteria, const std::pair< double, double > &global_min_and_max, const types::global_cell_index n_target_cells, const MPI_Comm mpi_communicator)
std::pair< number, number > compute_global_min_and_max_at_root(const ::Vector< number > &criteria, const MPI_Comm mpi_communicator)
void refine_and_coarsen_fixed_fraction(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_error, const double bottom_fraction_of_error, const VectorTools::NormType norm_type=VectorTools::NormType::L1_norm)
void refine_and_coarsen_fixed_number(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const types::global_cell_index max_n_cells=std::numeric_limits< types::global_cell_index >::max())
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Function &function, const unsigned int grainsize)
Definition: parallel.h:166
unsigned int global_cell_index
Definition: types.h:118
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::Triangulation< dim, spacedim > & tria