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 - 2023 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
22#include <deal.II/lac/vector.h>
23
24#ifdef DEAL_II_WITH_P4EST
25
27
30# include <deal.II/grid/tria.h>
33
34# include <algorithm>
35# include <functional>
36# include <limits>
37# include <numeric>
38
39
41
42
43namespace
44{
45 template <int dim, int spacedim>
46 unsigned int
47 n_locally_owned_active_cells(const Triangulation<dim, spacedim> &tria)
48 {
49 if (const auto parallel_tria =
51 &tria))
52 return parallel_tria->n_locally_owned_active_cells();
53 else
54 return tria.n_active_cells();
55 }
56
57 template <typename number>
58 inline number
59 max_element(const ::Vector<number> &criteria)
60 {
61 return (criteria.size() > 0) ?
62 (*std::max_element(criteria.begin(), criteria.end())) :
63 std::numeric_limits<number>::min();
64 }
65
66
67
68 template <typename number>
69 inline number
70 min_element(const ::Vector<number> &criteria)
71 {
72 return (criteria.size() > 0) ?
73 (*std::min_element(criteria.begin(), criteria.end())) :
74 std::numeric_limits<number>::max();
75 }
76
77
78
84 template <typename number>
85 double
86 compute_global_sum(const ::Vector<number> &criteria,
87 const MPI_Comm mpi_communicator)
88 {
89 double my_sum =
90 std::accumulate(criteria.begin(),
91 criteria.end(),
92 /* do accumulation in the correct data type: */
93 number());
94
95 double result = 0;
96 // compute the minimum on processor zero
97 const int ierr =
98 MPI_Reduce(&my_sum, &result, 1, MPI_DOUBLE, MPI_SUM, 0, mpi_communicator);
99 AssertThrowMPI(ierr);
100
101 // make sure only processor zero got something
102 if (Utilities::MPI::this_mpi_process(mpi_communicator) != 0)
103 Assert(result == 0, ExcInternalError());
104
105 return result;
106 }
107
108
109
114 template <int dim, int spacedim, typename Number>
115 void
116 get_locally_owned_indicators(const ::Triangulation<dim, spacedim> &tria,
117 const ::Vector<Number> &criteria,
118 Vector<Number> &locally_owned_indicators)
119 {
120 Assert(locally_owned_indicators.size() ==
121 n_locally_owned_active_cells(tria),
123
124 unsigned int owned_index = 0;
125 for (const auto &cell :
126 tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
127 {
128 locally_owned_indicators(owned_index) =
129 criteria(cell->active_cell_index());
130 ++owned_index;
131 }
132 Assert(owned_index == n_locally_owned_active_cells(tria),
134 }
135
136
137 // we compute refinement thresholds by bisection of the interval spanned by
138 // the smallest and largest error indicator. this leads to a small problem:
139 // if, for example, we want to coarsen zero per cent of the cells, then we
140 // need to pick a threshold equal to the smallest indicator, but of course
141 // the bisection algorithm can never find a threshold equal to one of the
142 // end points of the interval. So we slightly increase the interval before
143 // we even start
144 void
145 adjust_interesting_range(double (&interesting_range)[2])
146 {
147 Assert(interesting_range[0] <= interesting_range[1], ExcInternalError());
148
149 if (interesting_range[0] > 0)
150 {
151 // In this case, we calculate the first interval split point `m` in the
152 // `compute_threshold` functions in the optimized way: We exploit that
153 // the logarithms of all criteria are more uniformly distributed than
154 // their actual values, i.e. m=sqrt(b*e).
155 //
156 // Both factors will modify the split point only slightly by a factor of
157 // sqrt(1.01*0.99) = sqrt(0.9999) ~ 0.9950.
158 interesting_range[0] *= 0.99;
159 interesting_range[1] *= 1.01;
160 }
161 else
162 {
163 // In all other cases, we begin with an the arithmetic mean as the
164 // standard interval split point, i.e. m=(b+e)/2.
165 //
166 // Both increments will add up to zero when calculating the initial
167 // split point in the `compute_threshold` functions.
168 const double difference =
169 std::abs(interesting_range[1] - interesting_range[0]);
170 interesting_range[0] -= 0.01 * difference;
171 interesting_range[1] += 0.01 * difference;
172 }
173 }
174
175
176
182 template <int dim, int spacedim, typename Number>
183 void
185 const ::Vector<Number> & criteria,
186 const double top_threshold,
187 const double bottom_threshold)
188 {
189 ::GridRefinement::refine(tria, criteria, top_threshold);
190 ::GridRefinement::coarsen(tria, criteria, bottom_threshold);
191
192 // as a final good measure, delete all flags again from cells that we don't
193 // locally own
194 for (const auto &cell : tria.active_cell_iterators())
195 if (cell->subdomain_id() != tria.locally_owned_subdomain())
196 {
197 cell->clear_refine_flag();
198 cell->clear_coarsen_flag();
199 }
200 }
201
202
203
211 template <int dim, int spacedim, typename Number>
212 void
213 refine_and_coarsen_fixed_fraction_via_l1_norm(
215 const ::Vector<Number> & criteria,
216 const double top_fraction_of_error,
217 const double bottom_fraction_of_error)
218 {
219 // first extract from the vector of indicators the ones that correspond
220 // to cells that we locally own
221 Vector<Number> locally_owned_indicators(n_locally_owned_active_cells(tria));
222 get_locally_owned_indicators(tria, criteria, locally_owned_indicators);
223
224 MPI_Comm mpi_communicator = tria.get_communicator();
225
226 // figure out the global max and min of the indicators. we don't need it
227 // here, but it's a collective communication call
228 const std::pair<double, double> global_min_and_max =
230 compute_global_min_and_max_at_root(locally_owned_indicators,
231 mpi_communicator);
232
233 const double total_error =
234 compute_global_sum(locally_owned_indicators, mpi_communicator);
235
236 double top_target_error = top_fraction_of_error * total_error,
237 bottom_target_error = (1. - bottom_fraction_of_error) * total_error;
238
239 double top_threshold, bottom_threshold;
242 global_min_and_max,
243 top_target_error,
244 mpi_communicator);
245
246 // compute bottom threshold only if necessary. otherwise use the lowest
247 // threshold possible
248 if (bottom_fraction_of_error > 0)
249 bottom_threshold = ::internal::parallel::distributed::
251 locally_owned_indicators,
252 global_min_and_max,
253 bottom_target_error,
254 mpi_communicator);
255 else
256 bottom_threshold = std::numeric_limits<Number>::lowest();
257
258 // now refine the mesh
259 mark_cells(tria, criteria, top_threshold, bottom_threshold);
260 }
261} // namespace
262
263
264
265namespace internal
266{
267 namespace parallel
268 {
269 namespace distributed
270 {
271 namespace GridRefinement
272 {
273 template <typename number>
274 std::pair<number, number>
276 const ::Vector<number> &criteria,
277 const MPI_Comm mpi_communicator)
278 {
279 // we'd like to compute the global max and min from the local ones in
280 // one MPI communication. we can do that by taking the elementwise
281 // minimum of the local min and the negative maximum over all
282 // processors
283
284 const double local_min = min_element(criteria),
285 local_max = max_element(criteria);
286 double comp[2] = {local_min, -local_max};
287 double result[2] = {0, 0};
288
289 // compute the minimum on processor zero
290 const int ierr = MPI_Reduce(
291 comp, result, 2, MPI_DOUBLE, MPI_MIN, 0, mpi_communicator);
292 AssertThrowMPI(ierr);
293
294 // make sure only processor zero got something
295 if (Utilities::MPI::this_mpi_process(mpi_communicator) != 0)
296 Assert((result[0] == 0) && (result[1] == 0), ExcInternalError());
297
298 return std::make_pair(result[0], -result[1]);
299 }
300
301
302
303 namespace RefineAndCoarsenFixedNumber
304 {
305 template <typename number>
306 number
307 compute_threshold(const ::Vector<number> & criteria,
308 const std::pair<double, double> &global_min_and_max,
309 const types::global_cell_index n_target_cells,
310 const MPI_Comm mpi_communicator)
311 {
312 double interesting_range[2] = {global_min_and_max.first,
313 global_min_and_max.second};
314 adjust_interesting_range(interesting_range);
315
316 const unsigned int root_mpi_rank = 0;
317 unsigned int iteration = 0;
318
319 do
320 {
321 int ierr = MPI_Bcast(interesting_range,
322 2,
323 MPI_DOUBLE,
324 root_mpi_rank,
325 mpi_communicator);
326 AssertThrowMPI(ierr);
327
328 if (interesting_range[0] == interesting_range[1])
329 return interesting_range[0];
330
331 const double test_threshold =
332 (interesting_range[0] > 0 ?
333 std::sqrt(interesting_range[0] * interesting_range[1]) :
334 (interesting_range[0] + interesting_range[1]) / 2);
335
336 // Count how many of our own elements would be above this
337 // threshold. Use a 64bit result type if we are compiling with
338 // 64bit indices to avoid an overflow when computing the sum
339 // below.
340 const types::global_cell_index my_count =
341 std::count_if(criteria.begin(),
342 criteria.end(),
343 [test_threshold](const double c) {
344 return c > test_threshold;
345 });
346 const types::global_cell_index total_count =
347 Utilities::MPI::sum(my_count, mpi_communicator);
348
349 // now adjust the range. if we have too many cells, we take the
350 // upper half of the previous range, otherwise the lower half.
351 // if we have hit the right number, then set the range to the
352 // exact value. non-root nodes also update their own
353 // interesting_range, however their results are not significant
354 // since the values will be overwritten by MPI_Bcast from the
355 // root node in next loop.
356 if (total_count > n_target_cells)
357 interesting_range[0] = test_threshold;
358 else if (total_count < n_target_cells)
359 interesting_range[1] = test_threshold;
360 else
361 interesting_range[0] = interesting_range[1] = test_threshold;
362
363 // terminate the iteration after 25 go-arounds. this is
364 // necessary because oftentimes error indicators on cells have
365 // exactly the same value, and so there may not be a particular
366 // value that cuts the indicators in such a way that we can
367 // achieve the desired number of cells. using a maximal number
368 // of iterations means that we terminate the iteration after a
369 // fixed number N of steps if the indicators were perfectly
370 // badly distributed, and we make at most a mistake of 1/2^N in
371 // the number of cells flagged if indicators are perfectly
372 // equidistributed
373 ++iteration;
374 if (iteration == 25)
375 interesting_range[0] = interesting_range[1] = test_threshold;
376 }
377 while (true);
378
379 Assert(false, ExcInternalError());
380 return -1;
381 }
382 } // namespace RefineAndCoarsenFixedNumber
383
384
385
386 namespace RefineAndCoarsenFixedFraction
387 {
388 template <typename number>
389 number
390 compute_threshold(const ::Vector<number> & criteria,
391 const std::pair<double, double> &global_min_and_max,
392 const double target_error,
393 const MPI_Comm mpi_communicator)
394 {
395 double interesting_range[2] = {global_min_and_max.first,
396 global_min_and_max.second};
397 adjust_interesting_range(interesting_range);
398
399 const unsigned int root_mpi_rank = 0;
400 unsigned int iteration = 0;
401
402 do
403 {
404 int ierr = MPI_Bcast(interesting_range,
405 2,
406 MPI_DOUBLE,
407 root_mpi_rank,
408 mpi_communicator);
409 AssertThrowMPI(ierr);
410
411 if (interesting_range[0] == interesting_range[1])
412 {
413 // so we have found our threshold. since we adjust the range
414 // at the top of the function to be slightly larger than the
415 // actual extremes of the refinement criteria values, we can
416 // end up in a situation where the threshold is in fact
417 // larger than the maximal refinement indicator. in such
418 // cases, we get no refinement at all. thus, cap the
419 // threshold by the actual largest value
420 double final_threshold =
421 std::min(interesting_range[0], global_min_and_max.second);
422 ierr = MPI_Bcast(&final_threshold,
423 1,
424 MPI_DOUBLE,
425 root_mpi_rank,
426 mpi_communicator);
427 AssertThrowMPI(ierr);
428
429 return final_threshold;
430 }
431
432 const double test_threshold =
433 (interesting_range[0] > 0 ?
434 std::sqrt(interesting_range[0] * interesting_range[1]) :
435 (interesting_range[0] + interesting_range[1]) / 2);
436
437 // accumulate the error of those our own elements above this
438 // threshold and then add to it the number for all the others
439 double my_error = 0;
440 for (unsigned int i = 0; i < criteria.size(); ++i)
441 if (criteria(i) > test_threshold)
442 my_error += criteria(i);
443
444 double total_error = 0.;
445
446 ierr = MPI_Reduce(&my_error,
447 &total_error,
448 1,
449 MPI_DOUBLE,
450 MPI_SUM,
451 root_mpi_rank,
452 mpi_communicator);
453 AssertThrowMPI(ierr);
454
455 // now adjust the range. if we have too many cells, we take the
456 // upper half of the previous range, otherwise the lower half.
457 // if we have hit the right number, then set the range to the
458 // exact value. non-root nodes also update their own
459 // interesting_range, however their results are not significant
460 // since the values will be overwritten by MPI_Bcast from the
461 // root node in next loop.
462 if (total_error > target_error)
463 interesting_range[0] = test_threshold;
464 else if (total_error < target_error)
465 interesting_range[1] = test_threshold;
466 else
467 interesting_range[0] = interesting_range[1] = test_threshold;
468
469 // terminate the iteration after 25 go-arounds. this is
470 // necessary because oftentimes error indicators on cells
471 // have exactly the same value, and so there may not be a
472 // particular value that cuts the indicators in such a way
473 // that we can achieve the desired number of cells. using a
474 // max of 25 iterations means that we terminate the
475 // iteration after 25 steps if the indicators were perfectly
476 // badly distributed, and we make at most a mistake of
477 // 1/2^25 in the number of cells flagged if indicators are
478 // perfectly equidistributed
479 ++iteration;
480 if (iteration == 25)
481 interesting_range[0] = interesting_range[1] = test_threshold;
482 }
483 while (true);
484
485 Assert(false, ExcInternalError());
486 return -1;
487 }
488 } // namespace RefineAndCoarsenFixedFraction
489 } // namespace GridRefinement
490 } // namespace distributed
491 } // namespace parallel
492} // namespace internal
493
494
495
496namespace parallel
497{
498 namespace distributed
499 {
500 namespace GridRefinement
501 {
502 template <int dim, typename Number, int spacedim>
503 void
506 const ::Vector<Number> & criteria,
507 const double top_fraction_of_cells,
508 const double bottom_fraction_of_cells,
509 const types::global_cell_index max_n_cells)
510 {
511 Assert(criteria.size() == tria.n_active_cells(),
512 ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
513 Assert((top_fraction_of_cells >= 0) && (top_fraction_of_cells <= 1),
515 Assert((bottom_fraction_of_cells >= 0) &&
516 (bottom_fraction_of_cells <= 1),
518 Assert(top_fraction_of_cells + bottom_fraction_of_cells <= 1,
520 Assert(criteria.is_non_negative(),
522
523 const std::pair<double, double> adjusted_fractions =
526 max_n_cells,
527 top_fraction_of_cells,
528 bottom_fraction_of_cells);
529
530 // first extract from the vector of indicators the ones that correspond
531 // to cells that we locally own
532 Vector<Number> locally_owned_indicators(
533 n_locally_owned_active_cells(tria));
534 get_locally_owned_indicators(tria, criteria, locally_owned_indicators);
535
536 MPI_Comm mpi_communicator = tria.get_communicator();
537
538 // figure out the global max and min of the indicators. we don't need it
539 // here, but it's a collective communication call
540 const std::pair<Number, Number> global_min_and_max =
542 compute_global_min_and_max_at_root(locally_owned_indicators,
543 mpi_communicator);
544
545
546 double top_threshold, bottom_threshold;
549 locally_owned_indicators,
550 global_min_and_max,
551 static_cast<types::global_cell_index>(adjusted_fractions.first *
553 mpi_communicator);
554
555 // compute bottom threshold only if necessary. otherwise use the lowest
556 // threshold possible
557 if (adjusted_fractions.second > 0)
558 bottom_threshold = ::internal::parallel::distributed::
560 locally_owned_indicators,
561 global_min_and_max,
562 static_cast<types::global_cell_index>(
563 std::ceil((1. - adjusted_fractions.second) *
565 mpi_communicator);
566 else
567 bottom_threshold = std::numeric_limits<Number>::lowest();
568
569 // now refine the mesh
570 mark_cells(tria, criteria, top_threshold, bottom_threshold);
571 }
572
573
574
575 template <int dim, typename Number, int spacedim>
576 void
579 const ::Vector<Number> & criteria,
580 const double top_fraction_of_error,
581 const double bottom_fraction_of_error,
582 const VectorTools::NormType norm_type)
583 {
584 Assert(criteria.size() == tria.n_active_cells(),
585 ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
586 Assert((top_fraction_of_error >= 0) && (top_fraction_of_error <= 1),
588 Assert((bottom_fraction_of_error >= 0) &&
589 (bottom_fraction_of_error <= 1),
591 Assert(top_fraction_of_error + bottom_fraction_of_error <= 1,
593 Assert(criteria.is_non_negative(),
595
596 switch (norm_type)
597 {
599 // evaluate norms on subsets and compare them as
600 // c_0 + c_1 + ... < fraction * l1-norm(c)
601 refine_and_coarsen_fixed_fraction_via_l1_norm(
602 tria,
603 criteria,
604 top_fraction_of_error,
605 bottom_fraction_of_error);
606 break;
607
609 {
610 // we do not want to evaluate norms on subsets as:
611 // sqrt(c_0^2 + c_1^2 + ...) < fraction * l2-norm(c)
612 // instead take the square of both sides of the equation
613 // and evaluate:
614 // c_0^2 + c_1^2 + ... < fraction^2 * l1-norm(c.c)
615 // we adjust all parameters accordingly
616 Vector<Number> criteria_squared(criteria.size());
617 std::transform(criteria.begin(),
618 criteria.end(),
619 criteria_squared.begin(),
620 [](Number c) { return c * c; });
621
622 refine_and_coarsen_fixed_fraction_via_l1_norm(
623 tria,
624 criteria_squared,
625 top_fraction_of_error * top_fraction_of_error,
626 bottom_fraction_of_error * bottom_fraction_of_error);
627 }
628 break;
629
630 default:
631 Assert(false, ExcNotImplemented());
632 break;
633 }
634 }
635 } // namespace GridRefinement
636 } // namespace distributed
637} // namespace parallel
638
639
640// explicit instantiations
641# include "grid_refinement.inst"
642
644
645#endif
virtual types::global_cell_index n_global_active_cells() const
virtual MPI_Comm get_communicator() const
unsigned int n_active_cells() const
size_type size() const
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertThrowMPI(error_code)
static ::ExceptionBase & ExcNegativeCriteria()
static ::ExceptionBase & ExcInvalidParameterValue()
static ::ExceptionBase & ExcInternalError()
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 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)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161
number compute_threshold(const ::Vector< number > &criteria, const std::pair< double, double > &global_min_and_max, const double target_error, const MPI_Comm mpi_communicator)
number compute_threshold(const ::Vector< number > &criteria, const std::pair< double, double > &global_min_and_max, const types::global_cell_index n_target_cells, const MPI_Comm mpi_communicator)
std::pair< number, number > compute_global_min_and_max_at_root(const ::Vector< number > &criteria, const MPI_Comm mpi_communicator)
void refine_and_coarsen_fixed_fraction(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_error, const double bottom_fraction_of_error, const VectorTools::NormType norm_type=VectorTools::NormType::L1_norm)
void refine_and_coarsen_fixed_number(::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())
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int subdomain_id
Definition types.h:44
const ::Triangulation< dim, spacedim > & tria