Reference documentation for deal.II version GIT relicensing-1062-gc06da148b8 2024-07-15 19:20: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
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
25#include <deal.II/dofs/dof_accessor.templates.h>
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().size() == 0)
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().size() == 0)
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().size() == 0)
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().size() == 0)
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_communicator());
212 min_criterion_refine =
213 Utilities::MPI::min(min_criterion_refine,
214 parallel_tria->get_communicator());
215 max_criterion_coarsen =
216 Utilities::MPI::max(max_criterion_coarsen,
217 parallel_tria->get_communicator());
218 min_criterion_coarsen =
219 Utilities::MPI::min(min_criterion_coarsen,
220 parallel_tria->get_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().size() == 0)
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_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().size() == 0)
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().size() == 0)
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().size() == 0)
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#ifdef DEBUG
604 for (const auto &child : parent->child_iterators())
605 Assert(child->is_active() && child->coarsen_flag_set(),
606 typename Triangulation<
607 dim>::ExcInconsistentCoarseningFlags());
608#endif
609
610 parent_future_fe_index =
611 internal::hp::DoFHandlerImplementation::
612 dominated_future_fe_on_children<dim, spacedim>(parent);
613
614 future_fe_indices_on_coarsened_cells.insert(
615 {parent, parent_future_fe_index});
616 }
617 else
618 {
619 parent_future_fe_index =
620 future_fe_indices_on_coarsened_cells[parent];
621 }
622
623 future_fe_degree =
624 dof_handler.get_fe_collection()[parent_future_fe_index].degree;
625 }
626 else
627 {
628 // future finite element on current cell is already set
629 future_fe_degree =
630 dof_handler.get_fe_collection()[cell->future_fe_index()].degree;
631 }
632
633 // step 1: exponential decay with p-adaptation
634 if (cell->future_fe_index_set())
635 {
636 if (future_fe_degree > cell->get_fe().degree)
637 predicted_errors[cell->active_cell_index()] *=
638 Utilities::pow(gamma_p,
639 future_fe_degree - cell->get_fe().degree);
640 else if (future_fe_degree < cell->get_fe().degree)
641 predicted_errors[cell->active_cell_index()] /=
642 Utilities::pow(gamma_p,
643 cell->get_fe().degree - future_fe_degree);
644 else
645 {
646 // The two degrees are the same; we do not need to
647 // adapt the predicted error
648 }
649 }
650
651 // step 2: algebraic decay with h-adaptation
652 if (cell->refine_flag_set())
653 {
654 predicted_errors[cell->active_cell_index()] *=
655 (gamma_h * Utilities::pow(.5, future_fe_degree));
656
657 // predicted error will be split on children cells
658 // after adaptation via CellDataTransfer
659 }
660 else if (cell->coarsen_flag_set())
661 {
662 predicted_errors[cell->active_cell_index()] /=
663 (gamma_h * Utilities::pow(.5, future_fe_degree));
664
665 // predicted error will be summed up on parent cell
666 // after adaptation via CellDataTransfer
667 }
668 }
669 }
670
671
672
676 template <int dim, int spacedim>
677 void
679 {
680 if (dof_handler.get_fe_collection().size() == 0)
681 // nothing to do
682 return;
683
684 Assert(dof_handler.has_hp_capabilities(),
686
687 for (const auto &cell : dof_handler.active_cell_iterators())
688 if (cell->is_locally_owned() && cell->future_fe_index_set())
689 {
690 cell->clear_refine_flag();
691 cell->clear_coarsen_flag();
692 }
693 }
694
695
696
697 template <int dim, int spacedim>
698 void
700 {
701 if (dof_handler.get_fe_collection().size() == 0)
702 // nothing to do
703 return;
704
705 Assert(dof_handler.has_hp_capabilities(),
707
708 // Ghost siblings might occur on parallel Triangulation objects.
709 // We need information about refinement flags and future FE indices
710 // on all locally relevant cells here, and thus communicate them.
712 dynamic_cast<
714 const_cast<::Triangulation<dim, spacedim> *>(
715 &dof_handler.get_triangulation())))
716 {
719 }
720
722 const_cast<DoFHandler<dim, spacedim> &>(dof_handler));
723
724 // Now: choose p-adaptation over h-adaptation.
725 for (const auto &cell : dof_handler.active_cell_iterators())
726 if (cell->is_locally_owned() && cell->future_fe_index_set())
727 {
728 // This cell is flagged for p-adaptation.
729
730 // Remove any h-refinement flags.
731 cell->clear_refine_flag();
732
733 // A cell will only be coarsened into its parent if all of its
734 // siblings are flagged for h-coarsening as well. We must take this
735 // into account for our decision whether we would like to impose h-
736 // or p-adaptivity.
737 if (cell->coarsen_flag_set())
738 {
739 if (cell->level() == 0)
740 {
741 // This cell is a coarse cell and has neither parent nor
742 // siblings, thus it cannot be h-coarsened.
743 // Clear the flag and move on to the next cell.
744 cell->clear_coarsen_flag();
745 continue;
746 }
747
748 const auto &parent = cell->parent();
749 const unsigned int n_children = parent->n_children();
750
751 unsigned int h_flagged_children = 0, p_flagged_children = 0;
752 for (const auto &child : parent->child_iterators())
753 {
754 if (child->is_active())
755 {
756 Assert(child->is_artificial() == false,
758
759 if (child->coarsen_flag_set())
760 ++h_flagged_children;
761
762 // The public interface does not allow to access
763 // future FE indices on ghost cells. However, we
764 // need this information here and thus call the
765 // internal function that does not check for cell
766 // ownership.
767 if (internal::DoFCellAccessorImplementation::
768 Implementation::
769 future_fe_index_set<dim, spacedim, false>(
770 *child))
771 ++p_flagged_children;
772 }
773 }
774
775 if (h_flagged_children == n_children &&
776 p_flagged_children != n_children)
777 {
778 // Perform pure h-coarsening and
779 // drop all p-adaptation flags.
780 for (const auto &child : parent->child_iterators())
781 {
782 // h_flagged_children == n_children implies
783 // that all children are active
784 Assert(child->is_active(), ExcInternalError());
785 if (child->is_locally_owned())
786 child->clear_future_fe_index();
787 }
788 }
789 else
790 {
791 // Perform p-adaptation (if scheduled) and
792 // drop all h-coarsening flags.
793 for (const auto &child : parent->child_iterators())
794 {
795 if (child->is_active() && child->is_locally_owned())
796 child->clear_coarsen_flag();
797 }
798 }
799 }
800 }
801 }
802
803
804
808 template <int dim, int spacedim>
809 bool
811 const unsigned int max_difference,
812 const unsigned int contains_fe_index)
813 {
814 if (dof_handler.get_fe_collection().size() == 0)
815 // nothing to do
816 return false;
817
818 Assert(dof_handler.has_hp_capabilities(),
820 Assert(
821 max_difference > 0,
823 "This function does not serve any purpose for max_difference = 0."));
824 AssertIndexRange(contains_fe_index,
825 dof_handler.get_fe_collection().size());
826
827 //
828 // establish hierarchy
829 //
830 // - create bimap between hierarchy levels and FE indices
831
832 // there can be as many levels in the hierarchy as active FE indices are
833 // possible
834 using level_type = types::fe_index;
835 const auto invalid_level = static_cast<level_type>(-1);
836
837 // map from FE index to level in hierarchy
838 // FE indices that are not covered in the hierarchy are not in the map
839 const std::vector<unsigned int> fe_index_for_hierarchy_level =
841 contains_fe_index);
842
843 // map from level in hierarchy to FE index
844 // FE indices that are not covered in the hierarchy will be mapped to
845 // invalid_level
846 std::vector<unsigned int> hierarchy_level_for_fe_index(
847 dof_handler.get_fe_collection().size(), invalid_level);
848 for (unsigned int l = 0; l < fe_index_for_hierarchy_level.size(); ++l)
849 hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[l]] = l;
850
851
852 //
853 // parallelization
854 //
855 // - create distributed vector of level indices
856 // - update ghost values in each iteration (see later)
857 // - no need to compress, since the owning processor will have the correct
858 // level index
859
860 // HOTFIX: ::Vector does not accept integral types
862 if (const auto parallel_tria =
863 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
864 &(dof_handler.get_triangulation())))
865 {
866 future_levels.reinit(
867 parallel_tria->global_active_cell_index_partitioner().lock());
868 }
869 else
870 {
871 future_levels.reinit(
872 dof_handler.get_triangulation().n_active_cells());
873 }
874
875 for (const auto &cell : dof_handler.active_cell_iterators() |
877 future_levels[cell->global_active_cell_index()] =
878 hierarchy_level_for_fe_index[cell->future_fe_index()];
879
880
881 //
882 // limit level difference of neighboring cells
883 //
884 // - go over all locally relevant cells, and adjust the level indices of
885 // locally owned neighbors to match the level difference (as a
886 // consequence, indices on ghost cells will be updated only on the
887 // owning processor)
888 // - always raise levels to match criterion, never lower them
889 // - exchange level indices on ghost cells
890
891 // Function that updates the level of neighbor to fulfill difference
892 // criterion, and returns whether it was changed.
893 const auto update_neighbor_level =
894 [&future_levels, max_difference, invalid_level](
895 const auto &neighbor, const level_type cell_level) -> bool {
896 Assert(neighbor->is_active(), ExcInternalError());
897 // We only care about locally owned neighbors. If neighbor is a ghost
898 // cell, its future FE index will be updated on the owning process and
899 // communicated at the next loop iteration.
900 if (neighbor->is_locally_owned())
901 {
902 const level_type neighbor_level = static_cast<level_type>(
903 future_levels[neighbor->global_active_cell_index()]);
904
905 // ignore neighbors that are not part of the hierarchy
906 if (neighbor_level == invalid_level)
907 return false;
908
909 if ((cell_level - max_difference) > neighbor_level)
910 {
911 future_levels[neighbor->global_active_cell_index()] =
912 cell_level - max_difference;
913
914 return true;
915 }
916 }
917
918 return false;
919 };
920
921 // For cells to be h-coarsened, we need to determine a future FE for the
922 // parent cell, which will be the dominated FE among all children
923 // However, if we want to enforce the max_difference criterion on all
924 // cells on the updated mesh, we will need to simulate the updated mesh on
925 // the current mesh.
926 //
927 // As we are working on p-levels, we will set all siblings that will be
928 // coarsened to the highest p-level among them. The parent cell will be
929 // assigned exactly this level in form of the corresponding FE index in
930 // the adaptation process in
931 // Triangulation::execute_coarsening_and_refinement().
932 //
933 // This function takes a cell and sets all its siblings to the highest
934 // p-level among them. Returns whether any future levels have been
935 // changed.
936 const auto prepare_level_for_parent = [&](const auto &neighbor) -> bool {
937 Assert(neighbor->is_active(), ExcInternalError());
938 if (neighbor->coarsen_flag_set() && neighbor->is_locally_owned())
939 {
940 const auto parent = neighbor->parent();
941
942 std::vector<unsigned int> future_levels_children;
943 future_levels_children.reserve(parent->n_children());
944 for (const auto &child : parent->child_iterators())
945 {
946 Assert(child->is_active() && child->coarsen_flag_set(),
948 ExcInconsistentCoarseningFlags()));
949
950 const level_type child_level = static_cast<level_type>(
951 future_levels[child->global_active_cell_index()]);
952 Assert(child_level != invalid_level,
954 "The FiniteElement on one of the siblings of "
955 "a cell you are trying to coarsen is not part "
956 "of the registered p-adaptation hierarchy."));
957 future_levels_children.push_back(child_level);
958 }
959 Assert(!future_levels_children.empty(), ExcInternalError());
960
961 const unsigned int max_level_children =
962 *std::max_element(future_levels_children.begin(),
963 future_levels_children.end());
964
965 bool children_changed = false;
966 for (const auto &child : parent->child_iterators())
967 // We only care about locally owned children. If child is a ghost
968 // cell, its future FE index will be updated on the owning process
969 // and communicated at the next loop iteration.
970 if (child->is_locally_owned() &&
971 future_levels[child->global_active_cell_index()] !=
972 max_level_children)
973 {
974 future_levels[child->global_active_cell_index()] =
975 max_level_children;
976 children_changed = true;
977 }
978 return children_changed;
979 }
980
981 return false;
982 };
983
984 bool levels_changed = false;
985 bool levels_changed_in_cycle;
986 do
987 {
988 levels_changed_in_cycle = false;
989
990 future_levels.update_ghost_values();
991
992 for (const auto &cell : dof_handler.active_cell_iterators())
993 if (!cell->is_artificial())
994 {
995 const level_type cell_level = static_cast<level_type>(
996 future_levels[cell->global_active_cell_index()]);
997
998 // ignore cells that are not part of the hierarchy
999 if (cell_level == invalid_level)
1000 continue;
1001
1002 // ignore lowest levels of the hierarchy that always fulfill the
1003 // max_difference criterion
1004 if (cell_level <= max_difference)
1005 continue;
1006
1007 for (unsigned int f = 0; f < cell->n_faces(); ++f)
1008 if (cell->face(f)->at_boundary() == false)
1009 {
1010 if (cell->face(f)->has_children())
1011 {
1012 for (unsigned int sf = 0;
1013 sf < cell->face(f)->n_children();
1014 ++sf)
1015 {
1016 const auto neighbor =
1017 cell->neighbor_child_on_subface(f, sf);
1018
1019 levels_changed_in_cycle |=
1020 update_neighbor_level(neighbor, cell_level);
1021
1022 levels_changed_in_cycle |=
1023 prepare_level_for_parent(neighbor);
1024 }
1025 }
1026 else
1027 {
1028 const auto neighbor = cell->neighbor(f);
1029
1030 levels_changed_in_cycle |=
1031 update_neighbor_level(neighbor, cell_level);
1032
1033 levels_changed_in_cycle |=
1034 prepare_level_for_parent(neighbor);
1035 }
1036 }
1037 }
1038
1039 levels_changed_in_cycle =
1040 Utilities::MPI::logical_or(levels_changed_in_cycle,
1041 dof_handler.get_communicator());
1042 levels_changed |= levels_changed_in_cycle;
1043 }
1044 while (levels_changed_in_cycle);
1045
1046 // update future FE indices on locally owned cells
1047 for (const auto &cell : dof_handler.active_cell_iterators() |
1049 {
1050 const level_type cell_level = static_cast<level_type>(
1051 future_levels[cell->global_active_cell_index()]);
1052
1053 if (cell_level != invalid_level)
1054 {
1055 const unsigned int fe_index =
1056 fe_index_for_hierarchy_level[cell_level];
1057
1058 if (fe_index != cell->active_fe_index())
1059 cell->set_future_fe_index(fe_index);
1060 else
1061 cell->clear_future_fe_index();
1062 }
1063 }
1064
1065 return levels_changed;
1066 }
1067 } // namespace Refinement
1068} // namespace hp
1069
1070
1071// explicit instantiations
1072#include "refinement.inst"
1073
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
bool has_hp_capabilities() const
MPI_Comm get_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:264
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_communicator() const override
Definition tria_base.cc:160
#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
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)
#define DEAL_II_ASSERT_UNREACHABLE()
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:966
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:142
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
static const unsigned int invalid_unsigned_int
Definition types.h:220
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:59