Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
refinement.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2019 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#include <deal.II/base/config.h>
18
19#include <deal.II/base/mpi.h>
20
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 simplying
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
329 Assert(false, ExcInternalError());
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 predicted_errors[cell->active_cell_index()] *=
637 std::pow(gamma_p,
638 int(future_fe_degree) - int(cell->get_fe().degree));
639 }
640
641 // step 2: algebraic decay with h-adaptation
642 if (cell->refine_flag_set())
643 {
644 predicted_errors[cell->active_cell_index()] *=
645 (gamma_h * std::pow(.5, future_fe_degree));
646
647 // predicted error will be split on children cells
648 // after adaptation via CellDataTransfer
649 }
650 else if (cell->coarsen_flag_set())
651 {
652 predicted_errors[cell->active_cell_index()] /=
653 (gamma_h * std::pow(.5, future_fe_degree));
654
655 // predicted error will be summed up on parent cell
656 // after adaptation via CellDataTransfer
657 }
658 }
659 }
660
661
662
666 template <int dim, int spacedim>
667 void
669 {
670 if (dof_handler.get_fe_collection().size() == 0)
671 // nothing to do
672 return;
673
674 Assert(dof_handler.has_hp_capabilities(),
676
677 for (const auto &cell : dof_handler.active_cell_iterators())
678 if (cell->is_locally_owned() && cell->future_fe_index_set())
679 {
680 cell->clear_refine_flag();
681 cell->clear_coarsen_flag();
682 }
683 }
684
685
686
687 template <int dim, int spacedim>
688 void
690 {
691 if (dof_handler.get_fe_collection().size() == 0)
692 // nothing to do
693 return;
694
695 Assert(dof_handler.has_hp_capabilities(),
697
698 // Ghost siblings might occur on parallel::shared::Triangulation objects.
699 // We need information about future FE indices on all locally relevant
700 // cells here, and thus communicate them.
701 if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
702 &dof_handler.get_triangulation()) != nullptr)
704 const_cast<DoFHandler<dim, spacedim> &>(dof_handler));
705
706 for (const auto &cell : dof_handler.active_cell_iterators())
707 if (cell->is_locally_owned() && cell->future_fe_index_set())
708 {
709 cell->clear_refine_flag();
710
711 // A cell will only be coarsened into its parent if all of its
712 // siblings are flagged for h-coarsening as well. We must take this
713 // into account for our decision whether we would like to impose h-
714 // or p-adaptivity.
715 if (cell->coarsen_flag_set())
716 {
717 if (cell->level() == 0)
718 {
719 // This cell is a coarse cell and has neither parent nor
720 // siblings, thus it cannot be h-coarsened.
721 // Clear the flag and move on to the next cell.
722 cell->clear_coarsen_flag();
723 continue;
724 }
725
726 const auto & parent = cell->parent();
727 const unsigned int n_children = parent->n_children();
728
729 unsigned int h_flagged_children = 0, p_flagged_children = 0;
730 for (const auto &child : parent->child_iterators())
731 {
732 if (child->is_active())
733 {
734 if (child->is_locally_owned())
735 {
736 if (child->coarsen_flag_set())
737 ++h_flagged_children;
738 if (child->future_fe_index_set())
739 ++p_flagged_children;
740 }
741 else if (child->is_ghost())
742 {
743 // The case of siblings being owned by different
744 // processors can only occur for
745 // parallel::shared::Triangulation objects.
746 Assert(
747 (dynamic_cast<const parallel::shared::
748 Triangulation<dim, spacedim> *>(
749 &dof_handler.get_triangulation()) != nullptr),
751
752 if (child->coarsen_flag_set())
753 ++h_flagged_children;
754 // The public interface does not allow to access
755 // future FE indices on ghost cells. However, we
756 // need this information here and thus call the
757 // internal function that does not check for cell
758 // ownership.
759 if (internal::DoFCellAccessorImplementation::
760 Implementation::
761 future_fe_index_set<dim, spacedim, false>(
762 *child))
763 ++p_flagged_children;
764 }
765 else
766 {
767 // Siblings of locally owned cells are all
768 // either also locally owned or ghost cells.
769 Assert(false, ExcInternalError());
770 }
771 }
772 }
773
774 if (h_flagged_children == n_children &&
775 p_flagged_children != n_children)
776 {
777 // Perform pure h-coarsening and
778 // drop all p-adaptation flags.
779 for (const auto &child : parent->child_iterators())
780 {
781 // h_flagged_children == n_children implies
782 // that all children are active
783 Assert(child->is_active(), ExcInternalError());
784 if (child->is_locally_owned())
785 child->clear_future_fe_index();
786 }
787 }
788 else
789 {
790 // Perform p-adaptation on all children and
791 // drop all h-coarsening flags.
792 for (const auto &child : parent->child_iterators())
793 {
794 if (child->is_active() && child->is_locally_owned())
795 child->clear_coarsen_flag();
796 }
797 }
798 }
799 }
800 }
801
802
803
807 template <int dim, int spacedim>
808 bool
810 const unsigned int max_difference,
811 const unsigned int contains_fe_index)
812 {
813 if (dof_handler.get_fe_collection().size() == 0)
814 // nothing to do
815 return false;
816
817 Assert(dof_handler.has_hp_capabilities(),
819 Assert(
820 max_difference > 0,
822 "This function does not serve any purpose for max_difference = 0."));
823 AssertIndexRange(contains_fe_index,
824 dof_handler.get_fe_collection().size());
825
826 //
827 // establish hierarchy
828 //
829 // - create bimap between hierarchy levels and FE indices
830
831 // there can be as many levels in the hierarchy as active FE indices are
832 // possible
833 using level_type = types::fe_index;
834 const auto invalid_level = static_cast<level_type>(-1);
835
836 // map from FE index to level in hierarchy
837 // FE indices that are not covered in the hierarchy are not in the map
838 const std::vector<unsigned int> fe_index_for_hierarchy_level =
840 contains_fe_index);
841
842 // map from level in hierarchy to FE index
843 // FE indices that are not covered in the hierarchy will be mapped to
844 // invalid_level
845 std::vector<unsigned int> hierarchy_level_for_fe_index(
846 dof_handler.get_fe_collection().size(), invalid_level);
847 for (unsigned int l = 0; l < fe_index_for_hierarchy_level.size(); ++l)
848 hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[l]] = l;
849
850
851 //
852 // parallelization
853 //
854 // - create distributed vector of level indices
855 // - update ghost values in each iteration (see later)
856 // - no need to compress, since the owning processor will have the correct
857 // level index
858
859 // HOTFIX: ::Vector does not accept integral types
861 if (const auto parallel_tria =
862 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
863 &(dof_handler.get_triangulation())))
864 {
865 future_levels.reinit(
866 parallel_tria->global_active_cell_index_partitioner().lock());
867 }
868 else
869 {
870 future_levels.reinit(
871 dof_handler.get_triangulation().n_active_cells());
872 }
873
874 for (const auto &cell : dof_handler.active_cell_iterators() |
876 future_levels[cell->global_active_cell_index()] =
877 hierarchy_level_for_fe_index[cell->future_fe_index()];
878
879
880 //
881 // limit level difference of neighboring cells
882 //
883 // - go over all locally relevant cells, and adjust the level indices of
884 // locally owned neighbors to match the level difference (as a
885 // consequence, indices on ghost cells will be updated only on the
886 // owning processor)
887 // - always raise levels to match criterion, never lower them
888 // - exchange level indices on ghost cells
889
890 // Function that updates the level of neighbor to fulfill difference
891 // criterion, and returns whether it was changed.
892 const auto update_neighbor_level =
893 [&future_levels, max_difference, invalid_level](
894 const auto &neighbor, const level_type cell_level) -> bool {
895 Assert(neighbor->is_active(), ExcInternalError());
896 // We only care about locally owned neighbors. If neighbor is a ghost
897 // cell, its future FE index will be updated on the owning process and
898 // communicated at the next loop iteration.
899 if (neighbor->is_locally_owned())
900 {
901 const level_type neighbor_level = static_cast<level_type>(
902 future_levels[neighbor->global_active_cell_index()]);
903
904 // ignore neighbors that are not part of the hierarchy
905 if (neighbor_level == invalid_level)
906 return false;
907
908 if ((cell_level - max_difference) > neighbor_level)
909 {
910 future_levels[neighbor->global_active_cell_index()] =
911 cell_level - max_difference;
912
913 return true;
914 }
915 }
916
917 return false;
918 };
919
920 // For cells to be h-coarsened, we need to determine a future FE for the
921 // parent cell, which will be the dominated FE among all children
922 // However, if we want to enforce the max_difference criterion on all
923 // cells on the updated mesh, we will need to simulate the updated mesh on
924 // the current mesh.
925 //
926 // As we are working on p-levels, we will set all siblings that will be
927 // coarsened to the highest p-level among them. The parent cell will be
928 // assigned exactly this level in form of the corresponding FE index in
929 // the adaptation process in
930 // Triangulation::execute_coarsening_and_refinement().
931 //
932 // This function takes a cell and sets all its siblings to the highest
933 // p-level among them. Returns whether any future levels have been
934 // changed.
935 const auto prepare_level_for_parent = [&](const auto &neighbor) -> bool {
936 Assert(neighbor->is_active(), ExcInternalError());
937 if (neighbor->coarsen_flag_set() && neighbor->is_locally_owned())
938 {
939 const auto parent = neighbor->parent();
940
941 std::vector<unsigned int> future_levels_children;
942 future_levels_children.reserve(parent->n_children());
943 for (const auto &child : parent->child_iterators())
944 {
945 Assert(child->is_active() && child->coarsen_flag_set(),
947 ExcInconsistentCoarseningFlags()));
948
949 const level_type child_level = static_cast<level_type>(
950 future_levels[child->global_active_cell_index()]);
951 Assert(child_level != invalid_level,
953 "The FiniteElement on one of the siblings of "
954 "a cell you are trying to coarsen is not part "
955 "of the registered p-adaptation hierarchy."));
956 future_levels_children.push_back(child_level);
957 }
958 Assert(!future_levels_children.empty(), ExcInternalError());
959
960 const unsigned int max_level_children =
961 *std::max_element(future_levels_children.begin(),
962 future_levels_children.end());
963
964 bool children_changed = false;
965 for (const auto &child : parent->child_iterators())
966 // We only care about locally owned children. If child is a ghost
967 // cell, its future FE index will be updated on the owning process
968 // and communicated at the next loop iteration.
969 if (child->is_locally_owned() &&
970 future_levels[child->global_active_cell_index()] !=
971 max_level_children)
972 {
973 future_levels[child->global_active_cell_index()] =
974 max_level_children;
975 children_changed = true;
976 }
977 return children_changed;
978 }
979
980 return false;
981 };
982
983 bool levels_changed = false;
984 bool levels_changed_in_cycle;
985 do
986 {
987 levels_changed_in_cycle = false;
988
989 future_levels.update_ghost_values();
990
991 for (const auto &cell : dof_handler.active_cell_iterators())
992 if (!cell->is_artificial())
993 {
994 const level_type cell_level = static_cast<level_type>(
995 future_levels[cell->global_active_cell_index()]);
996
997 // ignore cells that are not part of the hierarchy
998 if (cell_level == invalid_level)
999 continue;
1000
1001 // ignore lowest levels of the hierarchy that always fulfill the
1002 // max_difference criterion
1003 if (cell_level <= max_difference)
1004 continue;
1005
1006 for (unsigned int f = 0; f < cell->n_faces(); ++f)
1007 if (cell->face(f)->at_boundary() == false)
1008 {
1009 if (cell->face(f)->has_children())
1010 {
1011 for (unsigned int sf = 0;
1012 sf < cell->face(f)->n_children();
1013 ++sf)
1014 {
1015 const auto neighbor =
1016 cell->neighbor_child_on_subface(f, sf);
1017
1018 levels_changed_in_cycle |=
1019 update_neighbor_level(neighbor, cell_level);
1020
1021 levels_changed_in_cycle |=
1022 prepare_level_for_parent(neighbor);
1023 }
1024 }
1025 else
1026 {
1027 const auto neighbor = cell->neighbor(f);
1028
1029 levels_changed_in_cycle |=
1030 update_neighbor_level(neighbor, cell_level);
1031
1032 levels_changed_in_cycle |=
1033 prepare_level_for_parent(neighbor);
1034 }
1035 }
1036 }
1037
1038 levels_changed_in_cycle =
1039 Utilities::MPI::logical_or(levels_changed_in_cycle,
1040 dof_handler.get_communicator());
1041 levels_changed |= levels_changed_in_cycle;
1042 }
1043 while (levels_changed_in_cycle);
1044
1045 // update future FE indices on locally owned cells
1046 for (const auto &cell : dof_handler.active_cell_iterators() |
1048 {
1049 const level_type cell_level = static_cast<level_type>(
1050 future_levels[cell->global_active_cell_index()]);
1051
1052 if (cell_level != invalid_level)
1053 {
1054 const unsigned int fe_index =
1055 fe_index_for_hierarchy_level[cell_level];
1056
1057 // only update if necessary
1058 if (fe_index != cell->active_fe_index())
1059 cell->set_future_fe_index(fe_index);
1060 }
1061 }
1062
1063 return levels_changed;
1064 }
1065 } // namespace Refinement
1066} // namespace hp
1067
1068
1069// explicit instantiations
1070#include "refinement.inst"
1071
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
iterator end()
size_type size() const
void grow_or_shrink(const size_type N)
iterator begin()
unsigned int size() const
Definition collection.h:265
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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)
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:118
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)
static const unsigned int invalid_unsigned_int
Definition types.h:213
typename type_identity< T >::type type_identity_t
Definition type_traits.h:96
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
unsigned short int fe_index
Definition types.h:60