Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2968-g5f01c80b02 2025-03-29 13:10:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
refinement.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2019 - 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
16#include <deal.II/base/config.h>
17
18#include <deal.II/base/mpi.h>
19
24
27
30
32
34#include <deal.II/lac/vector.h>
35
36#include <limits>
37
39
40namespace hp
41{
42 namespace Refinement
43 {
47 template <int dim, int spacedim>
48 void
50 {
51 if (dof_handler.get_fe_collection().empty())
52 // nothing to do
53 return;
54
55 Assert(dof_handler.has_hp_capabilities(),
57
58 std::vector<bool> p_flags(
59 dof_handler.get_triangulation().n_active_cells(), true);
60
61 p_adaptivity_from_flags(dof_handler, p_flags);
62 }
63
64
65
66 template <int dim, int spacedim>
67 void
69 const std::vector<bool> &p_flags)
70 {
71 if (dof_handler.get_fe_collection().empty())
72 // nothing to do
73 return;
74
75 Assert(dof_handler.has_hp_capabilities(),
78 p_flags.size());
79
80 for (const auto &cell : dof_handler.active_cell_iterators())
81 if (cell->is_locally_owned() && p_flags[cell->active_cell_index()])
82 {
83 if (cell->refine_flag_set())
84 {
85 const unsigned int super_fe_index =
87 cell->active_fe_index());
88
89 // Reject update if already most superordinate element.
90 if (super_fe_index != cell->active_fe_index())
91 cell->set_future_fe_index(super_fe_index);
92 }
93 else if (cell->coarsen_flag_set())
94 {
95 const unsigned int sub_fe_index =
97 cell->active_fe_index());
98
99 // Reject update if already least subordinate element.
100 if (sub_fe_index != cell->active_fe_index())
101 cell->set_future_fe_index(sub_fe_index);
102 }
103 }
104 }
105
106
107
108 template <int dim, typename Number, int spacedim>
109 void
111 const DoFHandler<dim, spacedim> &dof_handler,
112 const Vector<Number> &criteria,
113 const Number p_refine_threshold,
114 const Number p_coarsen_threshold,
116 &compare_refine,
118 &compare_coarsen)
119 {
120 if (dof_handler.get_fe_collection().empty())
121 // nothing to do
122 return;
123
124 Assert(dof_handler.has_hp_capabilities(),
127 criteria.size());
128
129 std::vector<bool> p_flags(
130 dof_handler.get_triangulation().n_active_cells(), false);
131
132 for (const auto &cell : dof_handler.active_cell_iterators())
133 if (cell->is_locally_owned() &&
134 ((cell->refine_flag_set() &&
135 compare_refine(criteria[cell->active_cell_index()],
136 p_refine_threshold)) ||
137 (cell->coarsen_flag_set() &&
138 compare_coarsen(criteria[cell->active_cell_index()],
139 p_coarsen_threshold))))
140 p_flags[cell->active_cell_index()] = true;
141
142 p_adaptivity_from_flags(dof_handler, p_flags);
143 }
144
145
146
147 template <int dim, typename Number, int spacedim>
148 void
150 const DoFHandler<dim, spacedim> &dof_handler,
151 const Vector<Number> &criteria,
152 const double p_refine_fraction,
153 const double p_coarsen_fraction,
155 &compare_refine,
157 &compare_coarsen)
158 {
159 if (dof_handler.get_fe_collection().empty())
160 // nothing to do
161 return;
162
163 Assert(dof_handler.has_hp_capabilities(),
166 criteria.size());
167 Assert((p_refine_fraction >= 0) && (p_refine_fraction <= 1),
169 Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
171
172 // We first have to determine the maximal and minimal values of the
173 // criteria of all flagged cells.
174 Number max_criterion_refine = std::numeric_limits<Number>::lowest(),
175 min_criterion_refine = std::numeric_limits<Number>::max();
176 Number max_criterion_coarsen = max_criterion_refine,
177 min_criterion_coarsen = min_criterion_refine;
178
179 for (const auto &cell : dof_handler.active_cell_iterators())
180 if (cell->is_locally_owned())
181 {
182 if (cell->refine_flag_set())
183 {
184 max_criterion_refine =
185 std::max(max_criterion_refine,
186 criteria(cell->active_cell_index()));
187 min_criterion_refine =
188 std::min(min_criterion_refine,
189 criteria(cell->active_cell_index()));
190 }
191 else if (cell->coarsen_flag_set())
192 {
193 max_criterion_coarsen =
194 std::max(max_criterion_coarsen,
195 criteria(cell->active_cell_index()));
196 min_criterion_coarsen =
197 std::min(min_criterion_coarsen,
198 criteria(cell->active_cell_index()));
199 }
200 }
201
202 const parallel::TriangulationBase<dim, spacedim> *parallel_tria =
203 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
204 &dof_handler.get_triangulation());
205 if (parallel_tria != nullptr &&
207 &dof_handler.get_triangulation()) == nullptr)
208 {
209 max_criterion_refine =
210 Utilities::MPI::max(max_criterion_refine,
211 parallel_tria->get_mpi_communicator());
212 min_criterion_refine =
213 Utilities::MPI::min(min_criterion_refine,
214 parallel_tria->get_mpi_communicator());
215 max_criterion_coarsen =
216 Utilities::MPI::max(max_criterion_coarsen,
217 parallel_tria->get_mpi_communicator());
218 min_criterion_coarsen =
219 Utilities::MPI::min(min_criterion_coarsen,
220 parallel_tria->get_mpi_communicator());
221 }
222
223 // Absent any better strategies, we will set the threshold by linear
224 // interpolation for both classes of cells individually.
225 const Number threshold_refine =
226 min_criterion_refine +
227 p_refine_fraction *
228 (max_criterion_refine - min_criterion_refine),
229 threshold_coarsen =
230 min_criterion_coarsen +
231 p_coarsen_fraction *
232 (max_criterion_coarsen - min_criterion_coarsen);
233
235 criteria,
236 threshold_refine,
237 threshold_coarsen,
238 compare_refine,
239 compare_coarsen);
240 }
241
242
243
244 template <int dim, typename Number, int spacedim>
245 void
247 const DoFHandler<dim, spacedim> &dof_handler,
248 const Vector<Number> &criteria,
249 const double p_refine_fraction,
250 const double p_coarsen_fraction,
252 &compare_refine,
254 &compare_coarsen)
255 {
256 if (dof_handler.get_fe_collection().empty())
257 // nothing to do
258 return;
259
260 Assert(dof_handler.has_hp_capabilities(),
263 criteria.size());
264 Assert((p_refine_fraction >= 0) && (p_refine_fraction <= 1),
266 Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
268
269 // ComparisonFunction returning 'true' or 'false' for any set of
270 // parameters. These will be used to overwrite user-provided comparison
271 // functions whenever no actual comparison is required in the decision
272 // process, i.e. when no or all cells will be refined or coarsened.
273 const ComparisonFunction<Number> compare_false =
274 [](const Number &, const Number &) { return false; };
275 const ComparisonFunction<Number> compare_true =
276 [](const Number &, const Number &) { return true; };
277
278 // 1.) First extract from the vector of indicators the ones that
279 // correspond to cells that we locally own.
280 unsigned int n_flags_refinement = 0;
281 unsigned int n_flags_coarsening = 0;
282 Vector<Number> indicators_refinement(
283 dof_handler.get_triangulation().n_active_cells());
284 Vector<Number> indicators_coarsening(
285 dof_handler.get_triangulation().n_active_cells());
286 for (const auto &cell :
288 if (!cell->is_artificial() && cell->is_locally_owned())
289 {
290 if (cell->refine_flag_set())
291 indicators_refinement(n_flags_refinement++) =
292 criteria(cell->active_cell_index());
293 else if (cell->coarsen_flag_set())
294 indicators_coarsening(n_flags_coarsening++) =
295 criteria(cell->active_cell_index());
296 }
297 indicators_refinement.grow_or_shrink(n_flags_refinement);
298 indicators_coarsening.grow_or_shrink(n_flags_coarsening);
299
300 // 2.) Determine the number of cells for p-refinement and p-coarsening on
301 // basis of the flagged cells.
302 //
303 // 3.) Find thresholds for p-refinement and p-coarsening on only those
304 // cells flagged for adaptation.
305 //
306 // For cases in which no or all cells flagged for refinement and/or
307 // coarsening are subject to p-adaptation, we usually pick thresholds
308 // that apply to all or none of the cells at once. However here, we
309 // do not know which threshold would suffice for this task because the
310 // user could provide any comparison function. Thus if necessary, we
311 // overwrite the user's choice with suitable functions simply
312 // returning 'true' and 'false' for any cell with reference wrappers.
313 // Thus, no function object copies are stored.
314 //
315 // 4.) Perform p-adaptation with absolute thresholds.
316 Number threshold_refinement = 0.;
317 Number threshold_coarsening = 0.;
318 auto reference_compare_refine = std::cref(compare_refine);
319 auto reference_compare_coarsen = std::cref(compare_coarsen);
320
321 const parallel::TriangulationBase<dim, spacedim> *parallel_tria =
322 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
323 &dof_handler.get_triangulation());
324 if (parallel_tria != nullptr &&
326 &dof_handler.get_triangulation()) == nullptr)
327 {
328#ifndef DEAL_II_WITH_P4EST
330#else
331 //
332 // parallel implementation with distributed memory
333 //
334
335 MPI_Comm mpi_communicator = parallel_tria->get_mpi_communicator();
336
337 // 2.) Communicate the number of cells scheduled for p-adaptation
338 // globally.
339 const unsigned int n_global_flags_refinement =
340 Utilities::MPI::sum(n_flags_refinement, mpi_communicator);
341 const unsigned int n_global_flags_coarsening =
342 Utilities::MPI::sum(n_flags_coarsening, mpi_communicator);
343
344 const unsigned int target_index_refinement =
345 static_cast<unsigned int>(
346 std::floor(p_refine_fraction * n_global_flags_refinement));
347 const unsigned int target_index_coarsening =
348 static_cast<unsigned int>(
349 std::ceil((1 - p_coarsen_fraction) * n_global_flags_coarsening));
350
351 // 3.) Figure out the global max and min of the criteria. We don't
352 // need it here, but it's a collective communication call.
353 const std::pair<Number, Number> global_min_max_refinement =
355 compute_global_min_and_max_at_root(indicators_refinement,
356 mpi_communicator);
357
358 const std::pair<Number, Number> global_min_max_coarsening =
360 compute_global_min_and_max_at_root(indicators_coarsening,
361 mpi_communicator);
362
363 // 3.) Compute thresholds if necessary.
364 if (target_index_refinement == 0)
365 reference_compare_refine = std::cref(compare_false);
366 else if (target_index_refinement == n_global_flags_refinement)
367 reference_compare_refine = std::cref(compare_true);
368 else
369 threshold_refinement = internal::parallel::distributed::
371 indicators_refinement,
372 global_min_max_refinement,
373 target_index_refinement,
374 mpi_communicator);
375
376 if (target_index_coarsening == n_global_flags_coarsening)
377 reference_compare_coarsen = std::cref(compare_false);
378 else if (target_index_coarsening == 0)
379 reference_compare_coarsen = std::cref(compare_true);
380 else
381 threshold_coarsening = internal::parallel::distributed::
383 indicators_coarsening,
384 global_min_max_coarsening,
385 target_index_coarsening,
386 mpi_communicator);
387#endif
388 }
389 else
390 {
391 //
392 // serial implementation (and parallel::shared implementation)
393 //
394
395 // 2.) Determine the number of cells scheduled for p-adaptation.
396 const unsigned int n_p_refine_cells = static_cast<unsigned int>(
397 std::floor(p_refine_fraction * n_flags_refinement));
398 const unsigned int n_p_coarsen_cells = static_cast<unsigned int>(
399 std::floor(p_coarsen_fraction * n_flags_coarsening));
400
401 // 3.) Compute thresholds if necessary.
402 if (n_p_refine_cells == 0)
403 reference_compare_refine = std::cref(compare_false);
404 else if (n_p_refine_cells == n_flags_refinement)
405 reference_compare_refine = std::cref(compare_true);
406 else
407 {
408 std::nth_element(indicators_refinement.begin(),
409 indicators_refinement.begin() +
410 n_p_refine_cells - 1,
411 indicators_refinement.end(),
412 std::greater<Number>());
413 threshold_refinement =
414 *(indicators_refinement.begin() + n_p_refine_cells - 1);
415 }
416
417 if (n_p_coarsen_cells == 0)
418 reference_compare_coarsen = std::cref(compare_false);
419 else if (n_p_coarsen_cells == n_flags_coarsening)
420 reference_compare_coarsen = std::cref(compare_true);
421 else
422 {
423 std::nth_element(indicators_coarsening.begin(),
424 indicators_coarsening.begin() +
425 n_p_coarsen_cells - 1,
426 indicators_coarsening.end(),
427 std::less<Number>());
428 threshold_coarsening =
429 *(indicators_coarsening.begin() + n_p_coarsen_cells - 1);
430 }
431 }
432
433 // 4.) Finally perform adaptation.
435 criteria,
436 threshold_refinement,
437 threshold_coarsening,
438 std::cref(reference_compare_refine),
439 std::cref(
440 reference_compare_coarsen));
441 }
442
443
444
445 template <int dim, typename Number, int spacedim>
446 void
448 const Vector<Number> &sobolev_indices)
449 {
450 if (dof_handler.get_fe_collection().empty())
451 // nothing to do
452 return;
453
454 Assert(dof_handler.has_hp_capabilities(),
457 sobolev_indices.size());
458
459 for (const auto &cell : dof_handler.active_cell_iterators())
460 if (cell->is_locally_owned())
461 {
462 if (cell->refine_flag_set())
463 {
464 const unsigned int super_fe_index =
466 cell->active_fe_index());
467
468 // Reject update if already most superordinate element.
469 if (super_fe_index != cell->active_fe_index())
470 {
471 const unsigned int super_fe_degree =
472 dof_handler.get_fe_collection()[super_fe_index].degree;
473
474 if (sobolev_indices[cell->active_cell_index()] >
475 super_fe_degree)
476 cell->set_future_fe_index(super_fe_index);
477 }
478 }
479 else if (cell->coarsen_flag_set())
480 {
481 const unsigned int sub_fe_index =
483 cell->active_fe_index());
484
485 // Reject update if already least subordinate element.
486 if (sub_fe_index != cell->active_fe_index())
487 {
488 const unsigned int sub_fe_degree =
489 dof_handler.get_fe_collection()[sub_fe_index].degree;
490
491 if (sobolev_indices[cell->active_cell_index()] <
492 sub_fe_degree)
493 cell->set_future_fe_index(sub_fe_index);
494 }
495 }
496 }
497 }
498
499
500
501 template <int dim, typename Number, int spacedim>
502 void
504 const DoFHandler<dim, spacedim> &dof_handler,
505 const Vector<Number> &criteria,
506 const Vector<Number> &references,
508 &compare_refine,
510 &compare_coarsen)
511 {
512 if (dof_handler.get_fe_collection().empty())
513 // nothing to do
514 return;
515
516 Assert(dof_handler.has_hp_capabilities(),
519 criteria.size());
521 references.size());
522
523 std::vector<bool> p_flags(
524 dof_handler.get_triangulation().n_active_cells(), false);
525
526 for (const auto &cell : dof_handler.active_cell_iterators())
527 if (cell->is_locally_owned() &&
528 ((cell->refine_flag_set() &&
529 compare_refine(criteria[cell->active_cell_index()],
530 references[cell->active_cell_index()])) ||
531 (cell->coarsen_flag_set() &&
532 compare_coarsen(criteria[cell->active_cell_index()],
533 references[cell->active_cell_index()]))))
534 p_flags[cell->active_cell_index()] = true;
535
536 p_adaptivity_from_flags(dof_handler, p_flags);
537 }
538
539
540
544 template <int dim, typename Number, int spacedim>
545 void
547 const Vector<Number> &error_indicators,
548 Vector<Number> &predicted_errors,
549 const double gamma_p,
550 const double gamma_h,
551 const double gamma_n)
552 {
553 if (dof_handler.get_fe_collection().empty())
554 // nothing to do
555 return;
556
558 error_indicators.size());
560 predicted_errors.size());
561 Assert(0 < gamma_p && gamma_p < 1,
565
566 // auxiliary variables
567 unsigned int future_fe_degree = numbers::invalid_unsigned_int;
568 unsigned int parent_future_fe_index = numbers::invalid_unsigned_int;
569 // store all determined future finite element indices on parent cells for
570 // coarsening
571 std::map<typename DoFHandler<dim, spacedim>::cell_iterator, unsigned int>
572 future_fe_indices_on_coarsened_cells;
573
574 // deep copy error indicators
575 predicted_errors = error_indicators;
576
577 for (const auto &cell : dof_handler.active_cell_iterators() |
579 {
580 // current cell will not be adapted
581 if (!(cell->future_fe_index_set()) && !(cell->refine_flag_set()) &&
582 !(cell->coarsen_flag_set()))
583 {
584 predicted_errors[cell->active_cell_index()] *= gamma_n;
585 continue;
586 }
587
588 // current cell will be adapted
589 // determine degree of its future finite element
590 if (cell->coarsen_flag_set())
591 {
592 Assert(cell->level() > 0,
593 ExcMessage("A coarse cell is flagged for coarsening. "
594 "Please read the note in the documentation "
595 "of predict_error()."));
596
597 // cell will be coarsened, thus determine future finite element
598 // on parent cell
599 const auto &parent = cell->parent();
600 if (future_fe_indices_on_coarsened_cells.find(parent) ==
601 future_fe_indices_on_coarsened_cells.end())
602 {
603 if constexpr (running_in_debug_mode())
604 {
605 for (const auto &child : parent->child_iterators())
606 Assert(child->is_active() && child->coarsen_flag_set(),
607 typename Triangulation<
608 dim>::ExcInconsistentCoarseningFlags());
609 }
610
611 parent_future_fe_index =
612 internal::hp::DoFHandlerImplementation::
613 dominated_future_fe_on_children<dim, spacedim>(parent);
614
615 future_fe_indices_on_coarsened_cells.insert(
616 {parent, parent_future_fe_index});
617 }
618 else
619 {
620 parent_future_fe_index =
621 future_fe_indices_on_coarsened_cells[parent];
622 }
623
624 future_fe_degree =
625 dof_handler.get_fe_collection()[parent_future_fe_index].degree;
626 }
627 else
628 {
629 // future finite element on current cell is already set
630 future_fe_degree =
631 dof_handler.get_fe_collection()[cell->future_fe_index()].degree;
632 }
633
634 // step 1: exponential decay with p-adaptation
635 if (cell->future_fe_index_set())
636 {
637 if (future_fe_degree > cell->get_fe().degree)
638 predicted_errors[cell->active_cell_index()] *=
639 Utilities::pow(gamma_p,
640 future_fe_degree - cell->get_fe().degree);
641 else if (future_fe_degree < cell->get_fe().degree)
642 predicted_errors[cell->active_cell_index()] /=
643 Utilities::pow(gamma_p,
644 cell->get_fe().degree - future_fe_degree);
645 else
646 {
647 // The two degrees are the same; we do not need to
648 // adapt the predicted error
649 }
650 }
651
652 // step 2: algebraic decay with h-adaptation
653 if (cell->refine_flag_set())
654 {
655 predicted_errors[cell->active_cell_index()] *=
656 (gamma_h * Utilities::pow(.5, future_fe_degree));
657
658 // predicted error will be split on children cells
659 // after adaptation via CellDataTransfer
660 }
661 else if (cell->coarsen_flag_set())
662 {
663 predicted_errors[cell->active_cell_index()] /=
664 (gamma_h * Utilities::pow(.5, future_fe_degree));
665
666 // predicted error will be summed up on parent cell
667 // after adaptation via CellDataTransfer
668 }
669 }
670 }
671
672
673
677 template <int dim, int spacedim>
678 void
680 {
681 if (dof_handler.get_fe_collection().empty())
682 // nothing to do
683 return;
684
685 Assert(dof_handler.has_hp_capabilities(),
687
688 for (const auto &cell : dof_handler.active_cell_iterators())
689 if (cell->is_locally_owned() && cell->future_fe_index_set())
690 {
691 cell->clear_refine_flag();
692 cell->clear_coarsen_flag();
693 }
694 }
695
696
697
698 template <int dim, int spacedim>
699 void
701 {
702 if (dof_handler.get_fe_collection().empty())
703 // nothing to do
704 return;
705
706 Assert(dof_handler.has_hp_capabilities(),
708
709 // Ghost siblings might occur on parallel Triangulation objects.
710 // We need information about refinement flags and future FE indices
711 // on all locally relevant cells here, and thus communicate them.
713 dynamic_cast<
715 const_cast<::Triangulation<dim, spacedim> *>(
716 &dof_handler.get_triangulation())))
717 {
720 }
721
723 const_cast<DoFHandler<dim, spacedim> &>(dof_handler));
724
725 // Now: choose p-adaptation over h-adaptation.
726 for (const auto &cell : dof_handler.active_cell_iterators())
727 if (cell->is_locally_owned() && cell->future_fe_index_set())
728 {
729 // This cell is flagged for p-adaptation.
730
731 // Remove any h-refinement flags.
732 cell->clear_refine_flag();
733
734 // A cell will only be coarsened into its parent if all of its
735 // siblings are flagged for h-coarsening as well. We must take this
736 // into account for our decision whether we would like to impose h-
737 // or p-adaptivity.
738 if (cell->coarsen_flag_set())
739 {
740 if (cell->level() == 0)
741 {
742 // This cell is a coarse cell and has neither parent nor
743 // siblings, thus it cannot be h-coarsened.
744 // Clear the flag and move on to the next cell.
745 cell->clear_coarsen_flag();
746 continue;
747 }
748
749 const auto &parent = cell->parent();
750 const unsigned int n_children = parent->n_children();
751
752 unsigned int h_flagged_children = 0, p_flagged_children = 0;
753 for (const auto &child : parent->child_iterators())
754 {
755 if (child->is_active())
756 {
757 Assert(child->is_artificial() == false,
759
760 if (child->coarsen_flag_set())
761 ++h_flagged_children;
762
763 // The public interface does not allow to access
764 // future FE indices on ghost cells. However, we
765 // need this information here and thus call the
766 // internal function that does not check for cell
767 // ownership.
769 Implementation::
770 future_fe_index_set<dim, spacedim, false>(
771 *child))
772 ++p_flagged_children;
773 }
774 }
775
776 if (h_flagged_children == n_children &&
777 p_flagged_children != n_children)
778 {
779 // Perform pure h-coarsening and
780 // drop all p-adaptation flags.
781 for (const auto &child : parent->child_iterators())
782 {
783 // h_flagged_children == n_children implies
784 // that all children are active
785 Assert(child->is_active(), ExcInternalError());
786 if (child->is_locally_owned())
787 child->clear_future_fe_index();
788 }
789 }
790 else
791 {
792 // Perform p-adaptation (if scheduled) and
793 // drop all h-coarsening flags.
794 for (const auto &child : parent->child_iterators())
795 {
796 if (child->is_active() && child->is_locally_owned())
797 child->clear_coarsen_flag();
798 }
799 }
800 }
801 }
802 }
803
804
805
809 template <int dim, int spacedim>
810 bool
812 const unsigned int max_difference,
813 const unsigned int contains_fe_index)
814 {
815 if (dof_handler.get_fe_collection().empty())
816 // nothing to do
817 return false;
818
819 Assert(dof_handler.has_hp_capabilities(),
821 Assert(
822 max_difference > 0,
824 "This function does not serve any purpose for max_difference = 0."));
825 AssertIndexRange(contains_fe_index,
826 dof_handler.get_fe_collection().size());
827
828 //
829 // establish hierarchy
830 //
831 // - create bimap between hierarchy levels and FE indices
832
833 // there can be as many levels in the hierarchy as active FE indices are
834 // possible
835 using level_type = types::fe_index;
836 const auto invalid_level = static_cast<level_type>(-1);
837
838 // map from FE index to level in hierarchy
839 // FE indices that are not covered in the hierarchy are not in the map
840 const std::vector<unsigned int> fe_index_for_hierarchy_level =
842 contains_fe_index);
843
844 // map from level in hierarchy to FE index
845 // FE indices that are not covered in the hierarchy will be mapped to
846 // invalid_level
847 std::vector<unsigned int> hierarchy_level_for_fe_index(
848 dof_handler.get_fe_collection().size(), invalid_level);
849 for (unsigned int l = 0; l < fe_index_for_hierarchy_level.size(); ++l)
850 hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[l]] = l;
851
852
853 //
854 // parallelization
855 //
856 // - create distributed vector of level indices
857 // - update ghost values in each iteration (see later)
858 // - no need to compress, since the owning processor will have the correct
859 // level index
860
861 // HOTFIX: ::Vector does not accept integral types
863 if (const auto parallel_tria =
864 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
865 &(dof_handler.get_triangulation())))
866 {
867 future_levels.reinit(
868 parallel_tria->global_active_cell_index_partitioner().lock());
869 }
870 else
871 {
872 future_levels.reinit(
873 dof_handler.get_triangulation().n_active_cells());
874 }
875
876 for (const auto &cell : dof_handler.active_cell_iterators() |
878 future_levels[cell->global_active_cell_index()] =
879 hierarchy_level_for_fe_index[cell->future_fe_index()];
880
881
882 //
883 // limit level difference of neighboring cells
884 //
885 // - go over all locally relevant cells, and adjust the level indices of
886 // locally owned neighbors to match the level difference (as a
887 // consequence, indices on ghost cells will be updated only on the
888 // owning processor)
889 // - always raise levels to match criterion, never lower them
890 // - exchange level indices on ghost cells
891
892 // Function that updates the level of neighbor to fulfill difference
893 // criterion, and returns whether it was changed.
894 const auto update_neighbor_level =
895 [&future_levels, max_difference, invalid_level](
896 const auto &neighbor, const level_type cell_level) -> bool {
897 Assert(neighbor->is_active(), ExcInternalError());
898 // We only care about locally owned neighbors. If neighbor is a ghost
899 // cell, its future FE index will be updated on the owning process and
900 // communicated at the next loop iteration.
901 if (neighbor->is_locally_owned())
902 {
903 const level_type neighbor_level = static_cast<level_type>(
904 future_levels[neighbor->global_active_cell_index()]);
905
906 // ignore neighbors that are not part of the hierarchy
907 if (neighbor_level == invalid_level)
908 return false;
909
910 if ((cell_level - max_difference) > neighbor_level)
911 {
912 future_levels[neighbor->global_active_cell_index()] =
913 cell_level - max_difference;
914
915 return true;
916 }
917 }
918
919 return false;
920 };
921
922 // For cells to be h-coarsened, we need to determine a future FE for the
923 // parent cell, which will be the dominated FE among all children
924 // However, if we want to enforce the max_difference criterion on all
925 // cells on the updated mesh, we will need to simulate the updated mesh on
926 // the current mesh.
927 //
928 // As we are working on p-levels, we will set all siblings that will be
929 // coarsened to the highest p-level among them. The parent cell will be
930 // assigned exactly this level in form of the corresponding FE index in
931 // the adaptation process in
932 // Triangulation::execute_coarsening_and_refinement().
933 //
934 // This function takes a cell and sets all its siblings to the highest
935 // p-level among them. Returns whether any future levels have been
936 // changed.
937 const auto prepare_level_for_parent = [&](const auto &neighbor) -> bool {
938 Assert(neighbor->is_active(), ExcInternalError());
939 if (neighbor->coarsen_flag_set() && neighbor->is_locally_owned())
940 {
941 const auto parent = neighbor->parent();
942
943 std::vector<unsigned int> future_levels_children;
944 future_levels_children.reserve(parent->n_children());
945 for (const auto &child : parent->child_iterators())
946 {
947 Assert(child->is_active() && child->coarsen_flag_set(),
949 ExcInconsistentCoarseningFlags()));
950
951 const level_type child_level = static_cast<level_type>(
952 future_levels[child->global_active_cell_index()]);
953 Assert(child_level != invalid_level,
955 "The FiniteElement on one of the siblings of "
956 "a cell you are trying to coarsen is not part "
957 "of the registered p-adaptation hierarchy."));
958 future_levels_children.push_back(child_level);
959 }
960 Assert(!future_levels_children.empty(), ExcInternalError());
961
962 const unsigned int max_level_children =
963 *std::max_element(future_levels_children.begin(),
964 future_levels_children.end());
965
966 bool children_changed = false;
967 for (const auto &child : parent->child_iterators())
968 // We only care about locally owned children. If child is a ghost
969 // cell, its future FE index will be updated on the owning process
970 // and communicated at the next loop iteration.
971 if (child->is_locally_owned() &&
972 future_levels[child->global_active_cell_index()] !=
973 max_level_children)
974 {
975 future_levels[child->global_active_cell_index()] =
976 max_level_children;
977 children_changed = true;
978 }
979 return children_changed;
980 }
981
982 return false;
983 };
984
985 bool levels_changed = false;
986 bool levels_changed_in_cycle;
987 do
988 {
989 levels_changed_in_cycle = false;
990
991 future_levels.update_ghost_values();
992
993 for (const auto &cell : dof_handler.active_cell_iterators())
994 if (!cell->is_artificial())
995 {
996 const level_type cell_level = static_cast<level_type>(
997 future_levels[cell->global_active_cell_index()]);
998
999 // ignore cells that are not part of the hierarchy
1000 if (cell_level == invalid_level)
1001 continue;
1002
1003 // ignore lowest levels of the hierarchy that always fulfill the
1004 // max_difference criterion
1005 if (cell_level <= max_difference)
1006 continue;
1007
1008 for (unsigned int f = 0; f < cell->n_faces(); ++f)
1009 if (cell->face(f)->at_boundary() == false)
1010 {
1011 if (cell->face(f)->has_children())
1012 {
1013 for (unsigned int sf = 0;
1014 sf < cell->face(f)->n_children();
1015 ++sf)
1016 {
1017 const auto neighbor =
1018 cell->neighbor_child_on_subface(f, sf);
1019
1020 levels_changed_in_cycle |=
1021 update_neighbor_level(neighbor, cell_level);
1022
1023 levels_changed_in_cycle |=
1024 prepare_level_for_parent(neighbor);
1025 }
1026 }
1027 else
1028 {
1029 const auto neighbor = cell->neighbor(f);
1030
1031 levels_changed_in_cycle |=
1032 update_neighbor_level(neighbor, cell_level);
1033
1034 levels_changed_in_cycle |=
1035 prepare_level_for_parent(neighbor);
1036 }
1037 }
1038 }
1039
1040 levels_changed_in_cycle =
1041 Utilities::MPI::logical_or(levels_changed_in_cycle,
1042 dof_handler.get_mpi_communicator());
1043 levels_changed |= levels_changed_in_cycle;
1044 }
1045 while (levels_changed_in_cycle);
1046
1047 // update future FE indices on locally owned cells
1048 for (const auto &cell : dof_handler.active_cell_iterators() |
1050 {
1051 const level_type cell_level = static_cast<level_type>(
1052 future_levels[cell->global_active_cell_index()]);
1053
1054 if (cell_level != invalid_level)
1055 {
1056 const unsigned int fe_index =
1057 fe_index_for_hierarchy_level[cell_level];
1058
1059 if (fe_index != cell->active_fe_index())
1060 cell->set_future_fe_index(fe_index);
1061 else
1062 cell->clear_future_fe_index();
1063 }
1064 }
1065
1066 return levels_changed;
1067 }
1068 } // namespace Refinement
1069} // namespace hp
1070
1071
1072// explicit instantiations
1073#include "hp/refinement.inst"
1074
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
bool has_hp_capabilities() const
MPI_Comm get_mpi_communicator() const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
unsigned int n_active_cells() const
virtual size_type size() const override
iterator end()
void grow_or_shrink(const size_type N)
iterator begin()
unsigned int size() const
Definition collection.h:316
bool empty() const
Definition collection.h:325
unsigned int previous_in_hierarchy(const unsigned int fe_index) const
std::vector< unsigned int > get_hierarchy_sequence(const unsigned int fe_index=0) const
unsigned int next_in_hierarchy(const unsigned int fe_index) const
virtual MPI_Comm get_mpi_communicator() const override
Definition tria_base.cc:160
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_ASSERT_UNREACHABLE()
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInvalidParameterValue()
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
T sum(const T &t, const MPI_Comm mpi_communicator)
T logical_or(const T &t, const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
void force_p_over_h(const DoFHandler< dim, spacedim > &dof_handler)
void p_adaptivity_from_reference(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const Vector< Number > &references, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen)
void full_p_adaptivity(const DoFHandler< dim, spacedim > &dof_handler)
Definition refinement.cc:49
void p_adaptivity_from_regularity(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &sobolev_indices)
bool limit_p_level_difference(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int max_difference=1, const unsigned int contains_fe_index=0)
void predict_error(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &error_indicators, Vector< Number > &predicted_errors, const double gamma_p=std::sqrt(0.4), const double gamma_h=2., const double gamma_n=1.)
void p_adaptivity_from_absolute_threshold(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const Number p_refine_threshold, const Number p_coarsen_threshold, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen=std::less_equal< Number >())
void choose_p_over_h(const DoFHandler< dim, spacedim > &dof_handler)
void p_adaptivity_fixed_number(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const double p_refine_fraction=0.5, const double p_coarsen_fraction=0.5, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen=std::less_equal< Number >())
std::function< bool(const Number &, const Number &)> ComparisonFunction
Definition refinement.h:143
void p_adaptivity_from_relative_threshold(const DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &criteria, const double p_refine_fraction=0.5, const double p_coarsen_fraction=0.5, const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_refine=std::greater_equal< Number >(), const ComparisonFunction< std_cxx20::type_identity_t< Number > > &compare_coarsen=std::less_equal< Number >())
void p_adaptivity_from_flags(const DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &p_flags)
Definition refinement.cc:68
Definition hp.h:117
void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
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 exchange_refinement_flags(::parallel::distributed::Triangulation< dim, spacedim > &tria)
Definition tria.cc:56
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
typename type_identity< T >::type type_identity_t
Definition type_traits.h:95
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:68