Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18:55:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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() |
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)
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
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
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() |
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() |
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 
Definition: vector.h:109
virtual MPI_Comm get_communicator() const override
Definition: tria_base.cc:144
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidParameterValue()
iterator end()
size_type size() const
void grow_or_shrink(const size_type N)
void reinit(const size_type size, const bool omit_zeroing_entries=false)
iterator begin()
Expression ceil(const Expression &x)
Expression floor(const Expression &x)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T logical_or(const T &t, const MPI_Comm &mpi_communicator)
T min(const T &t, const MPI_Comm &mpi_communicator)
T sum(const T &t, const MPI_Comm &mpi_communicator)
T max(const T &t, const MPI_Comm &mpi_communicator)
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
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
void choose_p_over_h(const ::DoFHandler< dim, spacedim > &dof_handler)
Definition: refinement.cc:689
void p_adaptivity_from_flags(const ::DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &p_flags)
Definition: refinement.cc:68
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
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
void force_p_over_h(const ::DoFHandler< dim, spacedim > &dof_handler)
Definition: refinement.cc:667
void full_p_adaptivity(const ::DoFHandler< dim, spacedim > &dof_handler)
Definition: refinement.cc:48
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
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
std::function< bool(const Number &, const Number &)> ComparisonFunction
Definition: refinement.h:139
void p_adaptivity_from_regularity(const ::DoFHandler< dim, spacedim > &dof_handler, const Vector< Number > &sobolev_indices)
Definition: refinement.cc:449
Definition: hp.h:118
void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
number compute_threshold(const ::Vector< number > &criteria, const std::pair< double, double > &global_min_and_max, const types::global_cell_index n_target_cells, const MPI_Comm &mpi_communicator)
std::pair< number, number > compute_global_min_and_max_at_root(const ::Vector< number > &criteria, const MPI_Comm &mpi_communicator)
static const unsigned int invalid_unsigned_int
Definition: types.h:201