Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
24 
25 #include <deal.II/dofs/dof_accessor.templates.h>
27 
30 
31 #include <deal.II/hp/refinement.h>
32 
34 #include <deal.II/lac/vector.h>
35 
37 
38 namespace 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
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
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,
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
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,
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),
170  Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
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 &&
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 
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
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,
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),
267  Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
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 &&
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 =
356  compute_global_min_and_max_at_root(indicators_refinement,
357  mpi_communicator);
358 
359  const std::pair<Number, Number> global_min_max_coarsening =
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::
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::
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.
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
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
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,
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,
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(),
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)
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
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
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 
static const unsigned int invalid_unsigned_int
Definition: types.h:196
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)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
void p_adaptivity_from_regularity(const ::DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &sobolev_indices)
Definition: refinement.cc:448
void p_adaptivity_from_flags(const ::DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &p_flags)
Definition: refinement.cc:67
T logical_or(const T &t, const MPI_Comm &mpi_communicator)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
void reinit(const size_type size, const bool omit_zeroing_entries=false)
void choose_p_over_h(const ::DoFHandler< dim, spacedim > &dof_handler)
Definition: refinement.cc:690
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
virtual MPI_Comm get_communicator() const override
Definition: tria_base.cc:140
static ::ExceptionBase & ExcMessage(std::string arg1)
void full_p_adaptivity(const ::DoFHandler< dim, spacedim > &dof_handler)
Definition: refinement.cc:47
Expression ceil(const Expression &x)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1465
void force_p_over_h(const ::DoFHandler< dim, spacedim > &dof_handler)
Definition: refinement.cc:668
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
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
TriangulationBase< dim, spacedim > Triangulation
Definition: tria_base.h:396
Definition: hp.h:117
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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
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
T min(const T &t, const MPI_Comm &mpi_communicator)
std::function< bool(const Number &, const Number &)> ComparisonFunction
Definition: refinement.h:139
std::pair< number, number > compute_global_min_and_max_at_root(const ::Vector< number > &criteria, const MPI_Comm &mpi_communicator)
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
static ::ExceptionBase & ExcInvalidParameterValue()
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
T max(const T &t, const MPI_Comm &mpi_communicator)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
Expression floor(const Expression &x)