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