Reference documentation for deal.II version Git 1dc1051882 2021-04-22 23:57:03 +0200
\(\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\}}\)
dof_handler_policy.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 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 
21 #include <deal.II/base/utilities.h>
23 
26 
30 
31 #include <deal.II/fe/fe.h>
32 
34 #include <deal.II/grid/tria.h>
36 
37 #include <algorithm>
38 #include <memory>
39 #include <numeric>
40 #include <set>
41 
43 
44 
45 namespace internal
46 {
47  namespace DoFHandlerImplementation
48  {
49  namespace Policy
50  {
51  // use class ::DoFHandler instead
52  // of namespace internal::DoFHandler in
53  // the following
54  using ::DoFHandler;
55 
56  namespace
57  {
64  const types::global_dof_index enumeration_dof_index =
66 
71  template <int dim, int spacedim>
72  void
73  update_all_active_cell_dof_indices_caches(
74  const DoFHandler<dim, spacedim> &dof_handler)
75  {
76  const auto worker = [](const auto &cell, void *, void *) {
77  if (!cell->is_artificial())
78  cell->update_cell_dof_indices_cache();
79  };
80 
81  // parallelize filling all of the cell caches. by using
82  // WorkStream, we make sure that we only run through the
83  // range of iterators once, whereas a parallel_for loop
84  // for example has to split the range multiple times,
85  // which is expensive because cell iterators are not
86  // random access iterators with a cheap operator-
87  WorkStream::run(dof_handler.begin_active(),
88  dof_handler.end(),
89  worker,
90  /* copier */ std::function<void(void *)>(),
91  /* scratch_data */ nullptr,
92  /* copy_data */ nullptr,
94  /* chunk_size = */ 32);
95  }
96 
97 
102  template <int dim, int spacedim>
103  void
104  update_all_level_cell_dof_indices_caches(
105  const DoFHandler<dim, spacedim> &dof_handler)
106  {
107  const auto worker = [](const auto &cell, void *, void *) {
108  if (cell->has_children() || !cell->is_artificial())
109  cell->update_cell_dof_indices_cache();
110  };
111 
112  // parallelize filling all of the cell caches. by using
113  // WorkStream, we make sure that we only run through the
114  // range of iterators once, whereas a parallel_for loop
115  // for example has to split the range multiple times,
116  // which is expensive because cell iterators are not
117  // random access iterators with a cheap operator-
118  WorkStream::run(dof_handler.begin(),
119  dof_handler.end(),
120  worker,
121  /* copier */ std::function<void(void *)>(),
122  /* scratch_data */ nullptr,
123  /* copy_data */ nullptr,
125  /* chunk_size = */ 32);
126  }
127 
128 
129  using DoFIdentities =
130  std::vector<std::pair<unsigned int, unsigned int>>;
131 
132 
143  template <int structdim, int dim, int spacedim>
144  const std::unique_ptr<DoFIdentities> &
145  ensure_existence_and_return_dof_identities(
146  const FiniteElement<dim, spacedim> &fe1,
147  const FiniteElement<dim, spacedim> &fe2,
148  std::unique_ptr<DoFIdentities> & identities,
149  const unsigned int face_no = numbers::invalid_unsigned_int)
150  {
151  Assert(structdim == 2 || face_no == numbers::invalid_unsigned_int,
152  ExcInternalError());
153 
154  // see if we need to fill this entry, or whether it already
155  // exists
156  if (identities.get() == nullptr)
157  {
158  switch (structdim)
159  {
160  case 0:
161  {
162  identities = std::make_unique<DoFIdentities>(
163  fe1.hp_vertex_dof_identities(fe2));
164  break;
165  }
166 
167  case 1:
168  {
169  identities = std::make_unique<DoFIdentities>(
170  fe1.hp_line_dof_identities(fe2));
171  break;
172  }
173 
174  case 2:
175  {
176  identities = std::make_unique<DoFIdentities>(
177  fe1.hp_quad_dof_identities(fe2, face_no));
178  break;
179  }
180 
181  default:
182  Assert(false, ExcNotImplemented());
183  }
184 
185  // double check whether the newly created entries make
186  // any sense at all
187  for (unsigned int i = 0; i < identities->size(); ++i)
188  {
189  Assert((*identities)[i].first <
190  fe1.template n_dofs_per_object<structdim>(face_no),
191  ExcInternalError());
192  Assert((*identities)[i].second <
193  fe2.template n_dofs_per_object<structdim>(face_no),
194  ExcInternalError());
195  }
196  }
197 
198  return identities;
199  }
200  } // namespace
201 
202 
203 
205  {
206  /* -------------- distribute_dofs functionality ------------- */
207 
212  template <int dim, int spacedim>
213  static std::map<types::global_dof_index, types::global_dof_index>
215  const DoFHandler<dim, spacedim> &dof_handler)
216  {
217  Assert(
218  dof_handler.hp_capability_enabled == true,
220 
221  std::map<types::global_dof_index, types::global_dof_index>
222  dof_identities;
223 
224  // Note: we may wish to have something here similar to what
225  // we do for lines and quads, namely that we only identify
226  // dofs for any FE towards the most dominating one. however,
227  // it is not clear whether this is actually necessary for
228  // vertices at all, I can't think of a finite element that
229  // would make that necessary...
231  vertex_dof_identities(dof_handler.get_fe_collection().size(),
232  dof_handler.get_fe_collection().size());
233 
234  // loop over all vertices and see which one we need to work on
235  for (unsigned int vertex_index = 0;
236  vertex_index < dof_handler.get_triangulation().n_vertices();
237  ++vertex_index)
238  if (dof_handler.get_triangulation()
239  .get_used_vertices()[vertex_index] == true)
240  {
241  const unsigned int n_active_fe_indices =
242  ::internal::DoFAccessorImplementation::Implementation::
243  n_active_fe_indices(dof_handler,
244  0,
245  vertex_index,
246  std::integral_constant<int, 0>());
247 
248  if (n_active_fe_indices > 1)
249  {
250  const std::set<unsigned int> fe_indices =
253  dof_handler,
254  0,
255  vertex_index,
256  std::integral_constant<int, 0>());
257 
258  // find out which is the most dominating finite
259  // element of the ones that are used on this vertex
260  unsigned int most_dominating_fe_index =
261  dof_handler.get_fe_collection().find_dominating_fe(
262  fe_indices,
263  /*codim*/ dim);
264 
265  // if we haven't found a dominating finite element,
266  // choose the very first one to be dominant
267  if (most_dominating_fe_index ==
269  most_dominating_fe_index =
270  ::internal::DoFAccessorImplementation::
271  Implementation::nth_active_fe_index(
272  dof_handler,
273  0,
274  vertex_index,
275  0,
276  std::integral_constant<int, 0>());
277 
278  // loop over the indices of all the finite
279  // elements that are not dominating, and
280  // identify their dofs to the most dominating
281  // one
282  for (const auto &other_fe_index : fe_indices)
283  if (other_fe_index != most_dominating_fe_index)
284  {
285  // make sure the entry in the equivalence
286  // table exists
287  const auto &identities =
288  *ensure_existence_and_return_dof_identities<0>(
289  dof_handler.get_fe(most_dominating_fe_index),
290  dof_handler.get_fe(other_fe_index),
291  vertex_dof_identities[most_dominating_fe_index]
292  [other_fe_index]);
293 
294  // then loop through the identities we
295  // have. first get the global numbers of the
296  // dofs we want to identify and make sure they
297  // are not yet constrained to anything else,
298  // except for to each other. use the rule that
299  // we will always constrain the dof with the
300  // higher FE index to the one with the lower,
301  // to avoid circular reasoning.
302  for (const auto &identity : identities)
303  {
304  const types::global_dof_index primary_dof_index =
305  ::internal::DoFAccessorImplementation::
306  Implementation::get_dof_index(
307  dof_handler,
308  0,
309  vertex_index,
310  most_dominating_fe_index,
311  identity.first,
312  std::integral_constant<int, 0>());
314  dependent_dof_index =
315  ::internal::DoFAccessorImplementation::
316  Implementation::get_dof_index(
317  dof_handler,
318  0,
319  vertex_index,
320  other_fe_index,
321  identity.second,
322  std::integral_constant<int, 0>());
323 
324  // on subdomain boundaries, we will
325  // encounter invalid DoFs on ghost cells,
326  // for which we have not yet distributed
327  // valid indices. depending on which finte
328  // element is dominating the other on this
329  // interface, we either have to constrain
330  // the valid to the invalid indices, or vice
331  // versa.
332  //
333  // we only store an identity if we are about
334  // to overwrite a valid DoF. we will skip
335  // constraining invalid DoFs for now, and
336  // consider them later in Phase 5.
337  if (dependent_dof_index !=
339  {
340  // if the DoF indices of both elements
341  // are already distributed, i.e., both
342  // of these 'fe_indices' are associated
343  // with a locally owned cell, then we
344  // should either not have a dof_identity
345  // yet, or it must come out here to be
346  // exactly as we had computed before
347  if (primary_dof_index !=
349  Assert(
350  (dof_identities.find(primary_dof_index) ==
351  dof_identities.end()) ||
352  (dof_identities[dependent_dof_index] ==
353  primary_dof_index),
354  ExcInternalError());
355 
356  dof_identities[dependent_dof_index] =
357  primary_dof_index;
358  }
359  }
360  }
361  }
362  }
363 
364  return dof_identities;
365  }
366 
367 
372  template <int spacedim>
373  static std::map<types::global_dof_index, types::global_dof_index>
375  {
376  (void)dof_handler;
377  Assert(dof_handler.hp_capability_enabled == true,
379 
380  return std::map<types::global_dof_index, types::global_dof_index>();
381  }
382 
383 
384  template <int dim, int spacedim>
385  static std::map<types::global_dof_index, types::global_dof_index>
387  const DoFHandler<dim, spacedim> &dof_handler)
388  {
389  Assert(
390  dof_handler.hp_capability_enabled == true,
392 
393  std::map<types::global_dof_index, types::global_dof_index>
394  dof_identities;
395 
396  // we will mark lines that we have already treated, so first save and
397  // clear the user flags on lines and later restore them
398  std::vector<bool> user_flags;
399  dof_handler.get_triangulation().save_user_flags_line(user_flags);
400  const_cast<::Triangulation<dim, spacedim> &>(
401  dof_handler.get_triangulation())
402  .clear_user_flags_line();
403 
404  // An implementation of the algorithm described in the hp-paper,
405  // including the modification mentioned later in the "complications in
406  // 3-d" subsections
407  //
408  // as explained there, we do something only if there are exactly 2
409  // finite elements associated with an object. if there is only one,
410  // then there is nothing to do anyway, and if there are 3 or more,
411  // then we can get into trouble. note that this only happens for lines
412  // in 3d and higher, and for quads only in 4d and higher, so this
413  // isn't a particularly frequent case
414  //
415  // there is one case, however, that we would like to handle (see, for
416  // example, the hp/crash_15 testcase): if we have
417  // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
418  // then we should be able to handle this because we can simply unify
419  // *all* dofs, not only a some. so what we do is to first treat all
420  // pairs of finite elements that have *identical* dofs, and then only
421  // deal with those that are not identical of which we can handle at
422  // most 2
423  ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
424  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
425 
426  for (const auto &cell : dof_handler.active_cell_iterators())
427  for (const auto l : cell->line_indices())
428  if (cell->line(l)->user_flag_set() == false)
429  {
430  const auto line = cell->line(l);
431  line->set_user_flag();
432 
433  unsigned int unique_sets_of_dofs =
434  line->n_active_fe_indices();
435 
436  // do a first loop over all sets of dofs and do identity
437  // uniquification
438  const unsigned int n_active_fe_indices =
439  line->n_active_fe_indices();
440  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
441  for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
442  {
443  const unsigned int fe_index_1 =
444  line->nth_active_fe_index(f),
445  fe_index_2 =
446  line->nth_active_fe_index(g);
447 
448  // as described in the hp-paper, we only unify on lines
449  // when there are at most two different FE objects
450  // assigned on it.
451  // however, more than two 'active_fe_indices' can be
452  // attached that still fulfill the above criterion,
453  // i.e. when two different FiniteElement objects are
454  // assigned to neighboring cells that map their degrees
455  // of freedom one-to-one.
456  // we cannot verify with certainty if two dofs each of
457  // separate FiniteElement objects actually map
458  // one-to-one. however, checking for the number of
459  // 'dofs_per_line' turned out to be a reasonable
460  // approach, that also works for e.g. two different
461  // FE_Q objects of the same order, from which one is
462  // enhanced by a bubble function that is zero on the
463  // boundary.
464  if ((dof_handler.get_fe(fe_index_1).n_dofs_per_line() ==
465  dof_handler.get_fe(fe_index_2)
466  .n_dofs_per_line()) &&
467  (dof_handler.get_fe(fe_index_1).n_dofs_per_line() >
468  0))
469  {
470  // the number of dofs per line is identical
471  const unsigned int dofs_per_line =
472  dof_handler.get_fe(fe_index_1).n_dofs_per_line();
473 
474  const auto &identities =
475  *ensure_existence_and_return_dof_identities<1>(
476  dof_handler.get_fe(fe_index_1),
477  dof_handler.get_fe(fe_index_2),
478  line_dof_identities[fe_index_1][fe_index_2]);
479  // see if these sets of dofs are identical. the
480  // first condition for this is that indeed there are
481  // n identities
482  if (identities.size() == dofs_per_line)
483  {
484  unsigned int i = 0;
485  for (; i < dofs_per_line; ++i)
486  if ((identities[i].first != i) &&
487  (identities[i].second != i))
488  // not an identity
489  break;
490 
491  if (i == dofs_per_line)
492  {
493  // The line dofs (i.e., the ones interior to
494  // a line) of these two finite elements are
495  // identical. Note that there could be
496  // situations when one element still
497  // dominates another, e.g.: FE_Q(2) x
498  // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
499 
500  --unique_sets_of_dofs;
501 
502  // determine which one of both finite
503  // elements is the dominating one.
504  const std::set<unsigned int> fe_indices{
505  fe_index_1, fe_index_2};
506 
507  unsigned int dominating_fe_index =
508  dof_handler.get_fe_collection()
509  .find_dominating_fe(fe_indices,
510  /*codim=*/dim - 1);
511  unsigned int other_fe_index =
513 
514  if (dominating_fe_index !=
516  other_fe_index =
517  (dominating_fe_index == fe_index_1) ?
518  fe_index_2 :
519  fe_index_1;
520  else
521  {
522  // if we haven't found a dominating
523  // finite element, choose the one with
524  // the lower index to be dominating
525  dominating_fe_index = fe_index_1;
526  other_fe_index = fe_index_2;
527  }
528 
529  for (unsigned int j = 0; j < dofs_per_line;
530  ++j)
531  {
533  primary_dof_index = line->dof_index(
534  j, dominating_fe_index);
536  dependent_dof_index =
537  line->dof_index(j, other_fe_index);
538 
539  // on subdomain boundaries, we will
540  // encounter invalid DoFs on ghost
541  // cells, for which we have not yet
542  // distributed valid indices. depending
543  // on which finte element is dominating
544  // the other on this interface, we
545  // either have to constrain the valid to
546  // the invalid indices, or vice versa.
547  //
548  // we only store an identity if we are
549  // about to overwrite a valid DoF. we
550  // will skip constraining invalid DoFs
551  // for now, and consider them later in
552  // Phase 5.
553  if (dependent_dof_index !=
555  {
556  if (primary_dof_index !=
558  {
559  // if primary dof was already
560  // constrained, constrain to
561  // that one, otherwise constrain
562  // dependent to primary
563  if (dof_identities.find(
564  primary_dof_index) !=
565  dof_identities.end())
566  {
567  // if the DoF indices of
568  // both elements are already
569  // distributed, i.e., both
570  // of these 'fe_indices' are
571  // associated with a locally
572  // owned cell, then we
573  // should either not have a
574  // dof_identity yet, or it
575  // must come out here to be
576  // exactly as we had
577  // computed before
578  Assert(
579  dof_identities.find(
580  dof_identities
581  [primary_dof_index]) ==
582  dof_identities.end(),
583  ExcInternalError());
584 
585  dof_identities
586  [dependent_dof_index] =
587  dof_identities
588  [primary_dof_index];
589  }
590  else
591  {
592  // see comment above for an
593  // explanation of this
594  // assertion
595  Assert(
596  (dof_identities.find(
597  primary_dof_index) ==
598  dof_identities.end()) ||
599  (dof_identities
600  [dependent_dof_index] ==
601  primary_dof_index),
602  ExcInternalError());
603 
604  dof_identities
605  [dependent_dof_index] =
606  primary_dof_index;
607  }
608  }
609  else
610  {
611  // set dependent_dof to
612  // primary_dof_index, which is
613  // invalid
614  dof_identities
615  [dependent_dof_index] =
617  }
618  }
619  }
620  }
621  }
622  }
623  }
624 
625  // if at this point, there is only one unique set of dofs
626  // left, then we have taken care of everything above. if there
627  // are two, then we need to deal with them here. if there are
628  // more, then we punt, as described in the paper (and
629  // mentioned above)
630  // TODO: The check for 'dim==2' was inserted by intuition. It
631  // fixes
632  // the previous problems with @ref step_27 "step-27" in 3D. But an
633  // explanation for this is still required, and what we do here
634  // is not what we describe in the paper!.
635  if ((unique_sets_of_dofs == 2) && (dim == 2))
636  {
637  const std::set<unsigned int> fe_indices =
638  line->get_active_fe_indices();
639 
640  // find out which is the most dominating finite element of
641  // the ones that are used on this line
642  const unsigned int most_dominating_fe_index =
643  dof_handler.get_fe_collection().find_dominating_fe(
644  fe_indices,
645  /*codim=*/dim - 1);
646 
647  // if we found the most dominating element, then use this
648  // to eliminate some of the degrees of freedom by
649  // identification. otherwise, the code that computes
650  // hanging node constraints will have to deal with it by
651  // computing appropriate constraints along this face/edge
652  if (most_dominating_fe_index !=
654  {
655  // loop over the indices of all the finite elements
656  // that are not dominating, and identify their dofs to
657  // the most dominating one
658  for (const auto &other_fe_index : fe_indices)
659  if (other_fe_index != most_dominating_fe_index)
660  {
661  const auto &identities =
662  *ensure_existence_and_return_dof_identities<
663  1>(dof_handler.get_fe(
664  most_dominating_fe_index),
665  dof_handler.get_fe(other_fe_index),
666  line_dof_identities
667  [most_dominating_fe_index]
668  [other_fe_index]);
669 
670  for (const auto &identity : identities)
671  {
673  primary_dof_index = line->dof_index(
674  identity.first,
675  most_dominating_fe_index);
677  dependent_dof_index =
678  line->dof_index(identity.second,
679  other_fe_index);
680 
681  // on subdomain boundaries, we will
682  // encounter invalid DoFs on ghost cells,
683  // for which we have not yet distributed
684  // valid indices. depending on which finte
685  // element is dominating the other on this
686  // interface, we either have to constrain
687  // the valid to the invalid indices, or vice
688  // versa.
689  //
690  // we only store an identity if we are about
691  // to overwrite a valid DoF. we will skip
692  // constraining invalid DoFs for now, and
693  // consider them later in Phase 5.
694  if (dependent_dof_index !=
696  {
697  // if the DoF indices of both elements
698  // are already distributed, i.e., both
699  // of these 'fe_indices' are associated
700  // with a locally owned cell, then we
701  // should either not have a dof_identity
702  // yet, or it must come out here to be
703  // exactly as we had computed before
704  if (primary_dof_index !=
706  Assert((dof_identities.find(
707  primary_dof_index) ==
708  dof_identities.end()) ||
709  (dof_identities
710  [dependent_dof_index] ==
711  primary_dof_index),
712  ExcInternalError());
713 
714  dof_identities[dependent_dof_index] =
715  primary_dof_index;
716  }
717  }
718  }
719  }
720  }
721  }
722 
723  // finally restore the user flags
724  const_cast<::Triangulation<dim, spacedim> &>(
725  dof_handler.get_triangulation())
726  .load_user_flags_line(user_flags);
727 
728  return dof_identities;
729  }
730 
731 
732 
737  template <int dim, int spacedim>
738  static std::map<types::global_dof_index, types::global_dof_index>
740  const DoFHandler<dim, spacedim> &dof_handler)
741  {
742  (void)dof_handler;
743  Assert(
744  dof_handler.hp_capability_enabled == true,
746 
747  // this function should only be called for dim<3 where there are
748  // no quad dof identies. for dim==3, the specialization below should
749  // take care of it
750  Assert(dim < 3, ExcInternalError());
751 
752  return std::map<types::global_dof_index, types::global_dof_index>();
753  }
754 
755 
756  template <int spacedim>
757  static std::map<types::global_dof_index, types::global_dof_index>
759  {
760  Assert(dof_handler.hp_capability_enabled == true,
762 
763  const int dim = 3;
764 
765  std::map<types::global_dof_index, types::global_dof_index>
766  dof_identities;
767 
768 
769  // we will mark quads that we have already treated, so first
770  // save and clear the user flags on quads and later restore
771  // them
772  std::vector<bool> user_flags;
773  dof_handler.get_triangulation().save_user_flags_quad(user_flags);
774  const_cast<::Triangulation<dim, spacedim> &>(
775  dof_handler.get_triangulation())
776  .clear_user_flags_quad();
777 
778  // An implementation of the algorithm described in the hp-
779  // paper, including the modification mentioned later in the
780  // "complications in 3-d" subsections
781  //
782  // as explained there, we do something only if there are
783  // exactly 2 finite elements associated with an object. if
784  // there is only one, then there is nothing to do anyway,
785  // and if there are 3 or more, then we can get into
786  // trouble. note that this only happens for lines in 3d and
787  // higher, and for quads only in 4d and higher, so this
788  // isn't a particularly frequent case
789  ::Table<3, std::unique_ptr<DoFIdentities>> quad_dof_identities(
790  dof_handler.fe_collection.size(),
791  dof_handler.fe_collection.size(),
792  2 /*triangle (0) or quadrilateral (1)*/);
793 
794  for (const auto &cell : dof_handler.active_cell_iterators())
795  for (const auto q : cell->face_indices())
796  if ((cell->quad(q)->user_flag_set() == false) &&
797  (cell->quad(q)->n_active_fe_indices() == 2))
798  {
799  const auto quad = cell->quad(q);
800  quad->set_user_flag();
801 
802  const std::set<unsigned int> fe_indices =
803  quad->get_active_fe_indices();
804 
805  // find out which is the most dominating finite
806  // element of the ones that are used on this quad
807  const unsigned int most_dominating_fe_index =
808  dof_handler.get_fe_collection().find_dominating_fe(
809  fe_indices,
810  /*codim=*/dim - 2);
811 
812  const unsigned int most_dominating_fe_index_face_no =
813  cell->active_fe_index() == most_dominating_fe_index ?
814  q :
815  cell->neighbor_face_no(q);
816 
817  // if we found the most dominating element, then use
818  // this to eliminate some of the degrees of freedom
819  // by identification. otherwise, the code that
820  // computes hanging node constraints will have to
821  // deal with it by computing appropriate constraints
822  // along this face/edge
823  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
824  {
825  // loop over the indices of all the finite
826  // elements that are not dominating, and
827  // identify their dofs to the most dominating
828  // one
829  for (const auto &other_fe_index : fe_indices)
830  if (other_fe_index != most_dominating_fe_index)
831  {
832  const auto &identities =
833  *ensure_existence_and_return_dof_identities<2>(
834  dof_handler.get_fe(most_dominating_fe_index),
835  dof_handler.get_fe(other_fe_index),
836  quad_dof_identities
837  [most_dominating_fe_index][other_fe_index]
838  [cell->quad(q)->reference_cell() ==
840  most_dominating_fe_index_face_no);
841 
842  for (const auto &identity : identities)
843  {
845  primary_dof_index =
846  quad->dof_index(identity.first,
847  most_dominating_fe_index);
849  dependent_dof_index =
850  quad->dof_index(identity.second,
851  other_fe_index);
852 
853  // we only store an identity if we are about to
854  // overwrite a valid degree of freedom. we will
855  // skip invalid degrees of freedom (that are
856  // associated with ghost cells) for now, and
857  // consider them later in phase 5.
858  if (dependent_dof_index !=
860  {
861  // if the DoF indices of both elements are
862  // already distributed, i.e., both of these
863  // 'fe_indices' are associated with a
864  // locally owned cell, then we should either
865  // not have a dof_identity yet, or it must
866  // come out here to be exactly as we had
867  // computed before
868  if (primary_dof_index !=
870  Assert((dof_identities.find(
871  primary_dof_index) ==
872  dof_identities.end()) ||
873  (dof_identities
874  [dependent_dof_index] ==
875  primary_dof_index),
876  ExcInternalError());
877 
878  dof_identities[dependent_dof_index] =
879  primary_dof_index;
880  }
881  }
882  }
883  }
884  }
885 
886  // finally restore the user flags
887  const_cast<::Triangulation<dim, spacedim> &>(
888  dof_handler.get_triangulation())
889  .load_user_flags_quad(user_flags);
890 
891  return dof_identities;
892  }
893 
894 
895 
900  template <int dim, int spacedim>
901  static void
904  &all_constrained_indices,
905  const DoFHandler<dim, spacedim> &dof_handler)
906  {
907  if (dof_handler.hp_capability_enabled == false)
908  return;
909 
910  Assert(all_constrained_indices.size() == dim, ExcInternalError());
911 
912  Threads::TaskGroup<> tasks;
913 
914  unsigned int i = 0;
915  tasks += Threads::new_task([&, i]() {
916  all_constrained_indices[i] =
917  compute_vertex_dof_identities(dof_handler);
918  });
919 
920  if (dim > 1)
921  {
922  ++i;
923  tasks += Threads::new_task([&, i]() {
924  all_constrained_indices[i] =
925  compute_line_dof_identities(dof_handler);
926  });
927  }
928 
929  if (dim > 2)
930  {
931  ++i;
932  tasks += Threads::new_task([&, i]() {
933  all_constrained_indices[i] =
934  compute_quad_dof_identities(dof_handler);
935  });
936  }
937 
938  tasks.join_all();
939  }
940 
941 
942 
962  template <int dim, int spacedim>
965  std::vector<types::global_dof_index> &new_dof_indices,
966  const std::vector<
967  std::map<types::global_dof_index, types::global_dof_index>>
968  &all_constrained_indices,
970  {
971  Assert(all_constrained_indices.size() == dim, ExcInternalError());
972 
973  // first preset the new DoF indices that are identities
974  for (const auto &constrained_dof_indices : all_constrained_indices)
975  for (const auto &p : constrained_dof_indices)
976  if (new_dof_indices[p.first] != numbers::invalid_dof_index)
977  {
978  Assert(new_dof_indices[p.first] == enumeration_dof_index,
979  ExcInternalError());
980 
981  new_dof_indices[p.first] = p.second;
982  }
983 
984  // then enumerate the rest
985  types::global_dof_index next_free_dof = 0;
986  for (auto &new_dof_index : new_dof_indices)
987  if (new_dof_index == enumeration_dof_index)
988  new_dof_index = next_free_dof++;
989 
990  // then loop over all those that are constrained and record the
991  // new dof number for those
992  for (const auto &constrained_dof_indices : all_constrained_indices)
993  for (const auto &p : constrained_dof_indices)
994  if (new_dof_indices[p.first] != numbers::invalid_dof_index)
995  {
996  Assert(new_dof_indices[p.first] != enumeration_dof_index,
997  ExcInternalError());
998 
999  if (p.second != numbers::invalid_dof_index)
1000  new_dof_indices[p.first] = new_dof_indices[p.second];
1001  }
1002 
1003  for (const types::global_dof_index new_dof_index : new_dof_indices)
1004  {
1005  (void)new_dof_index;
1006  Assert(new_dof_index != enumeration_dof_index,
1007  ExcInternalError());
1008  Assert(new_dof_index < next_free_dof ||
1009  new_dof_index == numbers::invalid_dof_index,
1010  ExcInternalError());
1011  }
1012 
1013  return next_free_dof;
1014  }
1015 
1016 
1017 
1026  template <int dim, int spacedim>
1029  const unsigned int n_dofs_before_identification,
1030  const bool check_validity)
1031  {
1032  if (dof_handler.hp_capability_enabled == false)
1033  return n_dofs_before_identification;
1034 
1035  std::vector<
1036  std::map<types::global_dof_index, types::global_dof_index>>
1037  all_constrained_indices(dim);
1038  compute_dof_identities(all_constrained_indices, dof_handler);
1039 
1040  std::vector<::types::global_dof_index> renumbering(
1041  n_dofs_before_identification, enumeration_dof_index);
1042  const types::global_dof_index n_dofs =
1044  all_constrained_indices,
1045  dof_handler);
1046 
1047  renumber_dofs(renumbering, IndexSet(0), dof_handler, check_validity);
1048 
1049  update_all_active_cell_dof_indices_caches(dof_handler);
1050 
1051  return n_dofs;
1052  }
1053 
1054 
1055 
1060  template <int dim, int spacedim>
1061  static void
1063  DoFHandler<dim, spacedim> &dof_handler)
1064  {
1065  Assert(
1066  dof_handler.hp_capability_enabled == true,
1068 
1069  // Note: we may wish to have something here similar to what
1070  // we do for lines and quads, namely that we only identify
1071  // dofs for any FE towards the most dominating one. however,
1072  // it is not clear whether this is actually necessary for
1073  // vertices at all, I can't think of a finite element that
1074  // would make that necessary...
1076  vertex_dof_identities(dof_handler.get_fe_collection().size(),
1077  dof_handler.get_fe_collection().size());
1078 
1079  // mark all vertices on ghost cells
1080  std::vector<bool> include_vertex(
1081  dof_handler.get_triangulation().n_vertices(), false);
1082  if (dynamic_cast<const ::parallel::
1083  DistributedTriangulationBase<dim, spacedim> *>(
1084  &dof_handler.get_triangulation()) != nullptr)
1085  for (const auto &cell : dof_handler.active_cell_iterators())
1086  if (cell->is_ghost())
1087  for (const unsigned int v : cell->vertex_indices())
1088  include_vertex[cell->vertex_index(v)] = true;
1089 
1090  // loop over all vertices and see which one we need to work on
1091  for (unsigned int vertex_index = 0;
1092  vertex_index < dof_handler.get_triangulation().n_vertices();
1093  ++vertex_index)
1094  if ((dof_handler.get_triangulation()
1095  .get_used_vertices()[vertex_index] == true) &&
1096  (include_vertex[vertex_index] == true))
1097  {
1098  const unsigned int n_active_fe_indices =
1099  ::internal::DoFAccessorImplementation::Implementation::
1100  n_active_fe_indices(dof_handler,
1101  0,
1102  vertex_index,
1103  std::integral_constant<int, 0>());
1104 
1105  if (n_active_fe_indices > 1)
1106  {
1107  const std::set<unsigned int> fe_indices =
1110  dof_handler,
1111  0,
1112  vertex_index,
1113  std::integral_constant<int, 0>());
1114 
1115  // find out which is the most dominating finite
1116  // element of the ones that are used on this vertex
1117  unsigned int most_dominating_fe_index =
1118  dof_handler.get_fe_collection().find_dominating_fe(
1119  fe_indices,
1120  /*codim=*/dim);
1121 
1122  // if we haven't found a dominating finite element,
1123  // choose the very first one to be dominant similar
1124  // to compute_vertex_dof_identities()
1125  if (most_dominating_fe_index ==
1127  most_dominating_fe_index =
1128  ::internal::DoFAccessorImplementation::
1129  Implementation::nth_active_fe_index(
1130  dof_handler,
1131  0,
1132  vertex_index,
1133  0,
1134  std::integral_constant<int, 0>());
1135 
1136  // loop over the indices of all the finite
1137  // elements that are not dominating, and
1138  // identify their dofs to the most dominating
1139  // one
1140  for (const auto &other_fe_index : fe_indices)
1141  if (other_fe_index != most_dominating_fe_index)
1142  {
1143  // make sure the entry in the equivalence
1144  // table exists
1145  const auto &identities =
1146  *ensure_existence_and_return_dof_identities<0>(
1147  dof_handler.get_fe(most_dominating_fe_index),
1148  dof_handler.get_fe(other_fe_index),
1149  vertex_dof_identities[most_dominating_fe_index]
1150  [other_fe_index]);
1151 
1152  // then loop through the identities we
1153  // have. first get the global numbers of the
1154  // dofs we want to identify and make sure they
1155  // are not yet constrained to anything else,
1156  // except for to each other. use the rule that
1157  // we will always constrain the dof with the
1158  // higher FE index to the one with the lower,
1159  // to avoid circular reasoning.
1160  for (const auto &identity : identities)
1161  {
1162  const types::global_dof_index primary_dof_index =
1163  ::internal::DoFAccessorImplementation::
1164  Implementation::get_dof_index(
1165  dof_handler,
1166  0,
1167  vertex_index,
1168  most_dominating_fe_index,
1169  identity.first,
1170  std::integral_constant<int, 0>());
1172  dependent_dof_index =
1173  ::internal::DoFAccessorImplementation::
1174  Implementation::get_dof_index(
1175  dof_handler,
1176  0,
1177  vertex_index,
1178  other_fe_index,
1179  identity.second,
1180  std::integral_constant<int, 0>());
1181 
1182  // check if we are on an interface between
1183  // a locally owned and a ghost cell on which
1184  // we need to work on.
1185  //
1186  // all degrees of freedom belonging to
1187  // dominating FE indices or to a processor
1188  // with a higher rank have been set at this
1189  // point (either in Phase 2, or after the
1190  // first ghost exchange in Phase 5). thus,
1191  // we only have to set the indices of
1192  // degrees of freedom that have been
1193  // previously flagged invalid.
1194  if ((dependent_dof_index ==
1196  (primary_dof_index !=
1198  ::internal::DoFAccessorImplementation::
1199  Implementation::set_dof_index(
1200  dof_handler,
1201  0,
1202  vertex_index,
1203  other_fe_index,
1204  identity.second,
1205  std::integral_constant<int, 0>(),
1206  primary_dof_index);
1207  }
1208  }
1209  }
1210  }
1211  }
1212 
1213 
1214 
1219  template <int spacedim>
1221  DoFHandler<1, spacedim> &dof_handler)
1222  {
1223  (void)dof_handler;
1224  Assert(dof_handler.hp_capability_enabled == true,
1226  }
1227 
1228 
1229  template <int dim, int spacedim>
1230  static void
1232  DoFHandler<dim, spacedim> &dof_handler)
1233  {
1234  Assert(
1235  dof_handler.hp_capability_enabled == true,
1237 
1238  // we will mark lines that we have already treated, so first save and
1239  // clear the user flags on lines and later restore them
1240  std::vector<bool> user_flags;
1241  dof_handler.get_triangulation().save_user_flags_line(user_flags);
1242  const_cast<::Triangulation<dim, spacedim> &>(
1243  dof_handler.get_triangulation())
1244  .clear_user_flags_line();
1245 
1246  // mark all lines on ghost cells
1247  for (const auto &cell : dof_handler.active_cell_iterators())
1248  if (cell->is_ghost())
1249  for (const auto l : cell->line_indices())
1250  cell->line(l)->set_user_flag();
1251 
1252  // An implementation of the algorithm described in the hp-paper,
1253  // including the modification mentioned later in the "complications in
1254  // 3-d" subsections
1255  //
1256  // as explained there, we do something only if there are exactly 2
1257  // finite elements associated with an object. if there is only one,
1258  // then there is nothing to do anyway, and if there are 3 or more,
1259  // then we can get into trouble. note that this only happens for lines
1260  // in 3d and higher, and for quads only in 4d and higher, so this
1261  // isn't a particularly frequent case
1262  //
1263  // there is one case, however, that we would like to handle (see, for
1264  // example, the hp/crash_15 testcase): if we have
1265  // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
1266  // then we should be able to handle this because we can simply unify
1267  // *all* dofs, not only a some. so what we do is to first treat all
1268  // pairs of finite elements that have *identical* dofs, and then only
1269  // deal with those that are not identical of which we can handle at
1270  // most 2
1271  ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
1272  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
1273 
1274  for (const auto &cell : dof_handler.active_cell_iterators())
1275  for (const auto l : cell->line_indices())
1276  if ((cell->is_locally_owned()) &&
1277  (cell->line(l)->user_flag_set() == true))
1278  {
1279  const auto line = cell->line(l);
1280  line->clear_user_flag();
1281 
1282  unsigned int unique_sets_of_dofs =
1283  line->n_active_fe_indices();
1284 
1285  // do a first loop over all sets of dofs and do identity
1286  // uniquification
1287  const unsigned int n_active_fe_indices =
1288  line->n_active_fe_indices();
1289  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1290  for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
1291  {
1292  const unsigned int fe_index_1 =
1293  line->nth_active_fe_index(f),
1294  fe_index_2 =
1295  line->nth_active_fe_index(g);
1296 
1297  if ((dof_handler.get_fe(fe_index_1).n_dofs_per_line() ==
1298  dof_handler.get_fe(fe_index_2)
1299  .n_dofs_per_line()) &&
1300  (dof_handler.get_fe(fe_index_1).n_dofs_per_line() >
1301  0))
1302  {
1303  // the number of dofs per line is identical
1304  const unsigned int dofs_per_line =
1305  dof_handler.get_fe(fe_index_1).n_dofs_per_line();
1306 
1307  const auto &identities =
1308  *ensure_existence_and_return_dof_identities<1>(
1309  dof_handler.get_fe(fe_index_1),
1310  dof_handler.get_fe(fe_index_2),
1311  line_dof_identities[fe_index_1][fe_index_2]);
1312  // see if these sets of dofs are identical. the
1313  // first condition for this is that indeed there are
1314  // n identities
1315  if (identities.size() == dofs_per_line)
1316  {
1317  unsigned int i = 0;
1318  for (; i < dofs_per_line; ++i)
1319  if ((identities[i].first != i) &&
1320  (identities[i].second != i))
1321  // not an identity
1322  break;
1323 
1324  if (i == dofs_per_line)
1325  {
1326  // The line dofs (i.e., the ones interior to
1327  // a line) of these two finite elements are
1328  // identical. Note that there could be
1329  // situations when one element still
1330  // dominates another, e.g.: FE_Q(2) x
1331  // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
1332 
1333  --unique_sets_of_dofs;
1334 
1335  // determine which one of both finite
1336  // elements is the dominating one.
1337  const std::set<unsigned int> fe_indices{
1338  fe_index_1, fe_index_2};
1339 
1340  unsigned int dominating_fe_index =
1341  dof_handler.get_fe_collection()
1342  .find_dominating_fe(fe_indices,
1343  /*codim*/ dim - 1);
1344  unsigned int other_fe_index =
1346 
1347  if (dominating_fe_index !=
1349  other_fe_index =
1350  (dominating_fe_index == fe_index_1) ?
1351  fe_index_2 :
1352  fe_index_1;
1353  else
1354  {
1355  // if we haven't found a dominating
1356  // finite element, choose the one with
1357  // the lower index to be dominating
1358  dominating_fe_index = fe_index_1;
1359  other_fe_index = fe_index_2;
1360  }
1361 
1362  for (unsigned int j = 0; j < dofs_per_line;
1363  ++j)
1364  {
1366  primary_dof_index = line->dof_index(
1367  j, dominating_fe_index);
1369  dependent_dof_index =
1370  line->dof_index(j, other_fe_index);
1371 
1372  // check if we are on an interface
1373  // between a locally owned and a ghost
1374  // cell on which we need to work on.
1375  //
1376  // all degrees of freedom belonging to
1377  // dominating fe_indices or to a
1378  // processor with a higher rank have
1379  // been set at this point (either in
1380  // Phase 2, or after the first ghost
1381  // exchange in Phase 5). thus, we only
1382  // have to set the indices of degrees
1383  // of freedom that have been previously
1384  // flagged invalid.
1385  if ((dependent_dof_index ==
1387  (primary_dof_index !=
1389  line->set_dof_index(j,
1390  primary_dof_index,
1391  fe_index_2);
1392  }
1393  }
1394  }
1395  }
1396  }
1397 
1398  // if at this point, there is only one unique set of dofs
1399  // left, then we have taken care of everything above. if there
1400  // are two, then we need to deal with them here. if there are
1401  // more, then we punt, as described in the paper (and
1402  // mentioned above)
1403  // TODO: The check for 'dim==2' was inserted by intuition. It
1404  // fixes
1405  // the previous problems with @ref step_27 "step-27" in 3D. But an
1406  // explanation for this is still required, and what we do here
1407  // is not what we describe in the paper!.
1408  if ((unique_sets_of_dofs == 2) && (dim == 2))
1409  {
1410  const std::set<unsigned int> fe_indices =
1411  line->get_active_fe_indices();
1412 
1413  // find out which is the most dominating finite element of
1414  // the ones that are used on this line
1415  const unsigned int most_dominating_fe_index =
1416  dof_handler.get_fe_collection().find_dominating_fe(
1417  fe_indices,
1418  /*codim=*/dim - 1);
1419 
1420  // if we found the most dominating element, then use this
1421  // to eliminate some of the degrees of freedom by
1422  // identification. otherwise, the code that computes
1423  // hanging node constraints will have to deal with it by
1424  // computing appropriate constraints along this face/edge
1425  if (most_dominating_fe_index !=
1427  {
1428  // loop over the indices of all the finite elements
1429  // that are not dominating, and identify their dofs to
1430  // the most dominating one
1431  for (const auto &other_fe_index : fe_indices)
1432  if (other_fe_index != most_dominating_fe_index)
1433  {
1434  const auto &identities =
1435  *ensure_existence_and_return_dof_identities<
1436  1>(dof_handler.get_fe(
1437  most_dominating_fe_index),
1438  dof_handler.get_fe(other_fe_index),
1439  line_dof_identities
1440  [most_dominating_fe_index]
1441  [other_fe_index]);
1442 
1443  for (const auto &identity : identities)
1444  {
1446  primary_dof_index = line->dof_index(
1447  identity.first,
1448  most_dominating_fe_index);
1450  dependent_dof_index =
1451  line->dof_index(identity.second,
1452  other_fe_index);
1453 
1454  // check if we are on an interface between
1455  // a locally owned and a ghost cell on which
1456  // we need to work on.
1457  //
1458  // all degrees of freedom belonging to
1459  // dominating FE indices or to a processor
1460  // with a higher rank have been set at this
1461  // point (either in Phase 2, or after the
1462  // first ghost exchange in Phase 5). thus,
1463  // we only have to set the indices of
1464  // degrees of freedom that have been
1465  // previously flagged invalid.
1466  if ((dependent_dof_index ==
1468  (primary_dof_index !=
1470  line->set_dof_index(identity.second,
1471  primary_dof_index,
1472  other_fe_index);
1473  }
1474  }
1475  }
1476  }
1477  }
1478 
1479  // finally restore the user flags
1480  const_cast<::Triangulation<dim, spacedim> &>(
1481  dof_handler.get_triangulation())
1482  .load_user_flags_line(user_flags);
1483  }
1484 
1485 
1486 
1491  template <int dim, int spacedim>
1492  static void
1494  DoFHandler<dim, spacedim> &dof_handler)
1495  {
1496  (void)dof_handler;
1497  Assert(
1498  dof_handler.hp_capability_enabled == true,
1500 
1501  // this function should only be called for dim<3 where there are
1502  // no quad dof identies. for dim>=3, the specialization below should
1503  // take care of it
1504  Assert(dim < 3, ExcInternalError());
1505  }
1506 
1507 
1508  template <int spacedim>
1510  DoFHandler<3, spacedim> &dof_handler)
1511  {
1512  Assert(dof_handler.hp_capability_enabled == true,
1514 
1515  const int dim = 3;
1516 
1517  // we will mark quads that we have already treated, so first
1518  // save and clear the user flags on quads and later restore
1519  // them
1520  std::vector<bool> user_flags;
1521  dof_handler.get_triangulation().save_user_flags_quad(user_flags);
1522  const_cast<::Triangulation<dim, spacedim> &>(
1523  dof_handler.get_triangulation())
1524  .clear_user_flags_quad();
1525 
1526  // mark all quads on ghost cells
1527  for (const auto &cell : dof_handler.active_cell_iterators())
1528  if (cell->is_ghost())
1529  for (const auto q : cell->face_indices())
1530  cell->quad(q)->set_user_flag();
1531 
1532  // An implementation of the algorithm described in the hp-
1533  // paper, including the modification mentioned later in the
1534  // "complications in 3-d" subsections
1535  //
1536  // as explained there, we do something only if there are
1537  // exactly 2 finite elements associated with an object. if
1538  // there is only one, then there is nothing to do anyway,
1539  // and if there are 3 or more, then we can get into
1540  // trouble. note that this only happens for lines in 3d and
1541  // higher, and for quads only in 4d and higher, so this
1542  // isn't a particularly frequent case
1543  ::Table<3, std::unique_ptr<DoFIdentities>> quad_dof_identities(
1544  dof_handler.fe_collection.size(),
1545  dof_handler.fe_collection.size(),
1546  2 /*triangle (0) or quadrilateral (1)*/);
1547 
1548  for (const auto &cell : dof_handler.active_cell_iterators())
1549  for (const auto q : cell->face_indices())
1550  if ((cell->is_locally_owned()) &&
1551  (cell->quad(q)->user_flag_set() == true) &&
1552  (cell->quad(q)->n_active_fe_indices() == 2))
1553  {
1554  const auto quad = cell->quad(q);
1555  quad->clear_user_flag();
1556 
1557  const std::set<unsigned int> fe_indices =
1558  quad->get_active_fe_indices();
1559 
1560  // find out which is the most dominating finite
1561  // element of the ones that are used on this quad
1562  const unsigned int most_dominating_fe_index =
1563  dof_handler.get_fe_collection().find_dominating_fe(
1564  fe_indices,
1565  /*codim=*/dim - 2);
1566 
1567  const unsigned int most_dominating_fe_index_face_no =
1568  cell->active_fe_index() == most_dominating_fe_index ?
1569  q :
1570  cell->neighbor_face_no(q);
1571 
1572  // if we found the most dominating element, then use
1573  // this to eliminate some of the degrees of freedom
1574  // by identification. otherwise, the code that
1575  // computes hanging node constraints will have to
1576  // deal with it by computing appropriate constraints
1577  // along this face/edge
1578  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
1579  {
1580  // loop over the indices of all the finite
1581  // elements that are not dominating, and
1582  // identify their dofs to the most dominating
1583  // one
1584  for (const auto &other_fe_index : fe_indices)
1585  if (other_fe_index != most_dominating_fe_index)
1586  {
1587  const auto &identities =
1588  *ensure_existence_and_return_dof_identities<2>(
1589  dof_handler.get_fe(most_dominating_fe_index),
1590  dof_handler.get_fe(other_fe_index),
1591  quad_dof_identities
1592  [most_dominating_fe_index][other_fe_index]
1593  [cell->quad(q)->reference_cell() ==
1595  most_dominating_fe_index_face_no);
1596 
1597  for (const auto &identity : identities)
1598  {
1600  primary_dof_index =
1601  quad->dof_index(identity.first,
1602  most_dominating_fe_index);
1604  dependent_dof_index =
1605  quad->dof_index(identity.second,
1606  other_fe_index);
1607 
1608  // check if we are on an interface between
1609  // a locally owned and a ghost cell on which
1610  // we need to work on.
1611  //
1612  // all degrees of freedom belonging to
1613  // dominating FE indices or to a processor with
1614  // a higher rank have been set at this point
1615  // (either in Phase 2, or after the first ghost
1616  // exchange in Phase 5). thus, we only have to
1617  // set the indices of degrees of freedom that
1618  // have been previously flagged invalid.
1619  if ((dependent_dof_index ==
1621  (primary_dof_index !=
1623  quad->set_dof_index(identity.second,
1624  primary_dof_index,
1625  other_fe_index);
1626  }
1627  }
1628  }
1629  }
1630 
1631  // finally restore the user flags
1632  const_cast<::Triangulation<dim, spacedim> &>(
1633  dof_handler.get_triangulation())
1634  .load_user_flags_quad(user_flags);
1635  }
1636 
1637 
1638 
1651  template <int dim, int spacedim>
1652  static void
1654  DoFHandler<dim, spacedim> &dof_handler)
1655  {
1656  if (dof_handler.hp_capability_enabled == false)
1657  return;
1658 
1659  {
1660  Threads::TaskGroup<> tasks;
1661 
1662  tasks += Threads::new_task([&]() {
1664  });
1665 
1666  if (dim > 1)
1667  {
1668  tasks += Threads::new_task([&]() {
1670  });
1671  }
1672 
1673  if (dim > 2)
1674  {
1675  tasks += Threads::new_task([&]() {
1677  });
1678  }
1679 
1680  tasks.join_all();
1681  }
1682 
1683  update_all_active_cell_dof_indices_caches(dof_handler);
1684  }
1685 
1686 
1687 
1694  template <int dim, int spacedim>
1697  DoFHandler<dim, spacedim> &dof_handler)
1698  {
1699  Assert(dof_handler.get_triangulation().n_levels() > 0,
1700  ExcMessage("Empty triangulation"));
1701 
1702  // Step 1: distribute dofs on all cells, but definitely
1703  // exclude artificial cells
1704  types::global_dof_index next_free_dof = 0;
1705 
1706  std::vector<types::global_dof_index> dof_indices;
1707 
1708  for (auto cell : dof_handler.active_cell_iterators())
1709  if (!cell->is_artificial())
1710  if ((subdomain_id == numbers::invalid_subdomain_id) ||
1711  (cell->subdomain_id() == subdomain_id))
1712  {
1713  dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1714 
1715  // circumvent cache
1716  internal::DoFAccessorImplementation::Implementation::
1717  get_dof_indices(*cell,
1718  dof_indices,
1719  cell->active_fe_index());
1720 
1721  for (auto &dof_index : dof_indices)
1722  if (dof_index == numbers::invalid_dof_index)
1723  dof_index = next_free_dof++;
1724 
1725  cell->set_dof_indices(dof_indices);
1726  }
1727 
1728  update_all_active_cell_dof_indices_caches(dof_handler);
1729 
1730  return next_free_dof;
1731  }
1732 
1733 
1734 
1748  template <int dim, int spacedim>
1749  static void
1751  std::vector<types::global_dof_index> &renumbering,
1753  const DoFHandler<dim, spacedim> & dof_handler)
1754  {
1755  std::vector<types::global_dof_index> local_dof_indices;
1756 
1757  for (const auto &cell : dof_handler.active_cell_iterators())
1758  if (cell->is_ghost() && (cell->subdomain_id() < subdomain_id))
1759  {
1760  // we found a neighboring ghost cell whose subdomain
1761  // is "stronger" than our own subdomain
1762 
1763  // delete all dofs that live there and that we have
1764  // previously assigned a number to (i.e. the ones on
1765  // the interface)
1766  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1767  cell->get_dof_indices(local_dof_indices);
1768  for (const auto &local_dof_index : local_dof_indices)
1769  if (local_dof_index != numbers::invalid_dof_index)
1770  renumbering[local_dof_index] = numbers::invalid_dof_index;
1771  }
1772  }
1773 
1774 
1775 
1776  /* -------------- distribute_mg_dofs functionality ------------- */
1777 
1778 
1779 
1780  template <int dim, int spacedim>
1783  DoFHandler<dim, spacedim> &dof_handler,
1784  const unsigned int level)
1785  {
1786  Assert(dof_handler.hp_capability_enabled == false,
1787  ExcInternalError());
1788 
1789  const ::Triangulation<dim, spacedim> &tria =
1790  dof_handler.get_triangulation();
1791  Assert(tria.n_levels() > 0, ExcMessage("Empty triangulation"));
1792  if (level >= tria.n_levels())
1793  return 0; // this is allowed for multigrid
1794 
1795  types::global_dof_index next_free_dof = 0;
1796 
1797  std::vector<types::global_dof_index> dof_indices;
1798 
1799  for (auto cell : dof_handler.cell_iterators_on_level(level))
1800  if ((level_subdomain_id == numbers::invalid_subdomain_id) ||
1801  (cell->level_subdomain_id() == level_subdomain_id))
1802  {
1803  dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1804 
1805  cell->get_mg_dof_indices(dof_indices);
1806 
1807  for (auto &dof_index : dof_indices)
1808  if (dof_index == numbers::invalid_dof_index)
1809  dof_index = next_free_dof++;
1810 
1811  cell->set_mg_dof_indices(dof_indices);
1812  }
1813 
1814  return next_free_dof;
1815  }
1816 
1817 
1818 
1819  /* --------------------- renumber_dofs functionality ---------------- */
1820 
1821 
1829  template <int dim, int spacedim>
1830  static void
1832  const std::vector<types::global_dof_index> &new_numbers,
1833  const IndexSet & indices_we_care_about,
1834  DoFHandler<dim, spacedim> & dof_handler)
1835  {
1836  for (unsigned int d = 1; d < dim; d++)
1837  for (auto &i : dof_handler.object_dof_indices[0][d])
1838  if (i != numbers::invalid_dof_index)
1839  i = ((indices_we_care_about.size() == 0) ?
1840  new_numbers[i] :
1841  new_numbers[indices_we_care_about.index_within_set(i)]);
1842  }
1843 
1844 
1845 
1846  template <int dim, int spacedim>
1847  static void
1849  const std::vector<types::global_dof_index> &new_numbers,
1850  const IndexSet & indices_we_care_about,
1851  DoFHandler<dim, spacedim> & dof_handler,
1852  const bool check_validity)
1853  {
1854  if (dof_handler.hp_capability_enabled == false)
1855  {
1856  // we can not use cell iterators in this function since then
1857  // we would renumber the dofs on the interface of two cells
1858  // more than once. Anyway, this way it's not only more
1859  // correct but also faster; note, however, that dof numbers
1860  // may be invalid_dof_index, namely when the appropriate
1861  // vertex/line/etc is unused
1862  for (std::vector<types::global_dof_index>::iterator i =
1863  dof_handler.object_dof_indices[0][0].begin();
1864  i != dof_handler.object_dof_indices[0][0].end();
1865  ++i)
1866  if (*i != numbers::invalid_dof_index)
1867  *i =
1868  (indices_we_care_about.size() == 0) ?
1869  (new_numbers[*i]) :
1870  (new_numbers[indices_we_care_about.index_within_set(*i)]);
1871  else if (check_validity)
1872  // if index is invalid_dof_index: check if this one
1873  // really is unused
1874  Assert(dof_handler.get_triangulation().vertex_used(
1875  (i - dof_handler.object_dof_indices[0][0].begin()) /
1876  dof_handler.get_fe().n_dofs_per_vertex()) == false,
1877  ExcInternalError());
1878  return;
1879  }
1880 
1881 
1882  for (unsigned int vertex_index = 0;
1883  vertex_index < dof_handler.get_triangulation().n_vertices();
1884  ++vertex_index)
1885  {
1886  const unsigned int n_active_fe_indices =
1887  ::internal::DoFAccessorImplementation::Implementation::
1888  n_active_fe_indices(dof_handler,
1889  0,
1890  vertex_index,
1891  std::integral_constant<int, 0>());
1892 
1893  // if this vertex is unused, then we really ought not to have
1894  // allocated any space for it, i.e., n_active_fe_indices should be
1895  // zero, and there is no space to actually store dof indices for
1896  // this vertex
1897  if (dof_handler.get_triangulation().vertex_used(vertex_index) ==
1898  false)
1899  Assert(n_active_fe_indices == 0, ExcInternalError());
1900 
1901  // otherwise the vertex is used; it may still not hold any dof
1902  // indices if it is located on an artificial cell and not adjacent
1903  // to a ghost cell, but in that case there is simply nothing for
1904  // us to do
1905  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1906  {
1907  const unsigned int fe_index =
1908  ::internal::DoFAccessorImplementation::
1909  Implementation::nth_active_fe_index(
1910  dof_handler,
1911  0,
1912  vertex_index,
1913  f,
1914  std::integral_constant<int, 0>());
1915 
1916  for (unsigned int d = 0;
1917  d < dof_handler.get_fe(fe_index).n_dofs_per_vertex();
1918  ++d)
1919  {
1920  const types::global_dof_index old_dof_index =
1921  ::internal::DoFAccessorImplementation::
1922  Implementation::get_dof_index(
1923  dof_handler,
1924  0,
1925  vertex_index,
1926  fe_index,
1927  d,
1928  std::integral_constant<int, 0>());
1929 
1930  // if check_validity was set, then we are to verify that
1931  // the previous indices were all valid. this really should
1932  // be the case: we allocated space for these vertex dofs,
1933  // i.e., at least one adjacent cell has a valid
1934  // active FE index, so there are DoFs that really live
1935  // on this vertex. if check_validity is set, then we
1936  // must make sure that they have been set to something
1937  // useful
1938  if (check_validity)
1939  Assert(old_dof_index != numbers::invalid_dof_index,
1940  ExcInternalError());
1941 
1942  if (old_dof_index != numbers::invalid_dof_index)
1943  {
1944  // In the following blocks, we first check whether
1945  // we were given an IndexSet of DoFs to touch. If not
1946  // (the first 'if' case here), then we are in the
1947  // sequential case and are allowed to touch all DoFs.
1948  //
1949  // If yes (the 'else' case), then we need to
1950  // distinguish whether the DoF whose number we want to
1951  // touch is in fact locally owned (i.e., is in the
1952  // index set) and then we can actually assign it a new
1953  // number; otherwise, we have encountered a
1954  // non-locally owned DoF for which we don't know the
1955  // new number yet and so set it to an invalid index.
1956  // This will later be fixed up after the first ghost
1957  // exchange phase when we unify hp-DoFs on neighboring
1958  // cells.
1959  if (indices_we_care_about.size() == 0)
1962  dof_handler,
1963  0,
1964  vertex_index,
1965  fe_index,
1966  d,
1967  std::integral_constant<int, 0>(),
1968  new_numbers[old_dof_index]);
1969  else
1970  {
1971  if (indices_we_care_about.is_element(
1972  old_dof_index))
1975  dof_handler,
1976  0,
1977  vertex_index,
1978  fe_index,
1979  d,
1980  std::integral_constant<int, 0>(),
1981  new_numbers[indices_we_care_about
1983  old_dof_index)]);
1984  else
1985  ::internal::DoFAccessorImplementation::
1986  Implementation::set_dof_index(
1987  dof_handler,
1988  0,
1989  vertex_index,
1990  fe_index,
1991  d,
1992  std::integral_constant<int, 0>(),
1994  }
1995  }
1996  }
1997  }
1998  }
1999  }
2000 
2001 
2002 
2003  template <int dim, int spacedim>
2004  static void
2006  const std::vector<types::global_dof_index> &new_numbers,
2007  const IndexSet & indices_we_care_about,
2008  DoFHandler<dim, spacedim> & dof_handler)
2009  {
2010  if (dof_handler.hp_capability_enabled == false)
2011  {
2012  for (unsigned int level = 0;
2013  level < dof_handler.object_dof_indices.size();
2014  ++level)
2015  for (auto &i : dof_handler.object_dof_indices[level][dim])
2016  if (i != numbers::invalid_dof_index)
2017  i = ((indices_we_care_about.size() == 0) ?
2018  new_numbers[i] :
2019  new_numbers[indices_we_care_about.index_within_set(
2020  i)]);
2021  return;
2022  }
2023 
2024  for (const auto &cell : dof_handler.active_cell_iterators())
2025  if (!cell->is_artificial())
2026  {
2027  const unsigned int fe_index = cell->active_fe_index();
2028 
2029  for (unsigned int d = 0;
2030  d < dof_handler.get_fe(fe_index)
2031  .template n_dofs_per_object<dim>();
2032  ++d)
2033  {
2034  const types::global_dof_index old_dof_index =
2035  cell->dof_index(d, fe_index);
2036  if (old_dof_index != numbers::invalid_dof_index)
2037  {
2038  // In the following blocks, we first check whether
2039  // we were given an IndexSet of DoFs to touch. If not
2040  // (the first 'if' case here), then we are in the
2041  // sequential case and are allowed to touch all DoFs.
2042  //
2043  // If yes (the 'else' case), then we need to distinguish
2044  // whether the DoF whose number we want to touch is in
2045  // fact locally owned (i.e., is in the index set) and
2046  // then we can actually assign it a new number;
2047  // otherwise, we have encountered a non-locally owned
2048  // DoF for which we don't know the new number yet and so
2049  // set it to an invalid index. This will later be fixed
2050  // up after the first ghost exchange phase when we unify
2051  // hp-DoFs on neighboring cells.
2052  if (indices_we_care_about.size() == 0)
2053  cell->set_dof_index(d,
2054  new_numbers[old_dof_index],
2055  fe_index);
2056  else
2057  {
2058  if (indices_we_care_about.is_element(old_dof_index))
2059  cell->set_dof_index(
2060  d,
2061  new_numbers[indices_we_care_about
2062  .index_within_set(old_dof_index)],
2063  fe_index);
2064  else
2065  cell->set_dof_index(d,
2067  fe_index);
2068  }
2069  }
2070  }
2071  }
2072  }
2073 
2074 
2075 
2076  template <int spacedim>
2077  static void
2079  const std::vector<types::global_dof_index> & /*new_numbers*/,
2080  const IndexSet & /*indices_we_care_about*/,
2081  DoFHandler<1, spacedim> & /*dof_handler*/)
2082  {
2083  // nothing to do in 1d since there are no separate faces -- we've
2084  // already taken care of this when dealing with the vertices
2085  }
2086 
2087 
2088 
2089  template <int spacedim>
2090  static void
2092  const std::vector<types::global_dof_index> &new_numbers,
2093  const IndexSet & indices_we_care_about,
2094  DoFHandler<2, spacedim> & dof_handler)
2095  {
2096  const unsigned int dim = 2;
2097 
2098  if (dof_handler.hp_capability_enabled == false)
2099  {
2100  for (unsigned int d = 1; d < dim; d++)
2101  for (auto &i : dof_handler.object_dof_indices[0][d])
2102  if (i != numbers::invalid_dof_index)
2103  i = ((indices_we_care_about.size() == 0) ?
2104  new_numbers[i] :
2105  new_numbers[indices_we_care_about.index_within_set(
2106  i)]);
2107  return;
2108  }
2109 
2110  // deal with DoFs on lines
2111  {
2112  // save user flags on lines so we can use them to mark lines
2113  // we've already treated
2114  std::vector<bool> saved_line_user_flags;
2115  const_cast<::Triangulation<dim, spacedim> &>(
2116  dof_handler.get_triangulation())
2117  .save_user_flags_line(saved_line_user_flags);
2118  const_cast<::Triangulation<dim, spacedim> &>(
2119  dof_handler.get_triangulation())
2120  .clear_user_flags_line();
2121 
2122  for (const auto &cell : dof_handler.active_cell_iterators())
2123  if (!cell->is_artificial())
2124  for (const auto l : cell->line_indices())
2125  if (cell->line(l)->user_flag_set() == false)
2126  {
2127  const auto line = cell->line(l);
2128  line->set_user_flag();
2129 
2130  const unsigned int n_active_fe_indices =
2131  line->n_active_fe_indices();
2132 
2133  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2134  {
2135  const unsigned int fe_index =
2136  line->nth_active_fe_index(f);
2137 
2138  for (unsigned int d = 0;
2139  d <
2140  dof_handler.get_fe(fe_index).n_dofs_per_line();
2141  ++d)
2142  {
2143  const types::global_dof_index old_dof_index =
2144  line->dof_index(d, fe_index);
2145  if (old_dof_index != numbers::invalid_dof_index)
2146  {
2147  // In the following blocks, we first check
2148  // whether we were given an IndexSet of DoFs
2149  // to touch. If not (the first 'if' case
2150  // here), then we are in the sequential case
2151  // and are allowed to touch all DoFs.
2152  //
2153  // If yes (the 'else' case), then we need to
2154  // distinguish whether the DoF whose number we
2155  // want to touch is in fact locally owned
2156  // (i.e., is in the index set) and then we can
2157  // actually assign it a new number; otherwise,
2158  // we have encountered a non-locally owned DoF
2159  // for which we don't know the new number yet
2160  // and so set it to an invalid index. This
2161  // will later be fixed up after the first
2162  // ghost exchange phase when we unify hp-DoFs
2163  // on neighboring cells.
2164  if (indices_we_care_about.size() == 0)
2165  line->set_dof_index(
2166  d, new_numbers[old_dof_index], fe_index);
2167  else
2168  {
2169  if (indices_we_care_about.is_element(
2170  old_dof_index))
2171  line->set_dof_index(
2172  d,
2173  new_numbers[indices_we_care_about
2175  old_dof_index)],
2176  fe_index);
2177  else
2178  line->set_dof_index(
2179  d,
2181  fe_index);
2182  }
2183  }
2184  }
2185  }
2186  }
2187 
2188  // at the end, restore the user
2189  // flags for the lines
2190  const_cast<::Triangulation<dim, spacedim> &>(
2191  dof_handler.get_triangulation())
2192  .load_user_flags_line(saved_line_user_flags);
2193  }
2194  }
2195 
2196 
2197 
2198  template <int spacedim>
2199  static void
2201  const std::vector<types::global_dof_index> &new_numbers,
2202  const IndexSet & indices_we_care_about,
2203  DoFHandler<3, spacedim> & dof_handler)
2204  {
2205  const unsigned int dim = 3;
2206 
2207  if (dof_handler.hp_capability_enabled == false)
2208  {
2209  for (unsigned int d = 1; d < dim; d++)
2210  for (auto &i : dof_handler.object_dof_indices[0][d])
2211  if (i != numbers::invalid_dof_index)
2212  i = ((indices_we_care_about.size() == 0) ?
2213  new_numbers[i] :
2214  new_numbers[indices_we_care_about.index_within_set(
2215  i)]);
2216  return;
2217  }
2218 
2219  // deal with DoFs on lines
2220  {
2221  // save user flags on lines so we can use them to mark lines
2222  // we've already treated
2223  std::vector<bool> saved_line_user_flags;
2224  const_cast<::Triangulation<dim, spacedim> &>(
2225  dof_handler.get_triangulation())
2226  .save_user_flags_line(saved_line_user_flags);
2227  const_cast<::Triangulation<dim, spacedim> &>(
2228  dof_handler.get_triangulation())
2229  .clear_user_flags_line();
2230 
2231  for (const auto &cell : dof_handler.active_cell_iterators())
2232  if (!cell->is_artificial())
2233  for (const auto l : cell->line_indices())
2234  if (cell->line(l)->user_flag_set() == false)
2235  {
2236  const auto line = cell->line(l);
2237  line->set_user_flag();
2238 
2239  const unsigned int n_active_fe_indices =
2240  line->n_active_fe_indices();
2241 
2242  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2243  {
2244  const unsigned int fe_index =
2245  line->nth_active_fe_index(f);
2246 
2247  for (unsigned int d = 0;
2248  d <
2249  dof_handler.get_fe(fe_index).n_dofs_per_line();
2250  ++d)
2251  {
2252  const types::global_dof_index old_dof_index =
2253  line->dof_index(d, fe_index);
2254  if (old_dof_index != numbers::invalid_dof_index)
2255  {
2256  // In the following blocks, we first check
2257  // whether we were given an IndexSet of DoFs
2258  // to touch. If not (the first 'if' case
2259  // here), then we are in the sequential case
2260  // and are allowed to touch all DoFs.
2261  //
2262  // If yes (the 'else' case), then we need to
2263  // distinguish whether the DoF whose number we
2264  // want to touch is in fact locally owned
2265  // (i.e., is in the index set) and then we can
2266  // actually assign it a new number; otherwise,
2267  // we have encountered a non-locally owned DoF
2268  // for which we don't know the new number yet
2269  // and so set it to an invalid index. This
2270  // will later be fixed up after the first
2271  // ghost exchange phase when we unify hp-DoFs
2272  // on neighboring cells.
2273  if (indices_we_care_about.size() == 0)
2274  line->set_dof_index(
2275  d, new_numbers[old_dof_index], fe_index);
2276  else if (indices_we_care_about.is_element(
2277  old_dof_index))
2278  line->set_dof_index(
2279  d,
2280  new_numbers[indices_we_care_about
2282  old_dof_index)],
2283  fe_index);
2284  else
2285  line->set_dof_index(
2286  d, numbers::invalid_dof_index, fe_index);
2287  }
2288  }
2289  }
2290  }
2291 
2292  // at the end, restore the user
2293  // flags for the lines
2294  const_cast<::Triangulation<dim, spacedim> &>(
2295  dof_handler.get_triangulation())
2296  .load_user_flags_line(saved_line_user_flags);
2297  }
2298 
2299  // then deal with dofs on quads
2300  {
2301  std::vector<bool> saved_quad_user_flags;
2302  const_cast<::Triangulation<dim, spacedim> &>(
2303  dof_handler.get_triangulation())
2304  .save_user_flags_quad(saved_quad_user_flags);
2305  const_cast<::Triangulation<dim, spacedim> &>(
2306  dof_handler.get_triangulation())
2307  .clear_user_flags_quad();
2308 
2309  for (const auto &cell : dof_handler.active_cell_iterators())
2310  if (!cell->is_artificial())
2311  for (const auto q : cell->face_indices())
2312  if (cell->quad(q)->user_flag_set() == false)
2313  {
2314  const auto quad = cell->quad(q);
2315  quad->set_user_flag();
2316 
2317  const unsigned int n_active_fe_indices =
2318  quad->n_active_fe_indices();
2319 
2320  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2321  {
2322  const unsigned int fe_index =
2323  quad->nth_active_fe_index(f);
2324 
2325  for (unsigned int d = 0;
2326  d <
2327  dof_handler.get_fe(fe_index).n_dofs_per_quad(q);
2328  ++d)
2329  {
2330  const types::global_dof_index old_dof_index =
2331  quad->dof_index(d, fe_index);
2332  if (old_dof_index != numbers::invalid_dof_index)
2333  {
2334  // In the following blocks, we first check
2335  // whether we were given an IndexSet of DoFs
2336  // to touch. If not (the first 'if' case
2337  // here), then we are in the sequential case
2338  // and are allowed to touch all DoFs.
2339  //
2340  // If yes (the 'else' case), then we need to
2341  // distinguish whether the DoF whose number we
2342  // want to touch is in fact locally owned
2343  // (i.e., is in the index set) and then we can
2344  // actually assign it a new number; otherwise,
2345  // we have encountered a non-locally owned DoF
2346  // for which we don't know the new number yet
2347  // and so set it to an invalid index. This
2348  // will later be fixed up after the first
2349  // ghost exchange phase when we unify hp-DoFs
2350  // on neighboring cells.
2351  if (indices_we_care_about.size() == 0)
2352  quad->set_dof_index(
2353  d, new_numbers[old_dof_index], fe_index);
2354  else
2355  {
2356  if (indices_we_care_about.is_element(
2357  old_dof_index))
2358  quad->set_dof_index(
2359  d,
2360  new_numbers[indices_we_care_about
2362  old_dof_index)],
2363  fe_index);
2364  else
2365  quad->set_dof_index(
2366  d,
2368  fe_index);
2369  }
2370  }
2371  }
2372  }
2373  }
2374 
2375  // at the end, restore the user flags for the quads
2376  const_cast<::Triangulation<dim, spacedim> &>(
2377  dof_handler.get_triangulation())
2378  .load_user_flags_quad(saved_quad_user_flags);
2379  }
2380  }
2381 
2382 
2383 
2395  template <int dim, int space_dim>
2396  static void
2397  renumber_dofs(const std::vector<types::global_dof_index> &new_numbers,
2398  const IndexSet & indices_we_care_about,
2399  const DoFHandler<dim, space_dim> &dof_handler,
2400  const bool check_validity)
2401  {
2402  if (dim == 1)
2403  Assert(indices_we_care_about == IndexSet(0), ExcNotImplemented());
2404 
2405  // renumber DoF indices on vertices, cells, and faces. this
2406  // can be done in parallel because the respective functions
2407  // work on separate data structures
2408  Threads::TaskGroup<> tasks;
2409  tasks += Threads::new_task([&]() {
2410  renumber_vertex_dofs(new_numbers,
2411  indices_we_care_about,
2412  const_cast<DoFHandler<dim, space_dim> &>(
2413  dof_handler),
2414  check_validity);
2415  });
2416  tasks += Threads::new_task([&]() {
2417  renumber_face_dofs(new_numbers,
2418  indices_we_care_about,
2419  const_cast<DoFHandler<dim, space_dim> &>(
2420  dof_handler));
2421  });
2422  tasks += Threads::new_task([&]() {
2423  renumber_cell_dofs(new_numbers,
2424  indices_we_care_about,
2425  const_cast<DoFHandler<dim, space_dim> &>(
2426  dof_handler));
2427  });
2428  tasks.join_all();
2429 
2430  // update the cache used for cell dof indices
2431  update_all_active_cell_dof_indices_caches(
2432  const_cast<DoFHandler<dim, space_dim> &>(dof_handler));
2433  }
2434 
2435 
2436 
2437  /* --------------------- renumber_mg_dofs functionality ----------------
2438  */
2439 
2447  template <int dim, int spacedim>
2448  static void
2450  const std::vector<::types::global_dof_index> &new_numbers,
2451  const IndexSet & indices_we_care_about,
2452  DoFHandler<dim, spacedim> &dof_handler,
2453  const unsigned int level,
2454  const bool check_validity)
2455  {
2456  (void)check_validity;
2457  Assert(level < dof_handler.get_triangulation().n_levels(),
2458  ExcInternalError());
2459 
2460  for (typename std::vector<
2461  typename DoFHandler<dim, spacedim>::MGVertexDoFs>::iterator i =
2462  dof_handler.mg_vertex_dofs.begin();
2463  i != dof_handler.mg_vertex_dofs.end();
2464  ++i)
2465  // if the present vertex lives on the current level
2466  if ((i->get_coarsest_level() <= level) &&
2467  (i->get_finest_level() >= level))
2468  for (unsigned int d = 0;
2469  d < dof_handler.get_fe().n_dofs_per_vertex();
2470  ++d)
2471  {
2473  i->get_index(level,
2474  d,
2475  dof_handler.get_fe().n_dofs_per_vertex());
2476 
2477  if (idx != numbers::invalid_dof_index)
2478  {
2479  Assert(check_validity == false ||
2480  (indices_we_care_about.size() > 0 ?
2481  indices_we_care_about.is_element(idx) :
2482  (idx < new_numbers.size())),
2483  ExcInternalError());
2484  i->set_index(level,
2485  d,
2486  dof_handler.get_fe().n_dofs_per_vertex(),
2487  (indices_we_care_about.size() == 0) ?
2488  (new_numbers[idx]) :
2489  (new_numbers[indices_we_care_about
2490  .index_within_set(idx)]));
2491  }
2492  }
2493  }
2494 
2495 
2496 
2504  template <int dim, int spacedim>
2505  static void
2507  const std::vector<::types::global_dof_index> &new_numbers,
2508  const IndexSet & indices_we_care_about,
2509  DoFHandler<dim, spacedim> &dof_handler,
2510  const unsigned int level)
2511  {
2512  for (std::vector<types::global_dof_index>::iterator i =
2513  dof_handler.mg_levels[level]->dof_object.dofs.begin();
2514  i != dof_handler.mg_levels[level]->dof_object.dofs.end();
2515  ++i)
2516  {
2517  if (*i != numbers::invalid_dof_index)
2518  {
2519  Assert((indices_we_care_about.size() > 0 ?
2520  indices_we_care_about.is_element(*i) :
2521  (*i < new_numbers.size())),
2522  ExcInternalError());
2523  *i =
2524  (indices_we_care_about.size() == 0) ?
2525  (new_numbers[*i]) :
2526  (new_numbers[indices_we_care_about.index_within_set(*i)]);
2527  }
2528  }
2529  }
2530 
2531 
2532 
2540  template <int spacedim>
2541  static void
2543  const std::vector<types::global_dof_index> & /*new_numbers*/,
2544  const IndexSet & /*indices_we_care_about*/,
2545  DoFHandler<1, spacedim> & /*dof_handler*/,
2546  const unsigned int /*level*/,
2547  const bool /*check_validity*/)
2548  {
2549  // nothing to do in 1d because there are no separate faces
2550  }
2551 
2552 
2553 
2554  template <int spacedim>
2555  static void
2557  const std::vector<::types::global_dof_index> &new_numbers,
2558  const IndexSet & indices_we_care_about,
2559  DoFHandler<2, spacedim> &dof_handler,
2560  const unsigned int level,
2561  const bool check_validity)
2562  {
2563  if (dof_handler.get_fe().n_dofs_per_line() > 0)
2564  {
2565  // save user flags as they will be modified
2566  std::vector<bool> user_flags;
2567  dof_handler.get_triangulation().save_user_flags(user_flags);
2568  const_cast<::Triangulation<2, spacedim> &>(
2569  dof_handler.get_triangulation())
2570  .clear_user_flags();
2571 
2572  // flag all lines adjacent to cells of the current
2573  // level, as those lines logically belong to the same
2574  // level as the cell, at least for for isotropic
2575  // refinement
2576  for (const auto &cell :
2577  dof_handler.cell_iterators_on_level(level))
2578  if (cell->level_subdomain_id() !=
2580  for (const unsigned int line : cell->face_indices())
2581  cell->face(line)->set_user_flag();
2582 
2583  for (const auto &cell :
2584  dof_handler.cell_iterators_on_level(level))
2585  for (const auto l : cell->line_indices())
2586  if (cell->line(l)->user_flag_set())
2587  {
2588  for (unsigned int d = 0;
2589  d < dof_handler.get_fe().n_dofs_per_line();
2590  ++d)
2591  {
2593  cell->line(l)->mg_dof_index(level, d);
2594  if (check_validity)
2596  ExcInternalError());
2597 
2598  if (idx != numbers::invalid_dof_index)
2599  cell->line(l)->set_mg_dof_index(
2600  level,
2601  d,
2602  ((indices_we_care_about.size() == 0) ?
2603  new_numbers[idx] :
2604  new_numbers[indices_we_care_about
2605  .index_within_set(idx)]));
2606  }
2607  cell->line(l)->clear_user_flag();
2608  }
2609  // finally, restore user flags
2610  const_cast<::Triangulation<2, spacedim> &>(
2611  dof_handler.get_triangulation())
2612  .load_user_flags(user_flags);
2613  }
2614  }
2615 
2616 
2617 
2618  template <int spacedim>
2619  static void
2621  const std::vector<::types::global_dof_index> &new_numbers,
2622  const IndexSet & indices_we_care_about,
2623  DoFHandler<3, spacedim> &dof_handler,
2624  const unsigned int level,
2625  const bool check_validity)
2626  {
2627  if (dof_handler.get_fe().n_dofs_per_line() > 0 ||
2628  dof_handler.get_fe().max_dofs_per_quad() > 0)
2629  {
2630  // save user flags as they will be modified
2631  std::vector<bool> user_flags;
2632  dof_handler.get_triangulation().save_user_flags(user_flags);
2633  const_cast<::Triangulation<3, spacedim> &>(
2634  dof_handler.get_triangulation())
2635  .clear_user_flags();
2636 
2637  // flag all lines adjacent to cells of the current
2638  // level, as those lines logically belong to the same
2639  // level as the cell, at least for isotropic refinement
2640  for (const auto &cell :
2641  dof_handler.cell_iterators_on_level(level))
2642  if (cell->level_subdomain_id() !=
2644  for (const auto line : cell->line_indices())
2645  cell->line(line)->set_user_flag();
2646 
2647  for (const auto &cell :
2648  dof_handler.cell_iterators_on_level(level))
2649  for (const auto l : cell->line_indices())
2650  if (cell->line(l)->user_flag_set())
2651  {
2652  for (unsigned int d = 0;
2653  d < dof_handler.get_fe().n_dofs_per_line();
2654  ++d)
2655  {
2657  cell->line(l)->mg_dof_index(level, d);
2658  if (check_validity)
2660  ExcInternalError());
2661 
2662  if (idx != numbers::invalid_dof_index)
2663  cell->line(l)->set_mg_dof_index(
2664  level,
2665  d,
2666  ((indices_we_care_about.size() == 0) ?
2667  new_numbers[idx] :
2668  new_numbers[indices_we_care_about
2669  .index_within_set(idx)]));
2670  }
2671  cell->line(l)->clear_user_flag();
2672  }
2673 
2674  // flag all quads adjacent to cells of the current level, as
2675  // those quads logically belong to the same level as the cell,
2676  // at least for isotropic refinement
2677  for (const auto &cell :
2678  dof_handler.cell_iterators_on_level(level))
2679  if (cell->level_subdomain_id() !=
2681  for (const auto quad : cell->face_indices())
2682  cell->quad(quad)->set_user_flag();
2683 
2684  for (const auto &cell : dof_handler.cell_iterators())
2685  for (const auto l : cell->face_indices())
2686  if (cell->quad(l)->user_flag_set())
2687  {
2688  for (unsigned int d = 0;
2689  d < dof_handler.get_fe().n_dofs_per_quad(l);
2690  ++d)
2691  {
2693  cell->quad(l)->mg_dof_index(level, d);
2694  if (check_validity)
2696  ExcInternalError());
2697 
2698  if (idx != numbers::invalid_dof_index)
2699  cell->quad(l)->set_mg_dof_index(
2700  level,
2701  d,
2702  ((indices_we_care_about.size() == 0) ?
2703  new_numbers[idx] :
2704  new_numbers[indices_we_care_about
2705  .index_within_set(idx)]));
2706  }
2707  cell->quad(l)->clear_user_flag();
2708  }
2709 
2710  // finally, restore user flags
2711  const_cast<::Triangulation<3, spacedim> &>(
2712  dof_handler.get_triangulation())
2713  .load_user_flags(user_flags);
2714  }
2715  }
2716 
2717 
2718 
2719  template <int dim, int spacedim>
2720  static void
2722  const std::vector<::types::global_dof_index> &new_numbers,
2723  const IndexSet & indices_we_care_about,
2724  DoFHandler<dim, spacedim> &dof_handler,
2725  const unsigned int level,
2726  const bool check_validity)
2727  {
2728  Assert(
2729  dof_handler.hp_capability_enabled == false,
2731 
2732  Assert(level < dof_handler.get_triangulation().n_global_levels(),
2733  ExcInternalError());
2734 
2735  // renumber DoF indices on vertices, cells, and faces. this
2736  // can be done in parallel because the respective functions
2737  // work on separate data structures
2738  Threads::TaskGroup<> tasks;
2739  tasks += Threads::new_task([&]() {
2740  renumber_vertex_mg_dofs(new_numbers,
2741  indices_we_care_about,
2742  dof_handler,
2743  level,
2744  check_validity);
2745  });
2746  tasks += Threads::new_task([&]() {
2747  renumber_face_mg_dofs(new_numbers,
2748  indices_we_care_about,
2749  dof_handler,
2750  level,
2751  check_validity);
2752  });
2753  tasks += Threads::new_task([&]() {
2754  renumber_cell_mg_dofs(new_numbers,
2755  indices_we_care_about,
2756  dof_handler,
2757  level);
2758  });
2759  tasks.join_all();
2760  }
2761  };
2762 
2763 
2764 
2765  /* --------------------- class Sequential ---------------- */
2766 
2767 
2768 
2769  template <int dim, int spacedim>
2771  DoFHandler<dim, spacedim> &dof_handler)
2772  : dof_handler(&dof_handler)
2773  {}
2774 
2775 
2776 
2777  template <int dim, int spacedim>
2778  NumberCache
2780  {
2781  const types::global_dof_index n_initial_dofs =
2783  *dof_handler);
2784 
2785  const types::global_dof_index n_dofs =
2786  Implementation::unify_dof_indices(*dof_handler,
2787  n_initial_dofs,
2788  /*check_validity=*/true);
2789 
2790  // return a sequential, complete index set
2791  return NumberCache(n_dofs);
2792  }
2793 
2794 
2795 
2796  template <int dim, int spacedim>
2797  std::vector<NumberCache>
2799  {
2800  std::vector<bool> user_flags;
2801  dof_handler->get_triangulation().save_user_flags(user_flags);
2802 
2803  const_cast<::Triangulation<dim, spacedim> &>(
2804  dof_handler->get_triangulation())
2805  .clear_user_flags();
2806 
2807  std::vector<NumberCache> number_caches;
2808  number_caches.reserve(dof_handler->get_triangulation().n_levels());
2809  for (unsigned int level = 0;
2810  level < dof_handler->get_triangulation().n_levels();
2811  ++level)
2812  {
2813  // first distribute dofs on this level
2814  const types::global_dof_index n_level_dofs =
2816  numbers::invalid_subdomain_id, *dof_handler, level);
2817 
2818  // then add a complete, sequential index set
2819  number_caches.emplace_back(n_level_dofs);
2820  }
2821 
2822  const_cast<::Triangulation<dim, spacedim> &>(
2823  dof_handler->get_triangulation())
2824  .load_user_flags(user_flags);
2825 
2826  return number_caches;
2827  }
2828 
2829 
2830 
2831  template <int dim, int spacedim>
2832  NumberCache
2834  const std::vector<types::global_dof_index> &new_numbers) const
2835  {
2836  Implementation::renumber_dofs(new_numbers,
2837  IndexSet(0),
2838  *dof_handler,
2839  /*check_validity=*/true);
2840 
2841  // return a sequential, complete index set. take into account that the
2842  // number of DoF indices may in fact be smaller than there were before
2843  // if some previously separately numbered dofs have been identified.
2844  // this is, for example, what we do when the DoFHandler has hp-
2845  // capabilities enabled: it first enumerates all DoFs on cells
2846  // independently, and then unifies some located at vertices or faces;
2847  // this leaves us with fewer DoFs than there were before, so use the
2848  // largest index as the one to determine the size of the index space
2849  return NumberCache(
2850  *std::max_element(new_numbers.begin(), new_numbers.end()) + 1);
2851  }
2852 
2853 
2854 
2855  template <int dim, int spacedim>
2856  NumberCache
2858  const unsigned int level,
2859  const std::vector<types::global_dof_index> &new_numbers) const
2860  {
2862  new_numbers, IndexSet(0), *dof_handler, level, true);
2863 
2864  // return a sequential, complete index set
2865  return NumberCache(new_numbers.size());
2866  }
2867 
2868 
2869  /* --------------------- class ParallelShared ---------------- */
2870 
2871 
2872  template <int dim, int spacedim>
2874  DoFHandler<dim, spacedim> &dof_handler)
2875  : dof_handler(&dof_handler)
2876  {}
2877 
2878 
2879 
2880  namespace
2881  {
2890  template <int dim, int spacedim>
2891  std::vector<types::subdomain_id>
2892  get_dof_subdomain_association(
2893  const DoFHandler<dim, spacedim> &dof_handler,
2894  const types::global_dof_index n_dofs,
2895  const unsigned int n_procs)
2896  {
2897  (void)n_procs;
2898  std::vector<types::subdomain_id> subdomain_association(
2900  std::vector<types::global_dof_index> local_dof_indices;
2901  local_dof_indices.reserve(
2902  dof_handler.get_fe_collection().max_dofs_per_cell());
2903 
2904  // loop over all cells and record which subdomain a DoF belongs to.
2905  // give to the smaller subdomain_id in case it is on an interface
2906  for (const auto &cell : dof_handler.active_cell_iterators())
2907  {
2908  // get the owner of the cell; note that we have made sure above
2909  // that all cells are either locally owned or ghosts (not
2910  // artificial), so this call will always yield the true owner
2911  const types::subdomain_id subdomain_id = cell->subdomain_id();
2912  const unsigned int dofs_per_cell =
2913  cell->get_fe().n_dofs_per_cell();
2914  local_dof_indices.resize(dofs_per_cell);
2915  cell->get_dof_indices(local_dof_indices);
2916 
2917  // set subdomain ids. if dofs already have their values set then
2918  // they must be on partition interfaces. In that case assign them
2919  // to the processor with the smaller subdomain id.
2920  for (unsigned int i = 0; i < dofs_per_cell; ++i)
2921  if (subdomain_association[local_dof_indices[i]] ==
2923  subdomain_association[local_dof_indices[i]] = subdomain_id;
2924  else if (subdomain_association[local_dof_indices[i]] >
2925  subdomain_id)
2926  {
2927  subdomain_association[local_dof_indices[i]] = subdomain_id;
2928  }
2929  }
2930 
2931  Assert(std::find(subdomain_association.begin(),
2932  subdomain_association.end(),
2934  subdomain_association.end(),
2935  ExcInternalError());
2936 
2937  Assert(*std::max_element(subdomain_association.begin(),
2938  subdomain_association.end()) < n_procs,
2939  ExcInternalError());
2940 
2941  return subdomain_association;
2942  }
2943 
2944 
2951  template <int dim, int spacedim>
2952  std::vector<types::subdomain_id>
2953  get_dof_level_subdomain_association(
2954  const DoFHandler<dim, spacedim> &dof_handler,
2955  const types::global_dof_index n_dofs_on_level,
2956  const unsigned int n_procs,
2957  const unsigned int level)
2958  {
2959  (void)n_procs;
2960  std::vector<types::subdomain_id> level_subdomain_association(
2961  n_dofs_on_level, numbers::invalid_subdomain_id);
2962  std::vector<types::global_dof_index> local_dof_indices;
2963  local_dof_indices.reserve(
2964  dof_handler.get_fe_collection().max_dofs_per_cell());
2965 
2966  // loop over all cells and record which subdomain a DoF belongs to.
2967  // interface goes to proccessor with smaller subdomain id
2968  for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2969  {
2970  // get the owner of the cell; note that we have made sure above
2971  // that all cells are either locally owned or ghosts (not
2972  // artificial), so this call will always yield the true owner
2973  const types::subdomain_id level_subdomain_id =
2974  cell->level_subdomain_id();
2975  const unsigned int dofs_per_cell =
2976  cell->get_fe().n_dofs_per_cell();
2977  local_dof_indices.resize(dofs_per_cell);
2978  cell->get_mg_dof_indices(local_dof_indices);
2979 
2980  // set level subdomain ids. if dofs already have their values set
2981  // then they must be on partition interfaces. In that case assign
2982  // them to the processor with the smaller subdomain id.
2983  for (unsigned int i = 0; i < dofs_per_cell; ++i)
2984  if (level_subdomain_association[local_dof_indices[i]] ==
2986  level_subdomain_association[local_dof_indices[i]] =
2987  level_subdomain_id;
2988  else if (level_subdomain_association[local_dof_indices[i]] >
2989  level_subdomain_id)
2990  {
2991  level_subdomain_association[local_dof_indices[i]] =
2992  level_subdomain_id;
2993  }
2994  }
2995 
2996  Assert(std::find(level_subdomain_association.begin(),
2997  level_subdomain_association.end(),
2999  level_subdomain_association.end(),
3000  ExcInternalError());
3001 
3002  Assert(*std::max_element(level_subdomain_association.begin(),
3003  level_subdomain_association.end()) < n_procs,
3004  ExcInternalError());
3005 
3006  return level_subdomain_association;
3007  }
3008  } // namespace
3009 
3010 
3011 
3012  template <int dim, int spacedim>
3013  NumberCache
3015  {
3016  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3017  (dynamic_cast<
3018  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3019  &this->dof_handler->get_triangulation()));
3020  Assert(tr != nullptr, ExcInternalError());
3021 
3022  const unsigned int n_procs =
3023  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3024 
3025  // If an underlying shared::Tria allows artificial cells, we need to
3026  // restore the true cell owners temporarily.
3027  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we
3028  // save the current set of subdomain ids, set subdomain ids to the
3029  // "true" owner of each cell upon construction of the
3030  // TemporarilyRestoreSubdomainIds object, and later restore these flags
3031  // when it is destroyed.
3032  const internal::parallel::shared::
3033  TemporarilyRestoreSubdomainIds<dim, spacedim>
3034  subdomain_modifier(*tr);
3035 
3036  // first let the sequential algorithm do its magic. it is going to
3037  // enumerate DoFs on all cells, regardless of owner
3038  const types::global_dof_index n_initial_dofs =
3040  *this->dof_handler);
3041 
3042  const types::global_dof_index n_dofs =
3043  Implementation::unify_dof_indices(*this->dof_handler,
3044  n_initial_dofs,
3045  /*check_validity=*/true);
3046 
3047  // then re-enumerate them based on their subdomain association.
3048  // for this, we first have to identify for each current DoF
3049  // index which subdomain they belong to. ideally, we would
3050  // like to call DoFRenumbering::subdomain_wise(), but
3051  // because the NumberCache of the current DoFHandler is not
3052  // fully set up yet, we can't quite do that. also, that
3053  // function has to deal with other kinds of triangulations as
3054  // well, whereas we here know what kind of triangulation
3055  // we have and can simplify the code accordingly
3056  std::vector<types::global_dof_index> new_dof_indices(
3057  n_dofs, enumeration_dof_index);
3058  {
3059  // first get the association of each dof with a subdomain and
3060  // determine the total number of subdomain ids used
3061  const std::vector<types::subdomain_id> subdomain_association =
3062  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3063 
3064  // then renumber the subdomains by first looking at those belonging
3065  // to subdomain 0, then those of subdomain 1, etc. note that the
3066  // algorithm is stable, i.e. if two dofs i,j have i<j and belong to
3067  // the same subdomain, then they will be in this order also after
3068  // reordering
3069  types::global_dof_index next_free_index = 0;
3070  for (types::subdomain_id subdomain = 0; subdomain < n_procs;
3071  ++subdomain)
3072  for (types::global_dof_index i = 0; i < n_dofs; ++i)
3073  if (subdomain_association[i] == subdomain)
3074  {
3075  Assert(new_dof_indices[i] == enumeration_dof_index,
3076  ExcInternalError());
3077  new_dof_indices[i] = next_free_index;
3078  ++next_free_index;
3079  }
3080 
3081  // we should have numbered all dofs
3082  Assert(next_free_index == n_dofs, ExcInternalError());
3083  Assert(std::find(new_dof_indices.begin(),
3084  new_dof_indices.end(),
3085  enumeration_dof_index) == new_dof_indices.end(),
3086  ExcInternalError());
3087  }
3088  // finally do the renumbering. we can use the sequential
3089  // version of the function because we do things on all
3090  // cells and all cells have their subdomain ids and DoFs
3091  // correctly set
3092  Implementation::renumber_dofs(new_dof_indices,
3093  IndexSet(0),
3094  *this->dof_handler,
3095  /*check_validity=*/true);
3096 
3097  // update the number cache. for this, we first have to find the
3098  // subdomain association for each DoF again following renumbering, from
3099  // which we can then compute the IndexSets of locally owned DoFs for all
3100  // processors. all other fields then follow from this
3101  //
3102  // given the way we enumerate degrees of freedom, the locally owned
3103  // ranges must all be contiguous and consecutive. this makes filling
3104  // the IndexSets cheap. an assertion at the top verifies that this
3105  // assumption is true
3106  const std::vector<types::subdomain_id> subdomain_association =
3107  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3108 
3109  for (unsigned int i = 1; i < n_dofs; ++i)
3110  Assert(subdomain_association[i] >= subdomain_association[i - 1],
3111  ExcInternalError());
3112 
3113  std::vector<IndexSet> locally_owned_dofs_per_processor(
3114  n_procs, IndexSet(n_dofs));
3115  {
3116  // we know that the set of subdomain indices is contiguous from
3117  // the assertion above; find the start and end index for each
3118  // processor, taking into account that sometimes a processor
3119  // may not in fact have any DoFs at all. we do the latter
3120  // by just identifying contiguous ranges of subdomain_ids
3121  // and filling IndexSets for those subdomains; subdomains
3122  // that don't appear will lead to IndexSets that are simply
3123  // never touched and remain empty as initialized above.
3124  unsigned int start_index = 0;
3125  unsigned int end_index = 0;
3126  while (start_index < n_dofs)
3127  {
3128  while ((end_index) < n_dofs &&
3129  (subdomain_association[end_index] ==
3130  subdomain_association[start_index]))
3131  ++end_index;
3132 
3133  // we've now identified a range of same indices. set that
3134  // range in the corresponding IndexSet
3135  if (end_index > start_index)
3136  {
3137  const unsigned int subdomain_owner =
3138  subdomain_association[start_index];
3139  locally_owned_dofs_per_processor[subdomain_owner].add_range(
3140  start_index, end_index);
3141  }
3142 
3143  // then move on to thinking about the next range
3144  start_index = end_index;
3145  }
3146  }
3147 
3148  // return a NumberCache object made up from the sets of locally
3149  // owned DoFs
3150  return NumberCache(
3151  locally_owned_dofs_per_processor,
3152  this->dof_handler->get_triangulation().locally_owned_subdomain());
3153  }
3154 
3155 
3156 
3157  template <int dim, int spacedim>
3158  std::vector<NumberCache>
3160  {
3161  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3162  (dynamic_cast<
3163  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3164  &this->dof_handler->get_triangulation()));
3165  Assert(tr != nullptr, ExcInternalError());
3166 
3167  const unsigned int n_procs =
3168  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3169  const unsigned int n_levels = tr->n_global_levels();
3170 
3171  std::vector<NumberCache> number_caches;
3172  number_caches.reserve(n_levels);
3173 
3174  // We create an index set for each level
3175  for (unsigned int lvl = 0; lvl < n_levels; ++lvl)
3176  {
3177  // If the underlying shared::Tria allows artificial cells,
3178  // then save the current set of level subdomain ids, and set
3179  // subdomain ids to the "true" owner of each cell. we later
3180  // restore these flags
3181  // Note: "allows_artificial_cells" is currently enforced for
3182  // MG computations.
3183  std::vector<types::subdomain_id> saved_level_subdomain_ids;
3184  saved_level_subdomain_ids.resize(tr->n_cells(lvl));
3185  {
3186  typename ::parallel::shared::Triangulation<dim, spacedim>::
3187  cell_iterator cell =
3188  this->dof_handler->get_triangulation().begin(
3189  lvl),
3190  endc =
3191  this->dof_handler->get_triangulation().end(lvl);
3192 
3193  const std::vector<types::subdomain_id> &true_level_subdomain_ids =
3194  tr->get_true_level_subdomain_ids_of_cells(lvl);
3195 
3196  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3197  {
3198  saved_level_subdomain_ids[index] = cell->level_subdomain_id();
3199  cell->set_level_subdomain_id(true_level_subdomain_ids[index]);
3200  }
3201  }
3202 
3203  // Next let the sequential algorithm do its magic. it is going to
3204  // enumerate DoFs on all cells on the given level, regardless of
3205  // owner
3206  const types::global_dof_index n_dofs_on_level =
3208  numbers::invalid_subdomain_id, *this->dof_handler, lvl);
3209 
3210  // then re-enumerate them based on their level subdomain
3211  // association. for this, we first have to identify for each current
3212  // DoF index which subdomain they belong to. ideally, we would like
3213  // to call DoFRenumbering::subdomain_wise(), but because the
3214  // NumberCache of the current DoFHandler is not fully set up yet, we
3215  // can't quite do that. also, that function has to deal with other
3216  // kinds of triangulations as well, whereas we here know what kind
3217  // of triangulation we have and can simplify the code accordingly
3218  std::vector<types::global_dof_index> new_dof_indices(
3219  n_dofs_on_level, numbers::invalid_dof_index);
3220  {
3221  // first get the association of each dof with a subdomain and
3222  // determine the total number of subdomain ids used
3223  const std::vector<types::subdomain_id>
3224  level_subdomain_association =
3225  get_dof_level_subdomain_association(*this->dof_handler,
3226  n_dofs_on_level,
3227  n_procs,
3228  lvl);
3229 
3230  // then renumber the subdomains by first looking at those
3231  // belonging to subdomain 0, then those of subdomain 1, etc. note
3232  // that the algorithm is stable, i.e. if two dofs i,j have i<j and
3233  // belong to the same subdomain, then they will be in this order
3234  // also after reordering
3235  types::global_dof_index next_free_index = 0;
3236  for (types::subdomain_id level_subdomain = 0;
3237  level_subdomain < n_procs;
3238  ++level_subdomain)
3239  for (types::global_dof_index i = 0; i < n_dofs_on_level; ++i)
3240  if (level_subdomain_association[i] == level_subdomain)
3241  {
3242  Assert(new_dof_indices[i] == numbers::invalid_dof_index,
3243  ExcInternalError());
3244  new_dof_indices[i] = next_free_index;
3245  ++next_free_index;
3246  }
3247 
3248  // we should have numbered all dofs
3249  Assert(next_free_index == n_dofs_on_level, ExcInternalError());
3250  Assert(std::find(new_dof_indices.begin(),
3251  new_dof_indices.end(),
3253  new_dof_indices.end(),
3254  ExcInternalError());
3255  }
3256 
3257  // finally do the renumbering. we can use the sequential
3258  // version of the function because we do things on all
3259  // cells and all cells have their subdomain ids and DoFs
3260  // correctly set
3262  new_dof_indices, IndexSet(0), *this->dof_handler, lvl, true);
3263 
3264  // update the number cache. for this, we first have to find the
3265  // level subdomain association for each DoF again following
3266  // renumbering, from which we can then compute the IndexSets of
3267  // locally owned DoFs for all processors. all other fields then
3268  // follow from this
3269  //
3270  // given the way we enumerate degrees of freedom, the locally owned
3271  // ranges must all be contiguous and consecutive. this makes filling
3272  // the IndexSets cheap. an assertion at the top verifies that this
3273  // assumption is true
3274  const std::vector<types::subdomain_id> level_subdomain_association =
3275  get_dof_level_subdomain_association(*this->dof_handler,
3276  n_dofs_on_level,
3277  n_procs,
3278  lvl);
3279 
3280  for (unsigned int i = 1; i < n_dofs_on_level; ++i)
3281  Assert(level_subdomain_association[i] >=
3282  level_subdomain_association[i - 1],
3283  ExcInternalError());
3284 
3285  std::vector<IndexSet> locally_owned_dofs_per_processor(
3286  n_procs, IndexSet(n_dofs_on_level));
3287  {
3288  // we know that the set of subdomain indices is contiguous from
3289  // the assertion above; find the start and end index for each
3290  // processor, taking into account that sometimes a processor
3291  // may not in fact have any DoFs at all. we do the latter
3292  // by just identifying contiguous ranges of level_subdomain_ids
3293  // and filling IndexSets for those subdomains; subdomains
3294  // that don't appear will lead to IndexSets that are simply
3295  // never touched and remain empty as initialized above.
3296  unsigned int start_index = 0;
3297  unsigned int end_index = 0;
3298  while (start_index < n_dofs_on_level)
3299  {
3300  while ((end_index) < n_dofs_on_level &&
3301  (level_subdomain_association[end_index] ==
3302  level_subdomain_association[start_index]))
3303  ++end_index;
3304 
3305  // we've now identified a range of same indices. set that
3306  // range in the corresponding IndexSet
3307  if (end_index > start_index)
3308  {
3309  const unsigned int level_subdomain_owner =
3310  level_subdomain_association[start_index];
3311  locally_owned_dofs_per_processor[level_subdomain_owner]
3312  .add_range(start_index, end_index);
3313  }
3314 
3315  // then move on to thinking about the next range
3316  start_index = end_index;
3317  }
3318  }
3319 
3320  // finally, restore current level subdomain ids
3321  {
3322  typename ::parallel::shared::Triangulation<dim, spacedim>::
3323  cell_iterator cell =
3324  this->dof_handler->get_triangulation().begin(
3325  lvl),
3326  endc =
3327  this->dof_handler->get_triangulation().end(lvl);
3328 
3329  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3330  cell->set_level_subdomain_id(saved_level_subdomain_ids[index]);
3331 
3332  // add NumberCache for current level
3333  number_caches.emplace_back(
3334  NumberCache(locally_owned_dofs_per_processor,
3335  this->dof_handler->get_triangulation()
3336  .locally_owned_subdomain()));
3337  }
3338  }
3339 
3340  return number_caches;
3341  }
3342 
3343 
3344 
3345  template <int dim, int spacedim>
3346  NumberCache
3348  const std::vector<types::global_dof_index> &new_numbers) const
3349  {
3350 #ifndef DEAL_II_WITH_MPI
3351  (void)new_numbers;
3352  Assert(false, ExcNotImplemented());
3353  return NumberCache();
3354 #else
3355  // Similar to distribute_dofs() we need to have a special treatment in
3356  // case artificial cells are present.
3357  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3358  (dynamic_cast<
3359  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3360  &this->dof_handler->get_triangulation()));
3361  Assert(tr != nullptr, ExcInternalError());
3362 
3363  // Set subdomain IDs to the "true" owner of each cell.
3364  const internal::parallel::shared::
3365  TemporarilyRestoreSubdomainIds<dim, spacedim>
3366  subdomain_modifier(*tr);
3367 
3368  std::vector<types::global_dof_index> global_gathered_numbers(
3369  this->dof_handler->n_dofs(), 0);
3370  // as we call DoFRenumbering::subdomain_wise (*dof_handler) from
3371  // distribute_dofs(), we need to support sequential-like input.
3372  // Distributed-like input from, for example, component_wise renumbering
3373  // is also supported.
3374  if (new_numbers.size() == this->dof_handler->n_dofs())
3375  {
3376  global_gathered_numbers = new_numbers;
3377  }
3378  else
3379  {
3380  Assert(new_numbers.size() ==
3381  this->dof_handler->locally_owned_dofs().n_elements(),
3382  ExcInternalError());
3383  const unsigned int n_cpu =
3384  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3385  std::vector<types::global_dof_index> gathered_new_numbers(
3386  this->dof_handler->n_dofs(), 0);
3387  Assert(Utilities::MPI::this_mpi_process(tr->get_communicator()) ==
3388  this->dof_handler->get_triangulation()
3389  .locally_owned_subdomain(),
3390  ExcInternalError())
3391 
3392  // gather new numbers among processors into one vector
3393  {
3394  std::vector<types::global_dof_index> new_numbers_copy(
3395  new_numbers);
3396 
3397  // store the number of elements that are to be received from each
3398  // process
3399  std::vector<int> rcounts(n_cpu);
3400 
3402  // set rcounts based on new_numbers:
3403  int cur_count = new_numbers_copy.size();
3404  int ierr = MPI_Allgather(&cur_count,
3405  1,
3406  MPI_INT,
3407  rcounts.data(),
3408  1,
3409  MPI_INT,
3410  tr->get_communicator());
3411  AssertThrowMPI(ierr);
3412 
3413  // compute the displacements (relative to recvbuf)
3414  // at which to place the incoming data from process i
3415  std::vector<int> displacements(n_cpu);
3416  for (unsigned int i = 0; i < n_cpu; i++)
3417  {
3418  displacements[i] = shift;
3419  shift += rcounts[i];
3420  }
3421  Assert(new_numbers_copy.size() ==
3422  static_cast<unsigned int>(
3424  tr->get_communicator())]),
3425  ExcInternalError());
3426  ierr = MPI_Allgatherv(new_numbers_copy.data(),
3427  new_numbers_copy.size(),
3429  gathered_new_numbers.data(),
3430  rcounts.data(),
3431  displacements.data(),
3433  tr->get_communicator());
3434  AssertThrowMPI(ierr);
3435  }
3436 
3437  // put new numbers according to the current
3438  // locally_owned_dofs_per_processor IndexSets
3440  // flag_1 and flag_2 are
3441  // used to control that there is a
3442  // one-to-one relation between old and new DoFs.
3443  std::vector<unsigned int> flag_1(this->dof_handler->n_dofs(), 0);
3444  std::vector<unsigned int> flag_2(this->dof_handler->n_dofs(), 0);
3445  for (unsigned int i = 0; i < n_cpu; i++)
3446  {
3447  const IndexSet iset =
3448  this->dof_handler->locally_owned_dofs_per_processor()[i];
3449  for (types::global_dof_index ind = 0; ind < iset.n_elements();
3450  ind++)
3451  {
3452  const types::global_dof_index target =
3453  iset.nth_index_in_set(ind);
3455  gathered_new_numbers[shift + ind];
3456  Assert(target < this->dof_handler->n_dofs(),
3457  ExcInternalError());
3458  Assert(value < this->dof_handler->n_dofs(),
3459  ExcInternalError());
3460  global_gathered_numbers[target] = value;
3461  flag_1[target]++;
3462  flag_2[value]++;
3463  }
3464  shift += iset.n_elements();
3465  }
3466 
3467  Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1,
3468  ExcInternalError());
3469  Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1,
3470  ExcInternalError());
3471  Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1,
3472  ExcInternalError());
3473  Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1,
3474  ExcInternalError());
3475  }
3476 
3477  // let the sequential algorithm do its magic; ignore the
3478  // return type, but reconstruct the number cache based on
3479  // which DoFs each process owns
3480  Implementation::renumber_dofs(global_gathered_numbers,
3481  IndexSet(0),
3482  *this->dof_handler,
3483  /*check_validity=*/true);
3484 
3485  const NumberCache number_cache(
3486  DoFTools::locally_owned_dofs_per_subdomain(*this->dof_handler),
3487  this->dof_handler->get_triangulation().locally_owned_subdomain());
3488 
3489  return number_cache;
3490 #endif
3491  }
3492 
3493 
3494 
3495  template <int dim, int spacedim>
3496  NumberCache
3498  const unsigned int /*level*/,
3499  const std::vector<types::global_dof_index> & /*new_numbers*/) const
3500  {
3501  // multigrid is not currently implemented for shared triangulations
3502  Assert(false, ExcNotImplemented());
3503 
3504  return {};
3505  }
3506 
3507 
3508 
3509  /* --------------------- class ParallelDistributed ---------------- */
3510 
3511 #ifdef DEAL_II_WITH_MPI
3512 
3513  namespace
3514  {
3515  template <int dim, int spacedim>
3516  void
3517  communicate_mg_ghost_cells(
3518  const typename ::parallel::
3519  DistributedTriangulationBase<dim, spacedim> &tria,
3520  DoFHandler<dim, spacedim> & dof_handler)
3521  {
3522  (void)tria;
3523  const auto pack = [&](const auto &cell) {
3524  // why would somebody request a cell that is not ours?
3525  Assert(cell->level_subdomain_id() == tria.locally_owned_subdomain(),
3526  ExcInternalError());
3527 
3528  std::vector<::types::global_dof_index> data(
3529  cell->get_fe().n_dofs_per_cell());
3530  cell->get_mg_dof_indices(data);
3531 
3532  return data;
3533  };
3534 
3535  const auto unpack = [](const auto &cell, const auto &dofs) {
3536  Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3537  ExcInternalError());
3538 
3539  Assert(cell->level_subdomain_id() !=
3541  ExcInternalError());
3542 
3543  std::vector<::types::global_dof_index> dof_indices(
3544  cell->get_fe().n_dofs_per_cell());
3545  cell->get_mg_dof_indices(dof_indices);
3546 
3547  bool complete = true;
3548  for (unsigned int i = 0; i < dof_indices.size(); ++i)
3549  if (dofs[i] != numbers::invalid_dof_index)
3550  {
3551  Assert((dof_indices[i] == (numbers::invalid_dof_index)) ||
3552  (dof_indices[i] == dofs[i]),
3553  ExcInternalError());
3554  dof_indices[i] = dofs[i];
3555  }
3556  else
3557  complete = false;
3558 
3559  if (!complete)
3560  const_cast<
3562  ->set_user_flag();
3563  else
3564  const_cast<
3566  ->clear_user_flag();
3567 
3568  const_cast<
3570  ->set_mg_dof_indices(dof_indices);
3571  };
3572 
3573  const auto filter = [](const auto &cell) {
3574  return cell->user_flag_set();
3575  };
3576 
3578  std::vector<types::global_dof_index>,
3580  }
3581 
3582 
3583 
3584  template <int spacedim>
3585  void
3586  communicate_mg_ghost_cells(const typename ::parallel::
3587  distributed::Triangulation<1, spacedim> &,
3589  {
3590  Assert(false, ExcNotImplemented());
3591  }
3592 
3593 
3594 
3613  template <int dim, int spacedim>
3614  void
3615  communicate_dof_indices_on_marked_cells(
3616  const DoFHandler<dim, spacedim> &dof_handler)
3617  {
3618 # ifndef DEAL_II_WITH_MPI
3619  (void)dof_handler;
3620  Assert(false, ExcNotImplemented());
3621 # else
3622 
3623  // define functions that pack data on cells that are ghost cells
3624  // somewhere else, and unpack data on cells where we get information
3625  // from elsewhere
3626  const auto pack = [](const auto &cell) {
3627  Assert(cell->is_locally_owned(), ExcInternalError());
3628 
3629  std::vector<::types::global_dof_index> data(
3630  cell->get_fe().n_dofs_per_cell());
3631  cell->get_dof_indices(data);
3632 
3633  return data;
3634  };
3635 
3636  const auto unpack = [](const auto &cell, const auto &dofs) {
3637  Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3638  ExcInternalError());
3639 
3640  Assert(cell->is_ghost(), ExcInternalError());
3641 
3642  std::vector<::types::global_dof_index> dof_indices(
3643  cell->get_fe().n_dofs_per_cell());
3644  cell->update_cell_dof_indices_cache();
3645  cell->get_dof_indices(dof_indices);
3646 
3647  bool complete = true;
3648  for (unsigned int i = 0; i < dof_indices.size(); ++i)
3649  if (dofs[i] != numbers::invalid_dof_index)
3650  {
3651  Assert((dof_indices[i] == (numbers::invalid_dof_index)) ||
3652  (dof_indices[i] == dofs[i]),
3653  ExcInternalError());
3654  dof_indices[i] = dofs[i];
3655  }
3656  else
3657  complete = false;
3658 
3659  if (!complete)
3660  const_cast<
3662  cell)
3663  ->set_user_flag();
3664  else
3665  const_cast<
3667  cell)
3668  ->clear_user_flag();
3669 
3670  const_cast<
3672  ->set_dof_indices(dof_indices);
3673  };
3674 
3675  const auto filter = [](const auto &cell) {
3676  return cell->user_flag_set();
3677  };
3678 
3680  std::vector<types::global_dof_index>,
3682 
3683  // finally update the cell DoF indices caches to make sure
3684  // our internal data structures are consistent
3685  update_all_active_cell_dof_indices_caches(dof_handler);
3686 
3687 
3688  // have a barrier so that sends between two calls to this
3689  // function are not mixed up.
3690  //
3691  // this is necessary because above we just see if there are
3692  // messages and then receive them, without discriminating
3693  // where they come from and whether they were sent in phase
3694  // 1 or 2 (the function is called twice in a row). the need
3695  // for a global communication step like this barrier could
3696  // be avoided by receiving messages specifically from those
3697  // processors from which we expect messages, and by using
3698  // different tags for phase 1 and 2, but the cost of a
3699  // barrier is negligible compared to everything else we do
3700  // here
3701  if (const auto *triangulation =
3702  dynamic_cast<const ::parallel::
3703  DistributedTriangulationBase<dim, spacedim> *>(
3704  &dof_handler.get_triangulation()))
3705  {
3706  const int ierr = MPI_Barrier(triangulation->get_communicator());
3707  AssertThrowMPI(ierr);
3708  }
3709  else
3710  {
3711  Assert(false,
3712  ExcMessage(
3713  "The function communicate_dof_indices_on_marked_cells() "
3714  "only works with parallel distributed triangulations."));
3715  }
3716 # endif
3717  }
3718 
3719 
3720 
3721  } // namespace
3722 
3723 #endif // DEAL_II_WITH_MPI
3724 
3725 
3726 
3727  template <int dim, int spacedim>
3729  DoFHandler<dim, spacedim> &dof_handler)
3730  : dof_handler(&dof_handler)
3731  {}
3732 
3733 
3734 
3735  template <int dim, int spacedim>
3736  NumberCache
3738  {
3739 #ifndef DEAL_II_WITH_MPI
3740  Assert(false, ExcNotImplemented());
3741  return NumberCache();
3742 #else
3743 
3745  *triangulation =
3746  (dynamic_cast<
3748  const_cast<::Triangulation<dim, spacedim> *>(
3749  &dof_handler->get_triangulation())));
3750  Assert(triangulation != nullptr, ExcInternalError());
3751 
3753  triangulation->locally_owned_subdomain();
3754 
3755 
3756  /*
3757  The following algorithm has a number of stages that are all
3758  documented in the paper that describes the parallel::distributed
3759  functionality:
3760 
3761  1/ locally enumerate dofs on locally owned cells
3762  2/ eliminate dof duplicates on all cells.
3763  un-numerate those that are on interfaces with ghost
3764  cells and that we don't own based on the tie-breaking
3765  criterion. unify dofs afterwards.
3766  3/ unify dofs and re-enumerate the remaining valid ones.
3767  the end result is that we only enumerate locally owned
3768  DoFs
3769  4/ shift indices so that each processor has a unique
3770  range of indices
3771  5/ for all locally owned cells that are ghost
3772  cells somewhere else, send our own DoF indices
3773  to the appropriate set of other processors.
3774  overwrite invalid DoF indices on ghost interfaces
3775  with the corresponding valid ones that we now know.
3776  6/ send DoF indices again to get the correct indices
3777  on ghost cells that we may not have known earlier
3778  */
3779 
3780  // --------- Phase 1: enumerate dofs on locally owned cells
3781  const types::global_dof_index n_initial_local_dofs =
3782  Implementation::distribute_dofs(subdomain_id, *dof_handler);
3783 
3784  // --------- Phase 2: eliminate dof duplicates on all cells:
3785  // - un-numerate dofs on interfaces to ghost cells
3786  // that we don't own
3787  // - in case of hp-support, unify dofs
3788  std::vector<::types::global_dof_index> renumbering(
3789  n_initial_local_dofs, enumeration_dof_index);
3790 
3791  // first, we invalidate degrees of freedom that belong to processors
3792  // of a lower rank, from which we will receive the final (and lower)
3793  // degrees of freedom later.
3796  renumbering, subdomain_id, *dof_handler);
3797 
3798  // then, we identify DoF duplicates if the DoFHandler has hp-
3799  // capabilities
3800  std::vector<std::map<types::global_dof_index, types::global_dof_index>>
3801  all_constrained_indices(dim);
3802  Implementation::compute_dof_identities(all_constrained_indices,
3803  *dof_handler);
3804 
3805  // --------- Phase 3: re-enumerate the valid degrees of freedom
3806  // consecutively. thus, we finally receive the
3807  // correct number of locally owned DoFs after
3808  // this step.
3809  //
3810  // the order in which we handle Phases 2 and 3 is important,
3811  // since we want to clarify ownership of degrees of freedom before
3812  // we actually unify and enumerate their indices. otherwise, we could
3813  // end up having a degee of freedom to which only invalid indices will
3814  // be assigned.
3815  const types::global_dof_index n_locally_owned_dofs =
3817  renumbering, all_constrained_indices, *dof_handler);
3818 
3819  // --------- Phase 4: shift indices so that each processor has a unique
3820  // range of indices
3821  ::types::global_dof_index my_shift = 0;
3822  const int ierr =
3823  MPI_Exscan(DEAL_II_MPI_CONST_CAST(&n_locally_owned_dofs),
3824  &my_shift,
3825  1,
3827  MPI_SUM,
3828  triangulation->get_communicator());
3829  AssertThrowMPI(ierr);
3830 
3831  // make dof indices globally consecutive
3832  for (auto &new_index : renumbering)
3833  if (new_index != numbers::invalid_dof_index)
3834  new_index += my_shift;
3835 
3836  // now re-enumerate all dofs to this shifted and condensed
3837  // numbering form. we renumber some dofs as invalid, so
3838  // choose the nocheck-version.
3839  Implementation::renumber_dofs(renumbering,
3840  IndexSet(0),
3841  *dof_handler,
3842  /*check_validity=*/false);
3843 
3844  // now a little bit of housekeeping
3845  const ::types::global_dof_index n_global_dofs =
3846  Utilities::MPI::sum(n_locally_owned_dofs,
3847  triangulation->get_communicator());
3848 
3849  NumberCache number_cache;
3850  number_cache.n_global_dofs = n_global_dofs;
3851  number_cache.n_locally_owned_dofs = n_locally_owned_dofs;
3852  number_cache.locally_owned_dofs = IndexSet(n_global_dofs);
3853  number_cache.locally_owned_dofs.add_range(my_shift,
3854  my_shift +
3855  n_locally_owned_dofs);
3856  number_cache.locally_owned_dofs.compress();
3857 
3858  // this ends the phase where we enumerate degrees of freedom on
3859  // each processor. what is missing is communicating DoF indices
3860  // on ghost cells
3861 
3862  // --------- Phase 5: for all locally owned cells that are ghost
3863  // cells somewhere else, send our own DoF indices
3864  // to the appropriate set of other processors
3865  {
3866  std::vector<bool> user_flags;
3867  triangulation->save_user_flags(user_flags);
3868  triangulation->clear_user_flags();
3869 
3870  // mark all cells that either have to send data (locally
3871  // owned cells that are adjacent to ghost neighbors in some
3872  // way) or receive data (all ghost cells) via the user flags
3873  for (const auto &cell : dof_handler->active_cell_iterators())
3874  if (cell->is_ghost())
3875  cell->set_user_flag();
3876 
3877  // Send and receive cells. After this, only the local cells
3878  // are marked, that received new data. This has to be
3879  // communicated in a second communication step.
3880  //
3881  // as explained in the 'distributed' paper, this has to be
3882  // done twice
3883  communicate_dof_indices_on_marked_cells(*dof_handler);
3884 
3885  // If the DoFHandler has hp-capabilities enabled, then we may have
3886  // received valid indices of degrees of freedom that are dominated
3887  // by a FE object adjacent to a ghost interface. thus, we overwrite
3888  // the remaining invalid indices with the valid ones in this step.
3890  *dof_handler);
3891 
3892  // --------- Phase 6: all locally owned cells have their correct
3893  // DoF indices set. however, some ghost cells
3894  // may still have invalid ones. thus, exchange
3895  // one more time.
3896  communicate_dof_indices_on_marked_cells(*dof_handler);
3897 
3898  // at this point, we must have taken care of the data transfer
3899  // on all cells we had previously marked. verify this
3900 # ifdef DEBUG
3901  for (const auto &cell : dof_handler->active_cell_iterators())
3902  Assert(cell->user_flag_set() == false, ExcInternalError());
3903 # endif
3904 
3905  triangulation->load_user_flags(user_flags);
3906  }
3907 
3908 # ifdef DEBUG
3909  // check that we are really done
3910  {
3911  std::vector<::types::global_dof_index> local_dof_indices;
3912 
3913  for (const auto &cell : dof_handler->active_cell_iterators())
3914  if (!cell->is_artificial())
3915  {
3916  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
3917  cell->get_dof_indices(local_dof_indices);
3918  if (local_dof_indices.end() !=
3919  std::find(local_dof_indices.begin(),
3920  local_dof_indices.end(),
3922  {
3923  if (cell->is_ghost())
3924  {
3925  Assert(false,
3926  ExcMessage(
3927  "A ghost cell ended up with incomplete "
3928  "DoF index information. This should not "
3929  "have happened!"));
3930  }
3931  else
3932  {
3933  Assert(
3934  false,
3935  ExcMessage(
3936  "A locally owned cell ended up with incomplete "
3937  "DoF index information. This should not "
3938  "have happened!"));
3939  }
3940  }
3941  }
3942  }
3943 # endif // DEBUG
3944  return number_cache;
3945 #endif // DEAL_II_WITH_MPI
3946  }
3947 
3948 
3949 
3950  template <int dim, int spacedim>
3951  std::vector<NumberCache>
3953  {
3954 #ifndef DEAL_II_WITH_MPI
3955  Assert(false, ExcNotImplemented());
3956  return std::vector<NumberCache>();
3957 #else
3958 
3960  *triangulation =
3961  (dynamic_cast<
3963  const_cast<::Triangulation<dim, spacedim> *>(
3964  &dof_handler->get_triangulation())));
3965  Assert(triangulation != nullptr, ExcInternalError());
3966 
3967  AssertThrow((triangulation->is_multilevel_hierarchy_constructed()),
3968  ExcMessage(
3969  "Multigrid DoFs can only be distributed on a parallel "
3970  "Triangulation if the flag construct_multigrid_hierarchy "
3971  "is set in the constructor."));
3972 
3973  // loop over all levels that exist globally (across all
3974  // processors), even if the current processor does not in fact
3975  // have any cells on that level or if the local part of the
3976  // Triangulation has fewer levels. we need to do this because
3977  // we need to communicate across all processors on all levels
3978  const unsigned int n_levels = triangulation->n_global_levels();
3979  std::vector<NumberCache> number_caches;
3980  number_caches.reserve(n_levels);
3981  for (unsigned int level = 0; level < n_levels; ++level)
3982  {
3983  NumberCache level_number_cache;
3984 
3985  //* 1. distribute on own subdomain
3986  const unsigned int n_initial_local_dofs =
3988  triangulation->locally_owned_subdomain(), *dof_handler, level);
3989 
3990  //* 2. iterate over ghostcells and kill dofs that are not
3991  // owned by us
3992  std::vector<::types::global_dof_index> renumbering(
3993  n_initial_local_dofs);
3994  for (::types::global_dof_index i = 0; i < renumbering.size();
3995  ++i)
3996  renumbering[i] = i;
3997 
3998  if (level < triangulation->n_levels())
3999  {
4000  std::vector<::types::global_dof_index> local_dof_indices;
4001 
4002  for (const auto &cell :
4003  dof_handler->cell_iterators_on_level(level))
4004  if (cell->level_subdomain_id() !=
4006  (cell->level_subdomain_id() <
4007  triangulation->locally_owned_subdomain()))
4008  {
4009  // we found a neighboring ghost cell whose
4010  // subdomain is "stronger" than our own
4011  // subdomain
4012 
4013  // delete all dofs that live there and that we
4014  // have previously assigned a number to
4015  // (i.e. the ones on the interface)
4016  local_dof_indices.resize(
4017  cell->get_fe().n_dofs_per_cell());
4018  cell->get_mg_dof_indices(local_dof_indices);
4019  for (unsigned int i = 0;
4020  i < cell->get_fe().n_dofs_per_cell();
4021  ++i)
4022  if (local_dof_indices[i] != numbers::invalid_dof_index)
4023  renumbering[local_dof_indices[i]] =
4025  }
4026  }
4027 
4028  level_number_cache.n_locally_owned_dofs = 0;
4029  for (types::global_dof_index &index : renumbering)
4030  if (index != numbers::invalid_dof_index)
4031  index = level_number_cache.n_locally_owned_dofs++;
4032 
4033  //* 3. communicate local dofcount and shift ids to make
4034  // them unique
4035  ::types::global_dof_index my_shift = 0;
4036  int ierr = MPI_Exscan(DEAL_II_MPI_CONST_CAST(
4037  &level_number_cache.n_locally_owned_dofs),
4038  &my_shift,
4039  1,
4041  MPI_SUM,
4042  triangulation->get_communicator());
4043  AssertThrowMPI(ierr);
4044 
4045  // The last processor knows about the total number of dofs, so we
4046  // can use a cheaper broadcast rather than an MPI_Allreduce via
4047  // MPI::sum().
4048  level_number_cache.n_global_dofs =
4049  my_shift + level_number_cache.n_locally_owned_dofs;
4050  ierr = MPI_Bcast(&level_number_cache.n_global_dofs,
4051  1,
4054  triangulation->get_communicator()) -
4055  1,
4056  triangulation->get_communicator());
4057  AssertThrowMPI(ierr);
4058 
4059  // shift indices
4060  for (types::global_dof_index &index : renumbering)
4061  if (index != numbers::invalid_dof_index)
4062  index += my_shift;
4063 
4064  // now re-enumerate all dofs to this shifted and condensed
4065  // numbering form. we renumber some dofs as invalid, so
4066  // choose the nocheck-version of the function
4067  //
4068  // of course there is nothing for us to renumber if the
4069  // level we are currently dealing with doesn't even exist
4070  // within the current triangulation, so skip renumbering
4071  // in that case
4072  if (level < triangulation->n_levels())
4074  renumbering, IndexSet(0), *dof_handler, level, false);
4075 
4076  // now a little bit of housekeeping
4077  level_number_cache.locally_owned_dofs =
4078  IndexSet(level_number_cache.n_global_dofs);
4079  level_number_cache.locally_owned_dofs.add_range(
4080  my_shift, my_shift + level_number_cache.n_locally_owned_dofs);
4081  level_number_cache.locally_owned_dofs.compress();
4082 
4083  number_caches.emplace_back(level_number_cache);
4084  }
4085 
4086 
4087  //* communicate ghost DoFs
4088  // We mark all ghost cells by setting the user_flag and then request
4089  // these cells from the corresponding owners. As this information
4090  // can be incomplete,
4091  {
4092  std::vector<bool> user_flags;
4093  triangulation->save_user_flags(user_flags);
4094  triangulation->clear_user_flags();
4095 
4096  // mark all ghost cells for transfer
4097  {
4098  for (const auto &cell : dof_handler->cell_iterators())
4099  if (cell->level_subdomain_id() !=
4101  !cell->is_locally_owned_on_level())
4102  cell->set_user_flag();
4103  }
4104 
4105  // Phase 1. Request all marked cells from corresponding owners. If we
4106  // managed to get every DoF, remove the user_flag, otherwise we
4107  // will request them again in the step below.
4108  communicate_mg_ghost_cells(*triangulation, *dof_handler);
4109 
4110  // have a barrier so that sends from above and below this
4111  // place are not mixed up.
4112  //
4113  // this is necessary because above we just see if there are
4114  // messages and then receive them, without discriminating
4115  // where they come from and whether they were sent in phase
4116  // 1 or 2 in communicate_mg_ghost_cells() on another
4117  // processor. the need for a global communication step like
4118  // this barrier could be avoided by receiving messages
4119  // specifically from those processors from which we expect
4120  // messages, and by using different tags for phase 1 and 2,
4121  // but the cost of a barrier is negligible compared to
4122  // everything else we do here
4123  const int ierr = MPI_Barrier(triangulation->get_communicator());
4124  AssertThrowMPI(ierr);
4125 
4126  // Phase 2, only request the cells that were not completed
4127  // in Phase 1.
4128  communicate_mg_ghost_cells(*triangulation, *dof_handler);
4129 
4130 # ifdef DEBUG
4131  // make sure we have removed all flags:
4132  {
4133  for (const auto &cell : dof_handler->cell_iterators())
4134  if (cell->level_subdomain_id() !=
4136  !cell->is_locally_owned_on_level())
4137  Assert(cell->user_flag_set() == false, ExcInternalError());
4138  }
4139 # endif
4140 
4141  triangulation->load_user_flags(user_flags);
4142  }
4143 
4144 
4145 
4146 # ifdef DEBUG
4147  // check that we are really done
4148  {
4149  std::vector<::types::global_dof_index> local_dof_indices;
4150  for (const auto &cell : dof_handler->cell_iterators())
4151  if (cell->level_subdomain_id() !=
4153  {
4154  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
4155  cell->get_mg_dof_indices(local_dof_indices);
4156  if (local_dof_indices.end() !=
4157  std::find(local_dof_indices.begin(),
4158  local_dof_indices.end(),
4160  {
4161  Assert(false, ExcMessage("not all DoFs got distributed!"));
4162  }
4163  }
4164  }
4165 # endif // DEBUG
4166 
4167  return number_caches;
4168 
4169 #endif // DEAL_II_WITH_MPI
4170  }
4171 
4172 
4173  template <int dim, int spacedim>
4174  NumberCache
4176  const std::vector<::types::global_dof_index> &new_numbers) const
4177  {
4178  (void)new_numbers;
4179 
4180  Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
4181  ExcInternalError());
4182 
4183 #ifndef DEAL_II_WITH_MPI
4184  Assert(false, ExcNotImplemented());
4185  return NumberCache();
4186 #else
4187 
4189  *triangulation =
4190  (dynamic_cast<
4192  const_cast<::Triangulation<dim, spacedim> *>(
4193  &dof_handler->get_triangulation())));
4194  Assert(triangulation != nullptr, ExcInternalError());
4195 
4196 
4197  // We start by checking whether only the numbering within the MPI
4198  // ranks changed. In that case, we can apply the renumbering with some
4199  // local renumbering only (this is similar to the renumber_mg_dofs()
4200  // function below)
4201  const bool locally_owned_set_changes =
4202  std::any_of(new_numbers.cbegin(),
4203  new_numbers.cend(),
4204  [this](const types::global_dof_index i) {
4205  return dof_handler->locally_owned_dofs().is_element(
4206  i) == false;
4207  });
4208 
4209  if (Utilities::MPI::sum(static_cast<unsigned int>(
4210  locally_owned_set_changes),
4211  triangulation->get_communicator()) == 0)
4212  {
4213  // Since only the order within the local subdomains has changed,
4214  // all we need to do is to propagate the knowledge about the
4215  // numbers from the locally owned dofs (given by the new_numbers
4216  // array) to all ghosted dofs on neighboring processors. We can do
4217  // this by ghost layer exchange routines as in parallel vectors:
4218  // We create an IndexSet for the relevant dofs and then export
4219  // into an array of those values via Utilities::MPI::Partitioner.
4220  IndexSet relevant_dofs;
4222  relevant_dofs);
4223  std::vector<types::global_dof_index> ghosted_new_numbers(
4224  relevant_dofs.n_elements());
4225  {
4226  Utilities::MPI::Partitioner partitioner(
4227  dof_handler->locally_owned_dofs(),
4228  relevant_dofs,
4229  triangulation->get_communicator());
4230 
4231  // choose some number that makes it unlikely to get conflicts
4232  // with other ongoing non-blocking communication (there
4233  // shouldn't be any at this place in most programs).
4234  const unsigned int communication_channel = 19;
4235  std::vector<types::global_dof_index> temp_array(
4236  partitioner.n_import_indices());
4237  std::vector<MPI_Request> requests;
4238  partitioner.export_to_ghosted_array_start(
4239  communication_channel,
4240  make_array_view(new_numbers),
4241  make_array_view(temp_array),
4243  ghosted_new_numbers.data() + new_numbers.size(),
4244  partitioner.n_ghost_indices()),
4245  requests);
4246  partitioner.export_to_ghosted_array_finish(
4248  ghosted_new_numbers.data() + new_numbers.size(),
4249  partitioner.n_ghost_indices()),
4250  requests);
4251 
4252  // we need to fill the indices of the locally owned part into
4253  // the new numbers array, which is not provided by the parallel
4254  // partitioner. their right position is somewhere in the middle
4255  // of the array, so we first copy the ghosted part from smaller
4256  // ranks to the front, then insert the data in the middle.
4257  unsigned int n_ghosts_on_smaller_ranks = 0;
4258  for (std::pair<unsigned int, unsigned int> t :
4259  partitioner.ghost_targets())
4260  {
4261  if (t.first > partitioner.this_mpi_process())
4262  break;
4263  n_ghosts_on_smaller_ranks += t.second;
4264  }
4265  if (n_ghosts_on_smaller_ranks > 0)
4266  {
4267  Assert(ghosted_new_numbers.data() != nullptr,
4268  ExcInternalError());
4269  std::memmove(ghosted_new_numbers.data(),
4270  ghosted_new_numbers.data() + new_numbers.size(),
4271  sizeof(types::global_dof_index) *
4272  n_ghosts_on_smaller_ranks);
4273  }
4274  if (new_numbers.size() > 0)
4275  {
4276  Assert(new_numbers.data() != nullptr, ExcInternalError());
4277  std::memcpy(ghosted_new_numbers.data() +
4278  n_ghosts_on_smaller_ranks,
4279  new_numbers.data(),
4280  sizeof(types::global_dof_index) *
4281  new_numbers.size());
4282  }
4283  }
4284 
4285  // In case we do not carry any relevant dof (but only some remote
4286  // processor), we do not need to call the renumbering. We call the
4287  // version without validity check because vertex dofs will be
4288  // set already in the artificial region.
4289  if (relevant_dofs.n_elements() > 0)
4290  Implementation::renumber_dofs(ghosted_new_numbers,
4291  relevant_dofs,
4292  *dof_handler,
4293  /*check_validity=*/false);
4294 
4295  NumberCache number_cache;
4296  number_cache.locally_owned_dofs = dof_handler->locally_owned_dofs();
4297  number_cache.n_global_dofs = dof_handler->n_dofs();
4298  number_cache.n_locally_owned_dofs =
4299  number_cache.locally_owned_dofs.n_elements();
4300  return number_cache;
4301  }
4302  else
4303  {
4304  // Now back to the more complicated case
4305  //
4306  // First figure out the new set of locally owned DoF indices.
4307  // If we own no DoFs, we still need to go through this function,
4308  // but we can skip this calculation.
4309  //
4310  // The IndexSet::add_indices() function is substantially more
4311  // efficient if the set of indices is already sorted because
4312  // it can then insert ranges instead of individual elements.
4313  // consequently, pre-sort the array of new indices
4314  IndexSet my_locally_owned_new_dof_indices(dof_handler->n_dofs());
4315  if (dof_handler->n_locally_owned_dofs() > 0)
4316  {
4317  std::vector<::types::global_dof_index>
4318  new_numbers_sorted = new_numbers;
4319  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4320 
4321  my_locally_owned_new_dof_indices.add_indices(
4322  new_numbers_sorted.begin(), new_numbers_sorted.end());
4323  my_locally_owned_new_dof_indices.compress();
4324 
4325  Assert(my_locally_owned_new_dof_indices.n_elements() ==
4326  new_numbers.size(),
4327  ExcInternalError());
4328  }
4329 
4330  // delete all knowledge of DoF indices that are not locally
4331  // owned. we do so by getting DoF indices on cells, checking
4332  // whether they are locally owned, if not, setting them to
4333  // an invalid value, and then setting them again on the current
4334  // cell
4335  //
4336  // DoFs we (i) know about, and (ii) don't own locally must be
4337  // located either on ghost cells, or on the interface between a
4338  // locally owned cell and a ghost cell. In any case, it is
4339  // sufficient to kill them only from the ghost side cell, so loop
4340  // only over ghost cells
4341  {
4342  std::vector<::types::global_dof_index> local_dof_indices;
4343 
4344  for (auto cell : dof_handler->active_cell_iterators())
4345  if (cell->is_ghost())
4346  {
4347  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
4348  cell->get_dof_indices(local_dof_indices);
4349 
4350  for (unsigned int i = 0;
4351  i < cell->get_fe().n_dofs_per_cell();
4352  ++i)
4353  // delete a DoF index if it has not already been deleted
4354  // (e.g., by visiting a neighboring cell, if it is on the
4355  // boundary), and if we don't own it
4356  if ((local_dof_indices[i] !=
4358  (!dof_handler->locally_owned_dofs().is_element(
4359  local_dof_indices[i])))
4360  local_dof_indices[i] = numbers::invalid_dof_index;
4361 
4362  cell->set_dof_indices(local_dof_indices);
4363  }
4364  }
4365 
4366 
4367  // renumber. Skip when there is nothing to do because we own no DoF.
4368  if (dof_handler->locally_owned_dofs().n_elements() > 0)
4369  Implementation::renumber_dofs(new_numbers,
4370  dof_handler->locally_owned_dofs(),
4371  *dof_handler,
4372  /*check_validity=*/false);
4373 
4374  // Communicate newly assigned DoF indices to other processors
4375  // and get the same information for our own ghost cells.
4376  //
4377  // This is the same as phase 5+6 in the distribute_dofs() algorithm,
4378  // taking into account that we have to unify a few DoFs in between
4379  // then communication phases if we do hp-numbering
4380  {
4381  std::vector<bool> user_flags;
4382  triangulation->save_user_flags(user_flags);
4383  triangulation->clear_user_flags();
4384 
4385  // mark all own cells for transfer
4386  for (const auto &cell : dof_handler->active_cell_iterators())
4387  if (cell->is_ghost())
4388  cell->set_user_flag();
4389 
4390 
4391  // Send and receive cells. After this, only the local cells
4392  // are marked, that received new data. This has to be
4393  // communicated in a second communication step.
4394  //
4395  // as explained in the 'distributed' paper, this has to be
4396  // done twice
4397  communicate_dof_indices_on_marked_cells(*dof_handler);
4398 
4399  // if the DoFHandler has hp-capabilities then we may have
4400  // received valid indices of degrees of freedom that are
4401  // dominated by a FE object adjacent to a ghost interface.
4402  // thus, we overwrite the remaining invalid indices with the
4403  // valid ones in this step.
4405  *dof_handler);
4406 
4407  communicate_dof_indices_on_marked_cells(*dof_handler);
4408 
4409  triangulation->load_user_flags(user_flags);
4410  }
4411 
4412  NumberCache number_cache;
4413  number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
4414  number_cache.n_global_dofs = dof_handler->n_dofs();
4415  number_cache.n_locally_owned_dofs =
4416  number_cache.locally_owned_dofs.n_elements();
4417  return number_cache;
4418  }
4419 #endif
4420  }
4421 
4422 
4423 
4424  template <int dim, int spacedim>
4425  NumberCache
4427  const unsigned int level,
4428  const std::vector<types::global_dof_index> &new_numbers) const
4429  {
4430  // we only implement the case where the multigrid numbers are
4431  // renumbered within the processor's partition, rather than the most
4432  // general case
4433  const IndexSet index_set = dof_handler->locally_owned_mg_dofs(level);
4434 
4435 #ifdef DEAL_II_WITH_MPI
4436 
4437  const ::parallel::TriangulationBase<dim, spacedim> *tr =
4438  (dynamic_cast<const ::parallel::TriangulationBase<dim, spacedim>
4439  *>(&this->dof_handler->get_triangulation()));
4440  Assert(tr != nullptr, ExcInternalError());
4441 
4442  const unsigned int my_rank =
4443  Utilities::MPI::this_mpi_process(tr->get_communicator());
4444 
4445 # ifdef DEBUG
4446  for (types::global_dof_index i : new_numbers)
4447  {
4448  Assert(index_set.is_element(i),
4450  "Renumberings that change the locally owned mg dofs "
4451  "partitioning are currently not implemented for "
4452  "the multigrid levels"));
4453  }
4454 # endif
4455 
4456  // we need to access all locally relevant degrees of freedom. we
4457  // use Utilities::MPI::Partitioner for handling the data exchange
4458  // of the new numbers, which is simply the extraction of ghost data
4459  IndexSet relevant_dofs;
4461  level,
4462  relevant_dofs);
4463  std::vector<types::global_dof_index> ghosted_new_numbers(
4464  relevant_dofs.n_elements());
4465  {
4466  Utilities::MPI::Partitioner partitioner(index_set,
4467  relevant_dofs,
4468  tr->get_communicator());
4469  std::vector<types::global_dof_index> temp_array(
4470  partitioner.n_import_indices());
4471  const unsigned int communication_channel = 17;
4472  std::vector<MPI_Request> requests;
4473  partitioner.export_to_ghosted_array_start(
4474  communication_channel,
4475  make_array_view(new_numbers),
4476  make_array_view(temp_array),
4477  ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
4478  new_numbers.size(),
4479  partitioner.n_ghost_indices()),
4480  requests);
4481  partitioner.export_to_ghosted_array_finish(
4482  ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
4483  new_numbers.size(),
4484  partitioner.n_ghost_indices()),
4485  requests);
4486 
4487  // we need to fill the indices of the locally owned part into the
4488  // new numbers array. their right position is somewhere in the
4489  // middle of the array, so we first copy the ghosted part from
4490  // smaller ranks to the front, then insert the data in the middle.
4491  unsigned int n_ghosts_on_smaller_ranks = 0;
4492  for (std::pair<unsigned int, unsigned int> t :
4493  partitioner.ghost_targets())
4494  {
4495  if (t.first > my_rank)
4496  break;
4497  n_ghosts_on_smaller_ranks += t.second;
4498  }
4499  if (n_ghosts_on_smaller_ranks > 0)
4500  {
4501  Assert(ghosted_new_numbers.data() != nullptr, ExcInternalError());
4502  std::memmove(ghosted_new_numbers.data(),
4503  ghosted_new_numbers.data() + new_numbers.size(),
4504  sizeof(types::global_dof_index) *
4505  n_ghosts_on_smaller_ranks);
4506  }
4507  if (new_numbers.size() > 0)
4508  {
4509  Assert(new_numbers.data() != nullptr, ExcInternalError());
4510  std::memcpy(ghosted_new_numbers.data() +
4511  n_ghosts_on_smaller_ranks,
4512  new_numbers.data(),
4513  sizeof(types::global_dof_index) * new_numbers.size());
4514  }
4515  }
4516 
4517  // in case we do not own any of the given level (but only some remote
4518  // processor), we do not need to call the renumbering
4519  if (level < this->dof_handler->get_triangulation().n_levels() &&
4520  relevant_dofs.n_elements() > 0)
4522  ghosted_new_numbers, relevant_dofs, *dof_handler, level, true);
4523 #else
4524  (void)new_numbers;
4525  Assert(false, ExcNotImplemented());
4526 #endif
4527 
4528  NumberCache number_cache;
4529  number_cache.locally_owned_dofs = index_set;
4530  number_cache.n_global_dofs = dof_handler->n_dofs(level);
4531  number_cache.n_locally_owned_dofs =
4532  number_cache.locally_owned_dofs.n_elements();
4533  return number_cache;
4534  }
4535  } // namespace Policy
4536  } // namespace DoFHandlerImplementation
4537 } // namespace internal
4538 
4539 
4540 
4541 /*-------------- Explicit Instantiations -------------------------------*/
4542 #include "dof_handler_policy.inst"
4543 
4544 
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const DoFHandler< dim, spacedim > &dof_handler)
ParallelShared(DoFHandler< dim, spacedim > &dof_handler)
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1639
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
static void renumber_vertex_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
static void renumber_face_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 2, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
const types::subdomain_id invalid_subdomain_id
Definition: types.h:276
void clear_user_flags()
Definition: tria.cc:11186
static types::global_dof_index unify_dof_indices(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int n_dofs_before_identification, const bool check_validity)
cell_iterator begin(const unsigned int level=0) const
static void renumber_vertex_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const bool check_validity)
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 3, spacedim > &dof_handler)
Task< RT > new_task(const std::function< RT()> &function)
static void compute_dof_identities(std::vector< std::map< types::global_dof_index, types::global_dof_index >> &all_constrained_indices, const DoFHandler< dim, spacedim > &dof_handler)
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
void load_user_flags(std::istream &in)
Definition: tria.cc:11246
cell_iterator end() const
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition: types.h:86
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
static types::global_dof_index enumerate_dof_indices_for_renumbering(std::vector< types::global_dof_index > &new_dof_indices, const std::vector< std::map< types::global_dof_index, types::global_dof_index >> &all_constrained_indices, const DoFHandler< dim, spacedim > &)
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
virtual std::vector< NumberCache > distribute_mg_dofs() const override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
Point< 2 > second
Definition: grid_out.cc:4588
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 2, spacedim > &dof_handler)
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition: dof_handler.h:1646
static types::global_dof_index distribute_dofs(const types::subdomain_id subdomain_id, DoFHandler< dim, spacedim > &dof_handler)
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1414
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
void extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1252
static void merge_invalid_vertex_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
static std::map< types::global_dof_index, types::global_dof_index > compute_vertex_dof_identities(const DoFHandler< dim, spacedim > &dof_handler)
static void renumber_face_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &)
size_type size() const
Definition: index_set.h:1634
static void merge_invalid_quad_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
active_cell_iterator begin_active(const unsigned int level=0) const
virtual MPI_Comm get_communicator() const override
Definition: tria_base.cc:140
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
virtual std::vector< NumberCache > distribute_mg_dofs() const override
static ::ExceptionBase & ExcMessage(std::string arg1)
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler)
unsigned int subdomain_id
Definition: types.h:43
static void merge_invalid_line_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
static void invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering(std::vector< types::global_dof_index > &renumbering, const types::subdomain_id subdomain_id, const DoFHandler< dim, spacedim > &dof_handler)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1473
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1550
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const DoFHandler< dim, spacedim > &dof_handler)
IteratorRange< active_cell_iterator > active_cell_iterators() const
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1921
types::global_dof_index n_dofs() const
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition: dof_handler.h:1595
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:697
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
unsigned int level
Definition: grid_out.cc:4590
static void renumber_cell_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler)
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:502
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
types::subdomain_id locally_owned_subdomain() const override
Definition: tria_base.cc:320
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:978
static void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, const DoFHandler< dim, space_dim > &dof_handler, const bool check_validity)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1218
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SmartPointer< DoFHandler< dim, spacedim > > dof_handler
Sequential(DoFHandler< dim, spacedim > &dof_handler)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
Point< 2 > first
Definition: grid_out.cc:4587
static void merge_invalid_line_dofs_on_ghost_interfaces(DoFHandler< 1, spacedim > &dof_handler)
#define DEAL_II_MPI_CONST_CAST(expr)
Definition: mpi.h:82
static types::global_dof_index distribute_dofs_on_level(const types::subdomain_id level_subdomain_id, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
const IndexSet & locally_owned_dofs() const
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1673
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const DoFHandler< 3, spacedim > &dof_handler)
unsigned int global_dof_index
Definition: types.h:76
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
void compress() const
Definition: index_set.h:1642
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:989
static void renumber_face_mg_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &, const unsigned int, const bool)
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1754
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
constexpr const ReferenceCell Quadrilateral
SmartPointer< DoFHandler< dim, spacedim > > dof_handler
static void merge_invalid_quad_dofs_on_ghost_interfaces(DoFHandler< 3, spacedim > &dof_handler)
const std::vector< IndexSet > & locally_owned_dofs_per_processor() const
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1399
void export_to_ghosted_array_start(const unsigned int communication_channel, const ArrayView< const Number, MemorySpaceType > &locally_owned_array, const ArrayView< Number, MemorySpaceType > &temporary_storage, const ArrayView< Number, MemorySpaceType > &ghost_array, std::vector< MPI_Request > &requests) const
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1321
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
void save_user_flags(std::ostream &out) const
Definition: tria.cc:11197
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcNotImplemented()
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1337
static unsigned int n_threads()
virtual std::vector< NumberCache > distribute_mg_dofs() const override
bool is_element(const size_type index) const
Definition: index_set.h:1765
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other, const unsigned int face_no=0) const
Definition: fe.cc:1000
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
static void set_dof_index(const DoFHandler< 1, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 1 >> &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 1 >> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const std::integral_constant< int, 1 >)
Definition: dof_handler.cc:801
void extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1210
static void renumber_cell_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
static void renumber_face_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 3, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
static void renumber_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
const types::global_dof_index invalid_dof_index
Definition: types.h:211
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const DoFHandler< 1, spacedim > &dof_handler)
size_type n_elements() const
Definition: index_set.h:1832
IteratorRange< cell_iterator > cell_iterators() const
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1880
bool hp_capability_enabled
Definition: dof_handler.h:1537
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:1005
types::global_dof_index n_locally_owned_dofs() const
static void merge_invalid_dof_indices_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()