Reference documentation for deal.II version 9.4.0
grid_refinement.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2021 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
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
40namespace
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>(
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
166template <int dim, typename Number, int spacedim>
167void
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(),
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
216template <int dim, typename Number, int spacedim>
217void
219 const Vector<Number> & criteria,
220 const double threshold)
221{
222 Assert(criteria.size() == tria.n_active_cells(),
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
237template <int dim>
238std::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
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
278 coarsen_cells =
279 (current_n_cells - max_n_cells) / cell_decrease_on_coarsen;
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 }
316}
317
318
319
320template <int dim, typename Number, int spacedim>
321void
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 <=
339
340 const std::pair<double, double> adjusted_fractions =
342 max_n_cells,
343 top_fraction,
344 bottom_fraction);
345
346 const int refine_cells =
348 const int coarsen_cells =
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>::lowest());
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())
372 else
373 {
374 std::nth_element(tmp.begin(),
375 tmp.begin() + tmp.size() - coarsen_cells,
376 tmp.end(),
377 std::greater<double>());
379 criteria,
380 *(tmp.begin() + tmp.size() - coarsen_cells));
381 }
382 }
383 }
384}
385
386
387
388template <int dim, typename Number, int spacedim>
389void
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 <=
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
449template <int dim, typename Number, int spacedim>
450void
452 const Vector<Number> &criteria,
453 const unsigned int order)
454{
455 Assert(criteria.size() == tria.n_active_cells(),
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
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:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNegativeCriteria()
static ::ExceptionBase & ExcInvalidParameterValue()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const value_type * const_iterator
Definition: vector.h:125
iterator end()
size_type size() const
bool is_non_negative() const
bool all_zero() const
real_type l1_norm() const
iterator begin()
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:13734
static const unsigned int invalid_unsigned_int
Definition: types.h:201
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Predicate &predicate, const unsigned int grainsize)
Definition: parallel.h:135
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::Triangulation< dim, spacedim > & tria