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