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