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