Reference documentation for deal.II version Git e7bb9ce7b3 2020-09-18 12:07:32 -0400
\(\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 
41 template <int dim, typename Number, int spacedim>
42 void
44  const Vector<Number> & criteria,
45  const double threshold,
46  const unsigned int max_to_mark)
47 {
48  Assert(criteria.size() == tria.n_active_cells(),
49  ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
50  Assert(criteria.is_non_negative(), ExcNegativeCriteria());
51 
52  // when all indicators are zero we
53  // do not need to refine but only
54  // to coarsen
55  if (criteria.all_zero())
56  return;
57 
58  const unsigned int n_cells = criteria.size();
59 
60  // TODO: This is undocumented, looks fishy and seems unnecessary
61 
62  double new_threshold = threshold;
63  // when threshold==0 find the
64  // smallest value in criteria
65  // greater 0
66  if (new_threshold == 0)
67  {
68  new_threshold = criteria(0);
69  for (unsigned int index = 1; index < n_cells; ++index)
70  if (criteria(index) > 0 && (criteria(index) < new_threshold))
71  new_threshold = criteria(index);
72  }
73 
74  unsigned int marked = 0;
75  for (const auto &cell : tria.active_cell_iterators())
77  &tria) == nullptr ||
78  cell->is_locally_owned()) &&
79  std::fabs(criteria(cell->active_cell_index())) >= new_threshold)
80  {
81  if (max_to_mark != numbers::invalid_unsigned_int &&
82  marked >= max_to_mark)
83  break;
84  ++marked;
85  cell->set_refine_flag();
86  }
87 }
88 
89 
90 
91 template <int dim, typename Number, int spacedim>
92 void
94  const Vector<Number> & criteria,
95  const double threshold)
96 {
97  Assert(criteria.size() == tria.n_active_cells(),
98  ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
99  Assert(criteria.is_non_negative(), ExcNegativeCriteria());
100 
101  for (const auto &cell : tria.active_cell_iterators())
103  &tria) == nullptr ||
104  cell->is_locally_owned()) &&
105  std::fabs(criteria(cell->active_cell_index())) <= threshold)
106  if (!cell->refine_flag_set())
107  cell->set_coarsen_flag();
108 }
109 
110 
111 
112 template <int dim>
113 std::pair<double, double>
115  const types::global_cell_index current_n_cells,
116  const types::global_cell_index max_n_cells,
117  const double top_fraction,
118  const double bottom_fraction)
119 {
120  Assert(top_fraction >= 0, ExcInvalidParameterValue());
121  Assert(top_fraction <= 1, ExcInvalidParameterValue());
122  Assert(bottom_fraction >= 0, ExcInvalidParameterValue());
123  Assert(bottom_fraction <= 1, ExcInvalidParameterValue());
124  Assert(top_fraction + bottom_fraction <=
127 
128  double refine_cells = current_n_cells * top_fraction;
129  double coarsen_cells = current_n_cells * bottom_fraction;
130 
131  const double cell_increase_on_refine =
133  const double cell_decrease_on_coarsen =
135 
136  std::pair<double, double> adjusted_fractions(top_fraction, bottom_fraction);
137  // first we have to see whether we
138  // currently already exceed the target
139  // number of cells
140  if (current_n_cells >= max_n_cells)
141  {
142  // if yes, then we need to stop
143  // refining cells and instead try to
144  // only coarsen as many as it would
145  // take to get to the target
146 
147  // as we have no information on cells
148  // being refined isotropically or
149  // anisotropically, assume isotropic
150  // refinement here, though that may
151  // result in a worse approximation
152  adjusted_fractions.first = 0;
153  coarsen_cells =
154  (current_n_cells - max_n_cells) / cell_decrease_on_coarsen;
155  adjusted_fractions.second =
156  std::min(coarsen_cells / current_n_cells, 1.0);
157  }
158  // otherwise, see if we would exceed the
159  // maximum desired number of cells with the
160  // number of cells that are likely going to
161  // result from refinement. here, each cell
162  // to be refined is replaced by
163  // C=GeometryInfo<dim>::max_children_per_cell
164  // new cells, i.e. there will be C-1 more
165  // cells than before. similarly, C cells
166  // will be replaced by 1
167 
168  // again, this is true for isotropically
169  // refined cells. we take this as an
170  // approximation of a mixed refinement.
171  else if (static_cast<types::global_cell_index>(
172  current_n_cells + refine_cells * cell_increase_on_refine -
173  coarsen_cells * cell_decrease_on_coarsen) > max_n_cells)
174  {
175  // we have to adjust the
176  // fractions. assume we want
177  // alpha*refine_fraction and
178  // alpha*coarsen_fraction as new
179  // fractions and the resulting number
180  // of cells to be equal to
181  // max_n_cells. this leads to the
182  // following equation for alpha
183  const double alpha = 1. * (max_n_cells - current_n_cells) /
184  (refine_cells * cell_increase_on_refine -
185  coarsen_cells * cell_decrease_on_coarsen);
186 
187  adjusted_fractions.first = alpha * top_fraction;
188  adjusted_fractions.second = alpha * bottom_fraction;
189  }
190  return (adjusted_fractions);
191 }
192 
193 
194 
195 template <int dim, typename Number, int spacedim>
196 void
199  const Vector<Number> & criteria,
200  const double top_fraction,
201  const double bottom_fraction,
202  const unsigned int max_n_cells)
203 {
204  // correct number of cells is
205  // checked in @p{refine}
206  Assert((top_fraction >= 0) && (top_fraction <= 1),
208  Assert((bottom_fraction >= 0) && (bottom_fraction <= 1),
210  Assert(top_fraction + bottom_fraction <=
213  Assert(criteria.is_non_negative(), ExcNegativeCriteria());
214 
215  const std::pair<double, double> adjusted_fractions =
216  adjust_refine_and_coarsen_number_fraction<dim>(criteria.size(),
217  max_n_cells,
218  top_fraction,
219  bottom_fraction);
220 
221  const int refine_cells =
222  static_cast<int>(adjusted_fractions.first * criteria.size());
223  const int coarsen_cells =
224  static_cast<int>(adjusted_fractions.second * criteria.size());
225 
226  if (refine_cells || coarsen_cells)
227  {
228  Vector<Number> tmp(criteria);
229  if (refine_cells)
230  {
231  if (static_cast<size_t>(refine_cells) == criteria.size())
232  refine(tria, criteria, -std::numeric_limits<double>::max());
233  else
234  {
235  std::nth_element(tmp.begin(),
236  tmp.begin() + refine_cells - 1,
237  tmp.end(),
238  std::greater<double>());
239  refine(tria, criteria, *(tmp.begin() + refine_cells - 1));
240  }
241  }
242 
243  if (coarsen_cells)
244  {
245  if (static_cast<size_t>(coarsen_cells) == criteria.size())
246  coarsen(tria, criteria, std::numeric_limits<double>::max());
247  else
248  {
249  std::nth_element(tmp.begin(),
250  tmp.begin() + tmp.size() - coarsen_cells,
251  tmp.end(),
252  std::greater<double>());
253  coarsen(tria,
254  criteria,
255  *(tmp.begin() + tmp.size() - coarsen_cells));
256  }
257  }
258  }
259 }
260 
261 
262 
263 template <int dim, typename Number, int spacedim>
264 void
267  const Vector<Number> & criteria,
268  const double top_fraction,
269  const double bottom_fraction,
270  const unsigned int max_n_cells)
271 {
272  // correct number of cells is
273  // checked in @p{refine}
274  Assert((top_fraction >= 0) && (top_fraction <= 1),
276  Assert((bottom_fraction >= 0) && (bottom_fraction <= 1),
278  Assert(top_fraction + bottom_fraction <=
281  Assert(criteria.is_non_negative(), ExcNegativeCriteria());
282 
283  // let tmp be the cellwise square of the
284  // error, which is what we have to sum
285  // up and compare with
286  // @p{fraction_of_error*total_error}.
287  Vector<Number> tmp;
288  tmp = criteria;
289  const double total_error = tmp.l1_norm();
290 
291  // sort the largest criteria to the
292  // beginning of the vector
293  std::sort(tmp.begin(), tmp.end(), std::greater<double>());
294 
295  // compute thresholds
296  typename Vector<Number>::const_iterator pp = tmp.begin();
297  for (double sum = 0; (sum < top_fraction * total_error) && (pp != tmp.end());
298  ++pp)
299  sum += *pp;
300  double top_threshold = (pp != tmp.begin() ? (*pp + *(pp - 1)) / 2 : *pp);
301  typename Vector<Number>::const_iterator qq = (tmp.end() - 1);
302  for (double sum = 0;
303  (sum < bottom_fraction * total_error) && (qq != tmp.begin() - 1);
304  --qq)
305  sum += *qq;
306  double bottom_threshold = (qq != (tmp.end() - 1) ? (*qq + *(qq + 1)) / 2 : 0);
307 
308  // we now have an idea how many cells we
309  // are going to refine and coarsen. we use
310  // this information to see whether we are
311  // over the limit and if so use a function
312  // that knows how to deal with this
313  // situation
314 
315  // note, that at this point, we have no
316  // information about anisotropically refined
317  // cells, thus use the situation of purely
318  // isotropic refinement as guess for a mixed
319  // refinemnt as well.
320  {
321  const unsigned int refine_cells = pp - tmp.begin(),
322  coarsen_cells = tmp.end() - 1 - qq;
323 
324  if (static_cast<unsigned int>(
325  tria.n_active_cells() +
326  refine_cells * (GeometryInfo<dim>::max_children_per_cell - 1) -
327  (coarsen_cells * (GeometryInfo<dim>::max_children_per_cell - 1) /
329  {
331  criteria,
332  1. * refine_cells / criteria.size(),
333  1. * coarsen_cells / criteria.size(),
334  max_n_cells);
335  return;
336  }
337  }
338 
339 
340  // in some rare cases it may happen that
341  // both thresholds are the same (e.g. if
342  // there are many cells with the same
343  // error indicator). That would mean that
344  // all cells will be flagged for
345  // refinement or coarsening, but some will
346  // be flagged for both, namely those for
347  // which the indicator equals the
348  // thresholds. This is forbidden, however.
349  //
350  // In some rare cases with very few cells
351  // we also could get integer round off
352  // errors and get problems with
353  // the top and bottom fractions.
354  //
355  // In these case we arbitrarily reduce the
356  // bottom threshold by one permille below
357  // the top threshold
358  //
359  // Finally, in some cases
360  // (especially involving symmetric
361  // solutions) there are many cells
362  // with the same error indicator
363  // values. if there are many with
364  // indicator equal to the top
365  // threshold, no refinement will
366  // take place below; to avoid this
367  // case, we also lower the top
368  // threshold if it equals the
369  // largest indicator and the
370  // top_fraction!=1
371  const auto minmax_element =
372  std::minmax_element(criteria.begin(), criteria.end());
373  if ((top_threshold == *minmax_element.second) && (top_fraction != 1))
374  top_threshold *= 0.999;
375 
376  if (bottom_threshold >= top_threshold)
377  bottom_threshold = 0.999 * top_threshold;
378 
379  // actually flag cells
380  if (top_threshold < *minmax_element.second)
381  refine(tria, criteria, top_threshold, pp - tmp.begin());
382 
383  if (bottom_threshold > *minmax_element.first)
384  coarsen(tria, criteria, bottom_threshold);
385 }
386 
387 
388 
389 template <int dim, typename Number, int spacedim>
390 void
392  const Vector<Number> &criteria,
393  const unsigned int order)
394 {
395  Assert(criteria.size() == tria.n_active_cells(),
396  ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
397  Assert(criteria.is_non_negative(), ExcNegativeCriteria());
398 
399  // get a decreasing order on the error indicator
400  std::vector<unsigned int> cell_indices(criteria.size());
401  std::iota(cell_indices.begin(), cell_indices.end(), 0u);
402 
403  std::sort(cell_indices.begin(),
404  cell_indices.end(),
405  [&criteria](const unsigned int left, const unsigned int right) {
406  return criteria[left] > criteria[right];
407  });
408 
409  double expected_error_reduction = 0;
410  const double original_error = criteria.l1_norm();
411 
412  const std::size_t N = criteria.size();
413 
414  // minimize the cost functional discussed in the documentation
415  double min_cost = std::numeric_limits<double>::max();
416  std::size_t min_arg = 0;
417 
418  for (std::size_t M = 0; M < criteria.size(); ++M)
419  {
420  expected_error_reduction +=
421  (1 - std::pow(2., -1. * order)) * criteria(cell_indices[M]);
422 
423  const double cost = std::pow(((std::pow(2., dim) - 1) * (1 + M) + N),
424  static_cast<double>(order) / dim) *
425  (original_error - expected_error_reduction);
426  if (cost <= min_cost)
427  {
428  min_cost = cost;
429  min_arg = M;
430  }
431  }
432 
433  refine(tria, criteria, criteria(cell_indices[min_arg]));
434 }
435 
436 
437 // explicit instantiations
438 #include "grid_refinement.inst"
439 
unsigned int n_active_cells() const
Definition: tria.cc:11801
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)
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:11311
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:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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())
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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:11744
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
T min(const T &t, const MPI_Comm &mpi_communicator)
static const char N
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())
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)