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