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