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