Reference documentation for deal.II version Git 73c87d96ef 2021-11-30 22:54:44 +0100
\(\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 
31 
32 #include <deal.II/hp/refinement.h>
33 
35 #include <deal.II/lac/vector.h>
36 
38 
39 namespace 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
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
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,
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
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,
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),
171  Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
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 &&
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 
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
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,
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),
268  Assert((p_coarsen_fraction >= 0) && (p_coarsen_fraction <= 1),
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 &&
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 =
357  compute_global_min_and_max_at_root(indicators_refinement,
358  mpi_communicator);
359 
360  const std::pair<Number, Number> global_min_max_coarsening =
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::
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::
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.
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
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
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,
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,
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  if (cell->is_locally_owned())
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]
624  .degree;
625  }
626  else
627  {
628  // future finite element on current cell is already set
629  future_fe_degree =
630  dof_handler.get_fe_collection()[cell->future_fe_index()]
631  .degree;
632  }
633 
634  // step 1: exponential decay with p-adaptation
635  if (cell->future_fe_index_set())
636  {
637  predicted_errors[cell->active_cell_index()] *=
638  std::pow(gamma_p,
639  int(future_fe_degree) - int(cell->get_fe().degree));
640  }
641 
642  // step 2: algebraic decay with h-adaptation
643  if (cell->refine_flag_set())
644  {
645  predicted_errors[cell->active_cell_index()] *=
646  (gamma_h * std::pow(.5, future_fe_degree));
647 
648  // predicted error will be split on children cells
649  // after adaptation via CellDataTransfer
650  }
651  else if (cell->coarsen_flag_set())
652  {
653  predicted_errors[cell->active_cell_index()] /=
654  (gamma_h * std::pow(.5, future_fe_degree));
655 
656  // predicted error will be summed up on parent cell
657  // after adaptation via CellDataTransfer
658  }
659  }
660  }
661 
662 
663 
667  template <int dim, int spacedim>
668  void
669  force_p_over_h(const ::DoFHandler<dim, spacedim> &dof_handler)
670  {
671  if (dof_handler.get_fe_collection().size() == 0)
672  // nothing to do
673  return;
674 
675  Assert(
676  dof_handler.has_hp_capabilities(),
677  (typename ::DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
678 
679  for (const auto &cell : dof_handler.active_cell_iterators())
680  if (cell->is_locally_owned() && cell->future_fe_index_set())
681  {
682  cell->clear_refine_flag();
683  cell->clear_coarsen_flag();
684  }
685  }
686 
687 
688 
689  template <int dim, int spacedim>
690  void
691  choose_p_over_h(const ::DoFHandler<dim, spacedim> &dof_handler)
692  {
693  if (dof_handler.get_fe_collection().size() == 0)
694  // nothing to do
695  return;
696 
697  Assert(
698  dof_handler.has_hp_capabilities(),
699  (typename ::DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
700 
701  // Ghost siblings might occur on parallel::shared::Triangulation objects.
702  // We need information about future FE indices on all locally relevant
703  // cells here, and thus communicate them.
704  if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
705  &dof_handler.get_triangulation()) != nullptr)
708  const_cast<::DoFHandler<dim, spacedim> &>(dof_handler));
709 
710  for (const auto &cell : dof_handler.active_cell_iterators())
711  if (cell->is_locally_owned() && cell->future_fe_index_set())
712  {
713  cell->clear_refine_flag();
714 
715  // A cell will only be coarsened into its parent if all of its
716  // siblings are flagged for h-coarsening as well. We must take this
717  // into account for our decision whether we would like to impose h-
718  // or p-adaptivity.
719  if (cell->coarsen_flag_set())
720  {
721  const auto & parent = cell->parent();
722  const unsigned int n_children = parent->n_children();
723 
724  unsigned int h_flagged_children = 0, p_flagged_children = 0;
725  for (const auto &child : parent->child_iterators())
726  {
727  if (child->is_active())
728  {
729  if (child->is_locally_owned())
730  {
731  if (child->coarsen_flag_set())
732  ++h_flagged_children;
733  if (child->future_fe_index_set())
734  ++p_flagged_children;
735  }
736  else if (child->is_ghost())
737  {
738  // The case of siblings being owned by different
739  // processors can only occur for
740  // parallel::shared::Triangulation objects.
741  Assert(
742  (dynamic_cast<const parallel::shared::
743  Triangulation<dim, spacedim> *>(
744  &dof_handler.get_triangulation()) != nullptr),
745  ExcInternalError());
746 
747  if (child->coarsen_flag_set())
748  ++h_flagged_children;
749  // The public interface does not allow to access
750  // future FE indices on ghost cells. However, we
751  // need this information here and thus call the
752  // internal function that does not check for cell
753  // ownership.
754  if (::internal::
755  DoFCellAccessorImplementation::
756  Implementation::
757  future_fe_index_set<dim, spacedim, false>(
758  *child))
759  ++p_flagged_children;
760  }
761  else
762  {
763  // Siblings of locally owned cells are all
764  // either also locally owned or ghost cells.
765  Assert(false, ExcInternalError());
766  }
767  }
768  }
769 
770  if (h_flagged_children == n_children &&
771  p_flagged_children != n_children)
772  {
773  // Perform pure h-coarsening and
774  // drop all p-adaptation flags.
775  for (const auto &child : parent->child_iterators())
776  {
777  // h_flagged_children == n_children implies
778  // that all children are active
779  Assert(child->is_active(), ExcInternalError());
780  if (child->is_locally_owned())
781  child->clear_future_fe_index();
782  }
783  }
784  else
785  {
786  // Perform p-adaptation on all children and
787  // drop all h-coarsening flags.
788  for (const auto &child : parent->child_iterators())
789  {
790  if (child->is_active() && child->is_locally_owned())
791  child->clear_coarsen_flag();
792  }
793  }
794  }
795  }
796  }
797 
798 
799 
803  template <int dim, int spacedim>
804  bool
806  const ::DoFHandler<dim, spacedim> &dof_handler,
807  const unsigned int max_difference,
808  const unsigned int contains_fe_index)
809  {
810  if (dof_handler.get_fe_collection().size() == 0)
811  // nothing to do
812  return false;
813 
814  Assert(
815  dof_handler.has_hp_capabilities(),
816  (typename ::DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
817  Assert(
818  max_difference > 0,
819  ExcMessage(
820  "This function does not serve any purpose for max_difference = 0."));
821  AssertIndexRange(contains_fe_index,
822  dof_handler.get_fe_collection().size());
823 
824  //
825  // establish hierarchy
826  //
827  // - create bimap between hierarchy levels and FE indices
828 
829  // there can be as many levels in the hierarchy as active FE indices are
830  // possible
831  using level_type =
832  typename ::DoFHandler<dim, spacedim>::active_fe_index_type;
833  const auto invalid_level = static_cast<level_type>(-1);
834 
835  // map from FE index to level in hierarchy
836  // FE indices that are not covered in the hierarchy are not in the map
837  const std::vector<unsigned int> fe_index_for_hierarchy_level =
838  dof_handler.get_fe_collection().get_hierarchy_sequence(
839  contains_fe_index);
840 
841  // map from level in hierarchy to FE index
842  // FE indices that are not covered in the hierarchy will be mapped to
843  // invalid_level
844  std::vector<unsigned int> hierarchy_level_for_fe_index(
845  dof_handler.get_fe_collection().size(), invalid_level);
846  for (unsigned int l = 0; l < fe_index_for_hierarchy_level.size(); ++l)
847  hierarchy_level_for_fe_index[fe_index_for_hierarchy_level[l]] = l;
848 
849 
850  //
851  // parallelization
852  //
853  // - create distributed vector of level indices
854  // - update ghost values in each iteration (see later)
855  // - no need to compress, since the owning processor will have the correct
856  // level index
857 
858  // HOTFIX: ::Vector does not accept integral types
860  if (const auto parallel_tria =
861  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
862  &(dof_handler.get_triangulation())))
863  {
864  future_levels.reinit(
865  parallel_tria->global_active_cell_index_partitioner().lock());
866  }
867  else
868  {
869  future_levels.reinit(
870  dof_handler.get_triangulation().n_active_cells());
871  }
872 
873  for (const auto &cell : dof_handler.active_cell_iterators() |
875  future_levels[cell->global_active_cell_index()] =
876  hierarchy_level_for_fe_index[cell->future_fe_index()];
877 
878 
879  //
880  // limit level difference of neighboring cells
881  //
882  // - go over all locally relevant cells, and adjust the level indices of
883  // locally owned neighbors to match the level difference (as a
884  // consequence, indices on ghost cells will be updated only on the
885  // owning processor)
886  // - always raise levels to match criterion, never lower them
887  // - exchange level indices on ghost cells
888 
889  // Function that updates the level of neighbor to fulfill difference
890  // criterion, and returns whether it was changed.
891  const auto update_neighbor_level =
892  [&future_levels, max_difference, invalid_level](
893  const auto &neighbor, const level_type cell_level) -> bool {
894  Assert(neighbor->is_active(), ExcInternalError());
895  // We only care about locally owned neighbors. If neighbor is a ghost
896  // cell, its future FE index will be updated on the owning process and
897  // communicated at the next loop iteration.
898  if (neighbor->is_locally_owned())
899  {
900  const level_type neighbor_level = static_cast<level_type>(
901  future_levels[neighbor->global_active_cell_index()]);
902 
903  // ignore neighbors that are not part of the hierarchy
904  if (neighbor_level == invalid_level)
905  return false;
906 
907  if ((cell_level - max_difference) > neighbor_level)
908  {
909  future_levels[neighbor->global_active_cell_index()] =
910  cell_level - max_difference;
911 
912  return true;
913  }
914  }
915 
916  return false;
917  };
918 
919  // For cells to be h-coarsened, we need to determine a future FE for the
920  // parent cell, which will be the dominated FE among all children
921  // However, if we want to enforce the max_difference criterion on all
922  // cells on the updated mesh, we will need to simulate the updated mesh on
923  // the current mesh.
924  //
925  // As we are working on p-levels, we will set all siblings that will be
926  // coarsened to the highest p-level among them. The parent cell will be
927  // assigned exactly this level in form of the corresponding FE index in
928  // the adaptation process in
929  // Triangulation::execute_coarsening_and_refinement().
930  //
931  // This function takes a cell and sets all its siblings to the highest
932  // p-level among them. Returns whether any future levels have been
933  // changed.
934  const auto prepare_level_for_parent = [&](const auto &neighbor) -> bool {
935  Assert(neighbor->is_active(), ExcInternalError());
936  if (neighbor->coarsen_flag_set() && neighbor->is_locally_owned())
937  {
938  const auto parent = neighbor->parent();
939 
940  std::vector<unsigned int> future_levels_children;
941  future_levels_children.reserve(parent->n_children());
942  for (const auto &child : parent->child_iterators())
943  {
944  Assert(child->is_active() && child->coarsen_flag_set(),
945  (typename ::Triangulation<dim, spacedim>::
946  ExcInconsistentCoarseningFlags()));
947 
948  const level_type child_level = static_cast<level_type>(
949  future_levels[child->global_active_cell_index()]);
950  Assert(child_level != invalid_level,
951  ExcMessage(
952  "The FiniteElement on one of the siblings of "
953  "a cell you are trying to coarsen is not part "
954  "of the registered p-adaptation hierarchy."));
955  future_levels_children.push_back(child_level);
956  }
957  Assert(!future_levels_children.empty(), ExcInternalError());
958 
959  const unsigned int max_level_children =
960  *std::max_element(future_levels_children.begin(),
961  future_levels_children.end());
962 
963  bool children_changed = false;
964  for (const auto &child : parent->child_iterators())
965  // We only care about locally owned children. If child is a ghost
966  // cell, its future FE index will be updated on the owning process
967  // and communicated at the next loop iteration.
968  if (child->is_locally_owned() &&
969  future_levels[child->global_active_cell_index()] !=
970  max_level_children)
971  {
972  future_levels[child->global_active_cell_index()] =
973  max_level_children;
974  children_changed = true;
975  }
976  return children_changed;
977  }
978 
979  return false;
980  };
981 
982  bool levels_changed = false;
983  bool levels_changed_in_cycle;
984  do
985  {
986  levels_changed_in_cycle = false;
987 
988  future_levels.update_ghost_values();
989 
990  for (const auto &cell : dof_handler.active_cell_iterators())
991  if (!cell->is_artificial())
992  {
993  const level_type cell_level = static_cast<level_type>(
994  future_levels[cell->global_active_cell_index()]);
995 
996  // ignore cells that are not part of the hierarchy
997  if (cell_level == invalid_level)
998  continue;
999 
1000  // ignore lowest levels of the hierarchy that always fulfill the
1001  // max_difference criterion
1002  if (cell_level <= max_difference)
1003  continue;
1004 
1005  for (unsigned int f = 0; f < cell->n_faces(); ++f)
1006  if (cell->face(f)->at_boundary() == false)
1007  {
1008  if (cell->face(f)->has_children())
1009  {
1010  for (unsigned int sf = 0;
1011  sf < cell->face(f)->n_children();
1012  ++sf)
1013  {
1014  const auto neighbor =
1015  cell->neighbor_child_on_subface(f, sf);
1016 
1017  levels_changed_in_cycle |=
1018  update_neighbor_level(neighbor, cell_level);
1019 
1020  levels_changed_in_cycle |=
1021  prepare_level_for_parent(neighbor);
1022  }
1023  }
1024  else
1025  {
1026  const auto neighbor = cell->neighbor(f);
1027 
1028  levels_changed_in_cycle |=
1029  update_neighbor_level(neighbor, cell_level);
1030 
1031  levels_changed_in_cycle |=
1032  prepare_level_for_parent(neighbor);
1033  }
1034  }
1035  }
1036 
1037  levels_changed_in_cycle =
1038  Utilities::MPI::logical_or(levels_changed_in_cycle,
1039  dof_handler.get_communicator());
1040  levels_changed |= levels_changed_in_cycle;
1041  }
1042  while (levels_changed_in_cycle);
1043 
1044  // update future FE indices on locally owned cells
1045  for (const auto &cell : dof_handler.active_cell_iterators() |
1047  {
1048  const level_type cell_level = static_cast<level_type>(
1049  future_levels[cell->global_active_cell_index()]);
1050 
1051  if (cell_level != invalid_level)
1052  {
1053  const unsigned int fe_index =
1054  fe_index_for_hierarchy_level[cell_level];
1055 
1056  // only update if necessary
1057  if (fe_index != cell->active_fe_index())
1058  cell->set_future_fe_index(fe_index);
1059  }
1060  }
1061 
1062  return levels_changed;
1063  }
1064  } // namespace Refinement
1065 } // namespace hp
1066 
1067 
1068 // explicit instantiations
1069 #include "refinement.inst"
1070 
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:1655
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:449
void p_adaptivity_from_flags(const ::DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &p_flags)
Definition: refinement.cc:68
T logical_or(const T &t, const MPI_Comm &mpi_communicator)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
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:691
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
virtual MPI_Comm get_communicator() const override
Definition: tria_base.cc:142
static ::ExceptionBase & ExcMessage(std::string arg1)
void full_p_adaptivity(const ::DoFHandler< dim, spacedim > &dof_handler)
Definition: refinement.cc:48
Expression ceil(const Expression &x)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1461
void force_p_over_h(const ::DoFHandler< dim, spacedim > &dof_handler)
Definition: refinement.cc:669
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
::DoFHandler< dim, spacedim > DoFHandler
Definition: dof_handler.h:124
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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:151
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
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
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:805
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:248
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)