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