Reference documentation for deal.II version GIT c00406db16 2023-09-28 18: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\}}\)
grid_refinement.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2022 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 
18 
20 
22 #include <deal.II/grid/tria.h>
25 
26 #include <deal.II/lac/vector.h>
27 
28 #include <algorithm>
29 #include <cmath>
30 #include <fstream>
31 #include <functional>
32 #include <limits>
33 #include <numeric>
34 
36 
37 namespace
38 {
46  template <int dim, int spacedim, typename Number>
47  void
48  refine_and_coarsen_fixed_fraction_via_l1_norm(
50  const Vector<Number> &criteria,
51  const double top_fraction,
52  const double bottom_fraction,
53  const unsigned int max_n_cells)
54  {
55  // sort the criteria in descending order in an auxiliary vector, which we
56  // have to sum up and compare with @p{fraction_of_error*total_error}
57  Vector<Number> criteria_sorted = criteria;
58  std::sort(criteria_sorted.begin(),
59  criteria_sorted.end(),
60  std::greater<double>());
61 
62  const double total_error = criteria_sorted.l1_norm();
63 
64  // compute thresholds
65  typename Vector<Number>::const_iterator pp = criteria_sorted.begin();
66  for (double sum = 0;
67  (sum < top_fraction * total_error) && (pp != criteria_sorted.end());
68  ++pp)
69  sum += *pp;
70  double top_threshold =
71  (pp != criteria_sorted.begin() ? (*pp + *(pp - 1)) / 2 : *pp);
72 
73  typename Vector<Number>::const_iterator qq = criteria_sorted.end() - 1;
74  for (double sum = 0; (sum < bottom_fraction * total_error) &&
75  (qq != criteria_sorted.begin() - 1);
76  --qq)
77  sum += *qq;
78  double bottom_threshold =
79  ((qq != criteria_sorted.end() - 1) ? (*qq + *(qq + 1)) / 2 : 0.);
80 
81  // we now have an idea how many cells we
82  // are going to refine and coarsen. we use
83  // this information to see whether we are
84  // over the limit and if so use a function
85  // that knows how to deal with this
86  // situation
87 
88  // note, that at this point, we have no
89  // information about anisotropically refined
90  // cells, thus use the situation of purely
91  // isotropic refinement as guess for a mixed
92  // refinemnt as well.
93  const unsigned int refine_cells = pp - criteria_sorted.begin(),
94  coarsen_cells = criteria_sorted.end() - 1 - qq;
95 
96  if (static_cast<unsigned int>(
98  refine_cells * (GeometryInfo<dim>::max_children_per_cell - 1) -
99  (coarsen_cells * (GeometryInfo<dim>::max_children_per_cell - 1) /
101  {
103  criteria,
104  1. * refine_cells /
105  criteria.size(),
106  1. * coarsen_cells /
107  criteria.size(),
108  max_n_cells);
109  return;
110  }
111 
112  // in some rare cases it may happen that
113  // both thresholds are the same (e.g. if
114  // there are many cells with the same
115  // error indicator). That would mean that
116  // all cells will be flagged for
117  // refinement or coarsening, but some will
118  // be flagged for both, namely those for
119  // which the indicator equals the
120  // thresholds. This is forbidden, however.
121  //
122  // In some rare cases with very few cells
123  // we also could get integer round off
124  // errors and get problems with
125  // the top and bottom fractions.
126  //
127  // In these case we arbitrarily reduce the
128  // bottom threshold by one permille below
129  // the top threshold
130  //
131  // Finally, in some cases
132  // (especially involving symmetric
133  // solutions) there are many cells
134  // with the same error indicator
135  // values. if there are many with
136  // indicator equal to the top
137  // threshold, no refinement will
138  // take place below; to avoid this
139  // case, we also lower the top
140  // threshold if it equals the
141  // largest indicator and the
142  // top_fraction!=1
143  const double max_criterion = *(criteria_sorted.begin()),
144  min_criterion = *(criteria_sorted.end() - 1);
145 
146  if ((top_threshold == max_criterion) && (top_fraction != 1))
147  top_threshold *= 0.999;
148 
149  if (bottom_threshold >= top_threshold)
150  bottom_threshold = 0.999 * top_threshold;
151 
152  // actually flag cells
153  if (top_threshold < max_criterion)
154  GridRefinement::refine(tria, criteria, top_threshold, refine_cells);
155 
156  if (bottom_threshold > min_criterion)
157  GridRefinement::coarsen(tria, criteria, bottom_threshold);
158  }
159 } // namespace
160 
161 
162 
163 template <int dim, typename Number, int spacedim>
164 void
166  const Vector<Number> &criteria,
167  const double threshold,
168  const unsigned int max_to_mark)
169 {
170  Assert(criteria.size() == tria.n_active_cells(),
173 
174  // when all indicators are zero we
175  // do not need to refine but only
176  // to coarsen
177  if (criteria.all_zero())
178  return;
179 
180  const unsigned int n_cells = criteria.size();
181 
182  // TODO: This is undocumented, looks fishy and seems unnecessary
183 
184  double new_threshold = threshold;
185  // when threshold==0 find the
186  // smallest value in criteria
187  // greater 0
188  if (new_threshold == 0)
189  {
190  new_threshold = criteria(0);
191  for (unsigned int index = 1; index < n_cells; ++index)
192  if (criteria(index) > 0 && (criteria(index) < new_threshold))
193  new_threshold = criteria(index);
194  }
195 
196  unsigned int marked = 0;
197  for (const auto &cell : tria.active_cell_iterators())
199  &tria) == nullptr ||
200  cell->is_locally_owned()) &&
201  std::fabs(criteria(cell->active_cell_index())) >= new_threshold)
202  {
203  if (max_to_mark != numbers::invalid_unsigned_int &&
204  marked >= max_to_mark)
205  break;
206  ++marked;
207  cell->set_refine_flag();
208  }
209 }
210 
211 
212 
213 template <int dim, typename Number, int spacedim>
214 void
216  const Vector<Number> &criteria,
217  const double threshold)
218 {
219  Assert(criteria.size() == tria.n_active_cells(),
222 
223  for (const auto &cell : tria.active_cell_iterators())
225  &tria) == nullptr ||
226  cell->is_locally_owned()) &&
227  std::fabs(criteria(cell->active_cell_index())) <= threshold)
228  if (!cell->refine_flag_set())
229  cell->set_coarsen_flag();
230 }
231 
232 
233 
234 template <int dim>
235 std::pair<double, double>
237  const types::global_cell_index current_n_cells,
238  const types::global_cell_index max_n_cells,
239  const double top_fraction,
240  const double bottom_fraction)
241 {
242  Assert(top_fraction >= 0, ExcInvalidParameterValue());
243  Assert(top_fraction <= 1, ExcInvalidParameterValue());
244  Assert(bottom_fraction >= 0, ExcInvalidParameterValue());
245  Assert(bottom_fraction <= 1, ExcInvalidParameterValue());
246  Assert(top_fraction + bottom_fraction <=
249 
250  double refine_cells = current_n_cells * top_fraction;
251  double coarsen_cells = current_n_cells * bottom_fraction;
252 
253  const double cell_increase_on_refine =
255  const double cell_decrease_on_coarsen =
257 
258  std::pair<double, double> adjusted_fractions(top_fraction, bottom_fraction);
259  // first we have to see whether we
260  // currently already exceed the target
261  // number of cells
262  if (current_n_cells >= max_n_cells)
263  {
264  // if yes, then we need to stop
265  // refining cells and instead try to
266  // only coarsen as many as it would
267  // take to get to the target
268 
269  // as we have no information on cells
270  // being refined isotropically or
271  // anisotropically, assume isotropic
272  // refinement here, though that may
273  // result in a worse approximation
274  adjusted_fractions.first = 0;
275  coarsen_cells =
276  (current_n_cells - max_n_cells) / cell_decrease_on_coarsen;
277  adjusted_fractions.second =
278  std::min(coarsen_cells / current_n_cells, 1.0);
279  }
280  // otherwise, see if we would exceed the
281  // maximum desired number of cells with the
282  // number of cells that are likely going to
283  // result from refinement. here, each cell
284  // to be refined is replaced by
285  // C=GeometryInfo<dim>::max_children_per_cell
286  // new cells, i.e. there will be C-1 more
287  // cells than before. similarly, C cells
288  // will be replaced by 1
289 
290  // again, this is true for isotropically
291  // refined cells. we take this as an
292  // approximation of a mixed refinement.
293  else if (static_cast<types::global_cell_index>(
294  current_n_cells + refine_cells * cell_increase_on_refine -
295  coarsen_cells * cell_decrease_on_coarsen) > max_n_cells)
296  {
297  // we have to adjust the
298  // fractions. assume we want
299  // alpha*refine_fraction and
300  // alpha*coarsen_fraction as new
301  // fractions and the resulting number
302  // of cells to be equal to
303  // max_n_cells. this leads to the
304  // following equation for alpha
305  const double alpha = 1. * (max_n_cells - current_n_cells) /
306  (refine_cells * cell_increase_on_refine -
307  coarsen_cells * cell_decrease_on_coarsen);
308 
309  adjusted_fractions.first = alpha * top_fraction;
310  adjusted_fractions.second = alpha * bottom_fraction;
311  }
312  return (adjusted_fractions);
313 }
314 
315 
316 
317 template <int dim, typename Number, int spacedim>
318 void
321  const Vector<Number> &criteria,
322  const double top_fraction,
323  const double bottom_fraction,
324  const unsigned int max_n_cells)
325 {
326  // correct number of cells is
327  // checked in @p{refine}
328  Assert((top_fraction >= 0) && (top_fraction <= 1),
330  Assert((bottom_fraction >= 0) && (bottom_fraction <= 1),
332  Assert(top_fraction + bottom_fraction <=
336 
337  const std::pair<double, double> adjusted_fractions =
338  adjust_refine_and_coarsen_number_fraction<dim>(criteria.size(),
339  max_n_cells,
340  top_fraction,
341  bottom_fraction);
342 
343  const int refine_cells =
344  static_cast<int>(adjusted_fractions.first * criteria.size());
345  const int coarsen_cells =
346  static_cast<int>(adjusted_fractions.second * criteria.size());
347 
348  if (refine_cells || coarsen_cells)
349  {
350  Vector<Number> tmp(criteria);
351  if (refine_cells)
352  {
353  if (static_cast<size_t>(refine_cells) == criteria.size())
354  refine(tria, criteria, std::numeric_limits<double>::lowest());
355  else
356  {
357  std::nth_element(tmp.begin(),
358  tmp.begin() + refine_cells - 1,
359  tmp.end(),
360  std::greater<double>());
361  refine(tria, criteria, *(tmp.begin() + refine_cells - 1));
362  }
363  }
364 
365  if (coarsen_cells)
366  {
367  if (static_cast<size_t>(coarsen_cells) == criteria.size())
369  else
370  {
371  std::nth_element(tmp.begin(),
372  tmp.begin() + tmp.size() - coarsen_cells,
373  tmp.end(),
374  std::greater<double>());
375  coarsen(tria,
376  criteria,
377  *(tmp.begin() + tmp.size() - coarsen_cells));
378  }
379  }
380  }
381 }
382 
383 
384 
385 template <int dim, typename Number, int spacedim>
386 void
389  const Vector<Number> &criteria,
390  const double top_fraction,
391  const double bottom_fraction,
392  const unsigned int max_n_cells,
393  const VectorTools::NormType norm_type)
394 {
395  // correct number of cells is checked in @p{refine}
396  Assert((top_fraction >= 0) && (top_fraction <= 1),
398  Assert((bottom_fraction >= 0) && (bottom_fraction <= 1),
400  Assert(top_fraction + bottom_fraction <=
404 
405  switch (norm_type)
406  {
408  // evaluate norms on subsets and compare them as
409  // c_0 + c_1 + ... < fraction * l1-norm(c)
410  refine_and_coarsen_fixed_fraction_via_l1_norm(
411  tria, criteria, top_fraction, bottom_fraction, max_n_cells);
412  break;
413 
415  {
416  // we do not want to evaluate norms on subsets as:
417  // sqrt(c_0^2 + c_1^2 + ...) < fraction * l2-norm(c)
418  // instead take the square of both sides of the equation
419  // and evaluate:
420  // c_0^2 + c_1^2 + ... < fraction^2 * l1-norm(c.c)
421  // we adjust all parameters accordingly
422  Vector<Number> criteria_squared(criteria.size());
423  std::transform(criteria.begin(),
424  criteria.end(),
425  criteria_squared.begin(),
426  [](Number c) { return c * c; });
427 
428  refine_and_coarsen_fixed_fraction_via_l1_norm(tria,
429  criteria_squared,
430  top_fraction *
431  top_fraction,
432  bottom_fraction *
433  bottom_fraction,
434  max_n_cells);
435  }
436  break;
437 
438  default:
439  Assert(false, ExcNotImplemented());
440  break;
441  }
442 }
443 
444 
445 
446 template <int dim, typename Number, int spacedim>
447 void
449  const Vector<Number> &criteria,
450  const unsigned int order)
451 {
452  Assert(criteria.size() == tria.n_active_cells(),
455 
456  // get a decreasing order on the error indicator
457  std::vector<unsigned int> cell_indices(criteria.size());
458  std::iota(cell_indices.begin(), cell_indices.end(), 0u);
459 
460  std::sort(cell_indices.begin(),
461  cell_indices.end(),
462  [&criteria](const unsigned int left, const unsigned int right) {
463  return criteria[left] > criteria[right];
464  });
465 
466  double expected_error_reduction = 0;
467  const double original_error = criteria.l1_norm();
468 
469  const std::size_t N = criteria.size();
470 
471  // minimize the cost functional discussed in the documentation
472  double min_cost = std::numeric_limits<double>::max();
473  std::size_t min_arg = 0;
474 
475  for (std::size_t M = 0; M < criteria.size(); ++M)
476  {
477  expected_error_reduction +=
478  (1 - std::pow(2., -1. * order)) * criteria(cell_indices[M]);
479 
480  const double cost = std::pow(((std::pow(2., dim) - 1) * (1 + M) + N),
481  static_cast<double>(order) / dim) *
482  (original_error - expected_error_reduction);
483  if (cost <= min_cost)
484  {
485  min_cost = cost;
486  min_arg = M;
487  }
488  }
489 
490  refine(tria, criteria, criteria(cell_indices[min_arg]));
491 }
492 
493 
494 // explicit instantiations
495 #include "grid_refinement.inst"
496 
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
unsigned int n_active_cells() const
Definition: vector.h:110
const value_type * const_iterator
Definition: vector.h:133
virtual size_type size() const override
iterator end()
bool is_non_negative() const
bool all_zero() const
real_type l1_norm() const
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 & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNegativeCriteria()
static ::ExceptionBase & ExcInvalidParameterValue()
Expression fabs(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 refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
void refine_and_coarsen_optimize(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const unsigned int order=2)
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)
void refine_and_coarsen_fixed_fraction(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double top_fraction, const double bottom_fraction, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max(), const VectorTools::NormType norm_type=VectorTools::NormType::L1_norm)
static const char N
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:14833
static const unsigned int invalid_unsigned_int
Definition: types.h:213
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
const ::Triangulation< dim, spacedim > & tria