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