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