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