Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_handler_policy.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2023 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/types.h>
22 #include <deal.II/base/utilities.h>
24 
27 
31 
32 #include <deal.II/fe/fe.h>
33 
35 #include <deal.II/grid/tria.h>
37 
38 #include <algorithm>
39 #include <limits>
40 #include <memory>
41 #include <numeric>
42 #include <set>
43 
45 
46 
47 namespace internal
48 {
49  namespace DoFHandlerImplementation
50  {
51  namespace Policy
52  {
53  namespace
54  {
61  const types::global_dof_index enumeration_dof_index =
63 
64 
65  using DoFIdentities =
66  std::vector<std::pair<unsigned int, unsigned int>>;
67 
68 
79  template <int structdim, int dim, int spacedim>
80  const std::unique_ptr<DoFIdentities> &
81  ensure_existence_and_return_dof_identities(
82  const ::hp::FECollection<dim, spacedim> &fes,
83  const types::fe_index fe_index_1,
84  const types::fe_index fe_index_2,
85  std::unique_ptr<DoFIdentities> &identities,
86  const unsigned int face_no = numbers::invalid_unsigned_int)
87  {
88  Assert(structdim == 2 || face_no == numbers::invalid_unsigned_int,
90 
91  // see if we need to fill this entry, or whether it already
92  // exists
93  if (identities.get() == nullptr)
94  {
95  // TODO: Change to
96  // std::vector<std::map<types::fe_index, unsigned int>>
97  std::vector<std::map<unsigned int, unsigned int>>
98  complete_identities;
99 
100  switch (structdim)
101  {
102  case 0:
103  {
104  // TODO: Change set to types::fe_index
105  complete_identities = fes.hp_vertex_dof_identities(
106  std::set<unsigned int>{fe_index_1, fe_index_2});
107  break;
108  }
109 
110  case 1:
111  {
112  // TODO: Change set to types::fe_index
113  complete_identities = fes.hp_line_dof_identities(
114  std::set<unsigned int>{fe_index_1, fe_index_2});
115  break;
116  }
117 
118  case 2:
119  {
120  // TODO: Change set to types::fe_index
121  complete_identities = fes.hp_quad_dof_identities(
122  std::set<unsigned int>{fe_index_1, fe_index_2},
123  face_no);
124  break;
125  }
126 
127  default:
128  Assert(false, ExcNotImplemented());
129  }
130 
131 #ifdef DEBUG
132  // Each entry of 'complete_identities' contains a set of
133  // pairs (fe_index,dof_index). Because we put in exactly
134  // two fe indices, we know that each entry of the outer
135  // vector needs to contain a set of exactly two such
136  // pairs. Check this. While there, also check that
137  // the two entries actually reference fe_index_1 and
138  // fe_index_2:
139  for (const auto &complete_identity : complete_identities)
140  {
141  Assert(complete_identity.size() == 2, ExcInternalError());
142  Assert(complete_identity.find(fe_index_1) !=
143  complete_identity.end(),
144  ExcInternalError());
145  Assert(complete_identity.find(fe_index_2) !=
146  complete_identity.end(),
147  ExcInternalError());
148  }
149 #endif
150 
151  // Next reduce these sets of two pairs by removing the
152  // fe_index parts: We know which indices we have. But we
153  // have to make sure in which order we consider the
154  // pair, by considering whether the fe_index part we are
155  // throwing away matched fe_index_1 or fe_index_2. Fortunately,
156  // this is easy to do because we can ask the std::map for the
157  // dof_index that matches a given fe_index:
158  DoFIdentities reduced_identities;
159  for (const auto &complete_identity : complete_identities)
160  {
161  const unsigned int dof_index_1 =
162  complete_identity.at(fe_index_1);
163  const unsigned int dof_index_2 =
164  complete_identity.at(fe_index_2);
165 
166  reduced_identities.emplace_back(dof_index_1, dof_index_2);
167  }
168 
169 #ifdef DEBUG
170  // double check whether the newly created entries make
171  // any sense at all
172  for (const auto &identity : reduced_identities)
173  {
174  Assert(identity.first <
175  fes[fe_index_1]
176  .template n_dofs_per_object<structdim>(face_no),
177  ExcInternalError());
178  Assert(identity.second <
179  fes[fe_index_2]
180  .template n_dofs_per_object<structdim>(face_no),
181  ExcInternalError());
182  }
183 #endif
184 
185  identities =
186  std::make_unique<DoFIdentities>(std::move(reduced_identities));
187  }
188 
189  return identities;
190  }
191  } // namespace
192 
193 
194 
196  {
197  /* -------------- distribute_dofs functionality ------------- */
198 
203  template <int dim, int spacedim>
204  static std::map<types::global_dof_index, types::global_dof_index>
206  const DoFHandler<dim, spacedim> &dof_handler)
207  {
208  Assert(
209  dof_handler.hp_capability_enabled == true,
211 
212  std::map<types::global_dof_index, types::global_dof_index>
213  dof_identities;
214 
215  // Note: we may wish to have something here similar to what
216  // we do for lines and quads, namely that we only identify
217  // dofs for any FE towards the most dominating one. however,
218  // it is not clear whether this is actually necessary for
219  // vertices at all, I can't think of a finite element that
220  // would make that necessary...
222  vertex_dof_identities(dof_handler.get_fe_collection().size(),
223  dof_handler.get_fe_collection().size());
224 
225  // loop over all vertices and see which one we need to work on
226  for (unsigned int vertex_index = 0;
227  vertex_index < dof_handler.get_triangulation().n_vertices();
228  ++vertex_index)
229  if (dof_handler.get_triangulation()
230  .get_used_vertices()[vertex_index] == true)
231  {
232  const unsigned int n_active_fe_indices =
233  ::internal::DoFAccessorImplementation::Implementation::
234  n_active_fe_indices(dof_handler,
235  0,
236  vertex_index,
237  std::integral_constant<int, 0>());
238 
239  if (n_active_fe_indices > 1)
240  {
241  const std::set<types::fe_index> fe_indices =
244  dof_handler,
245  0,
246  vertex_index,
247  std::integral_constant<int, 0>());
248 
249  // find out which is the most dominating finite
250  // element of the ones that are used on this vertex
251  // TODO: Change set to types::fe_index
252  types::fe_index most_dominating_fe_index =
254  {fe_indices.begin(), fe_indices.end()},
255  /*codim*/ dim);
256 
257  // if we haven't found a dominating finite element,
258  // choose the very first one to be dominant
259  // TODO: Change assert to numbers::invalid_fe_index
260  if (most_dominating_fe_index == numbers::invalid_fe_index)
261  most_dominating_fe_index =
262  ::internal::DoFAccessorImplementation::
263  Implementation::nth_active_fe_index(
264  dof_handler,
265  0,
266  vertex_index,
267  0,
268  std::integral_constant<int, 0>());
269 
270  // loop over the indices of all the finite
271  // elements that are not dominating, and
272  // identify their dofs to the most dominating
273  // one
274  for (const auto &other_fe_index : fe_indices)
275  if (other_fe_index != most_dominating_fe_index)
276  {
277  // make sure the entry in the equivalence
278  // table exists
279  const auto &identities =
280  *ensure_existence_and_return_dof_identities<0>(
281  dof_handler.get_fe_collection(),
282  most_dominating_fe_index,
283  other_fe_index,
284  vertex_dof_identities[most_dominating_fe_index]
285  [other_fe_index]);
286 
287  // then loop through the identities we
288  // have. first get the global numbers of the
289  // dofs we want to identify and make sure they
290  // are not yet constrained to anything else,
291  // except for to each other. use the rule that
292  // we will always constrain the dof with the
293  // higher FE index to the one with the lower,
294  // to avoid circular reasoning.
295  for (const auto &identity : identities)
296  {
297  const types::global_dof_index primary_dof_index =
298  ::internal::DoFAccessorImplementation::
299  Implementation::get_dof_index(
300  dof_handler,
301  0,
302  vertex_index,
303  most_dominating_fe_index,
304  identity.first,
305  std::integral_constant<int, 0>());
307  dependent_dof_index =
308  ::internal::DoFAccessorImplementation::
309  Implementation::get_dof_index(
310  dof_handler,
311  0,
312  vertex_index,
313  other_fe_index,
314  identity.second,
315  std::integral_constant<int, 0>());
316 
317  // on subdomain boundaries, we will
318  // encounter invalid DoFs on ghost cells,
319  // for which we have not yet distributed
320  // valid indices. depending on which finte
321  // element is dominating the other on this
322  // interface, we either have to constrain
323  // the valid to the invalid indices, or vice
324  // versa.
325  //
326  // we only store an identity if we are about
327  // to overwrite a valid DoF. we will skip
328  // constraining invalid DoFs for now, and
329  // consider them later in Phase 5.
330  if (dependent_dof_index !=
332  {
333  // if the DoF indices of both elements
334  // are already distributed, i.e., both
335  // of these 'fe_indices' are associated
336  // with a locally owned cell, then we
337  // should either not have a dof_identity
338  // yet, or it must come out here to be
339  // exactly as we had computed before
340  if (primary_dof_index !=
342  Assert(
343  (dof_identities.find(primary_dof_index) ==
344  dof_identities.end()) ||
345  (dof_identities[dependent_dof_index] ==
346  primary_dof_index),
347  ExcInternalError());
348 
349  dof_identities[dependent_dof_index] =
350  primary_dof_index;
351  }
352  }
353  }
354  }
355  }
356 
357  return dof_identities;
358  }
359 
360 
365  template <int spacedim>
366  static std::map<types::global_dof_index, types::global_dof_index>
368  {
369  (void)dof_handler;
370  Assert(dof_handler.hp_capability_enabled == true,
372 
373  return std::map<types::global_dof_index, types::global_dof_index>();
374  }
375 
376 
377  template <int dim, int spacedim>
378  static std::map<types::global_dof_index, types::global_dof_index>
380  const DoFHandler<dim, spacedim> &dof_handler)
381  {
382  Assert(
383  dof_handler.hp_capability_enabled == true,
385 
386  std::map<types::global_dof_index, types::global_dof_index>
387  dof_identities;
388 
389  // An implementation of the algorithm described in the hp-paper,
390  // including the modification mentioned later in the "complications in
391  // 3-d" subsections
392  //
393  // as explained there, we do something only if there are exactly 2
394  // finite elements associated with an object. if there is only one,
395  // then there is nothing to do anyway, and if there are 3 or more,
396  // then we can get into trouble. note that this only happens for lines
397  // in 3d and higher, and for quads only in 4d and higher, so this
398  // isn't a particularly frequent case
399  //
400  // there is one case, however, that we would like to handle (see, for
401  // example, the hp/crash_15 testcase): if we have
402  // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
403  // then we should be able to handle this because we can simply unify
404  // *all* dofs, not only a some. so what we do is to first treat all
405  // pairs of finite elements that have *identical* dofs, and then only
406  // deal with those that are not identical of which we can handle at
407  // most 2
408  ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
409  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
410 
411  std::vector<bool> line_touched(
412  dof_handler.get_triangulation().n_raw_lines());
413  for (const auto &cell : dof_handler.active_cell_iterators())
414  for (const auto l : cell->line_indices())
415  if (!line_touched[cell->line(l)->index()])
416  {
417  const auto line = cell->line(l);
418  line_touched[line->index()] = true;
419 
420  unsigned int unique_sets_of_dofs =
421  line->n_active_fe_indices();
422 
423  // do a first loop over all sets of dofs and do identity
424  // uniquification
425  const unsigned int n_active_fe_indices =
426  line->n_active_fe_indices();
427  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
428  for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
429  {
430  const types::fe_index fe_index_1 =
431  line->nth_active_fe_index(f),
432  fe_index_2 =
433  line->nth_active_fe_index(g);
434 
435  // as described in the hp-paper, we only unify on lines
436  // when there are at most two different FE objects
437  // assigned on it.
438  // however, more than two 'active_fe_indices' can be
439  // attached that still fulfill the above criterion,
440  // i.e. when two different FiniteElement objects are
441  // assigned to neighboring cells that map their degrees
442  // of freedom one-to-one.
443  // we cannot verify with certainty if two dofs each of
444  // separate FiniteElement objects actually map
445  // one-to-one. however, checking for the number of
446  // 'dofs_per_line' turned out to be a reasonable
447  // approach, that also works for e.g. two different
448  // FE_Q objects of the same order, from which one is
449  // enhanced by a bubble function that is zero on the
450  // boundary.
451  if ((dof_handler.get_fe(fe_index_1).n_dofs_per_line() ==
452  dof_handler.get_fe(fe_index_2)
453  .n_dofs_per_line()) &&
454  (dof_handler.get_fe(fe_index_1).n_dofs_per_line() >
455  0))
456  {
457  // the number of dofs per line is identical
458  const unsigned int dofs_per_line =
459  dof_handler.get_fe(fe_index_1).n_dofs_per_line();
460 
461  const auto &identities =
462  *ensure_existence_and_return_dof_identities<1>(
463  dof_handler.get_fe_collection(),
464  fe_index_1,
465  fe_index_2,
466  line_dof_identities[fe_index_1][fe_index_2]);
467  // see if these sets of dofs are identical. the
468  // first condition for this is that indeed there are
469  // n identities
470  if (identities.size() == dofs_per_line)
471  {
472  unsigned int i = 0;
473  for (; i < dofs_per_line; ++i)
474  if ((identities[i].first != i) &&
475  (identities[i].second != i))
476  // not an identity
477  break;
478 
479  if (i == dofs_per_line)
480  {
481  // The line dofs (i.e., the ones interior to
482  // a line) of these two finite elements are
483  // identical. Note that there could be
484  // situations when one element still
485  // dominates another, e.g.: FE_Q(2) x
486  // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
487 
488  --unique_sets_of_dofs;
489 
490  // determine which one of both finite
491  // elements is the dominating one.
492  const std::set<types::fe_index> fe_indices{
493  fe_index_1, fe_index_2};
494 
495  // TODO: Change set to types::fe_index
496  types::fe_index dominating_fe_index =
497  dof_handler.get_fe_collection()
498  .find_dominating_fe({fe_indices.begin(),
499  fe_indices.end()},
500  /*codim=*/dim - 1);
501  types::fe_index other_fe_index =
503 
504  if (dominating_fe_index !=
506  other_fe_index =
507  (dominating_fe_index == fe_index_1) ?
508  fe_index_2 :
509  fe_index_1;
510  else
511  {
512  // if we haven't found a dominating
513  // finite element, choose the one with
514  // the lower index to be dominating
515  dominating_fe_index = fe_index_1;
516  other_fe_index = fe_index_2;
517  }
518 
519  for (unsigned int j = 0; j < dofs_per_line;
520  ++j)
521  {
523  primary_dof_index = line->dof_index(
524  j, dominating_fe_index);
526  dependent_dof_index =
527  line->dof_index(j, other_fe_index);
528 
529  // on subdomain boundaries, we will
530  // encounter invalid DoFs on ghost
531  // cells, for which we have not yet
532  // distributed valid indices. depending
533  // on which finte element is dominating
534  // the other on this interface, we
535  // either have to constrain the valid to
536  // the invalid indices, or vice versa.
537  //
538  // we only store an identity if we are
539  // about to overwrite a valid DoF. we
540  // will skip constraining invalid DoFs
541  // for now, and consider them later in
542  // Phase 5.
543  if (dependent_dof_index !=
545  {
546  if (primary_dof_index !=
548  {
549  // if primary dof was already
550  // constrained, constrain to
551  // that one, otherwise constrain
552  // dependent to primary
553  if (dof_identities.find(
554  primary_dof_index) !=
555  dof_identities.end())
556  {
557  // if the DoF indices of
558  // both elements are already
559  // distributed, i.e., both
560  // of these 'fe_indices' are
561  // associated with a locally
562  // owned cell, then we
563  // should either not have a
564  // dof_identity yet, or it
565  // must come out here to be
566  // exactly as we had
567  // computed before
568  Assert(
569  dof_identities.find(
570  dof_identities
571  [primary_dof_index]) ==
572  dof_identities.end(),
573  ExcInternalError());
574 
575  dof_identities
576  [dependent_dof_index] =
577  dof_identities
578  [primary_dof_index];
579  }
580  else
581  {
582  // see comment above for an
583  // explanation of this
584  // assertion
585  Assert(
586  (dof_identities.find(
587  primary_dof_index) ==
588  dof_identities.end()) ||
589  (dof_identities
590  [dependent_dof_index] ==
591  primary_dof_index),
592  ExcInternalError());
593 
594  dof_identities
595  [dependent_dof_index] =
596  primary_dof_index;
597  }
598  }
599  else
600  {
601  // set dependent_dof to
602  // primary_dof_index, which is
603  // invalid
604  dof_identities
605  [dependent_dof_index] =
607  }
608  }
609  }
610  }
611  }
612  }
613  }
614 
615  // if at this point, there is only one unique set of dofs
616  // left, then we have taken care of everything above. if there
617  // are two, then we need to deal with them here. if there are
618  // more, then we punt, as described in the paper (and
619  // mentioned above)
620  // TODO: The check for 'dim==2' was inserted by intuition. It
621  // fixes
622  // the previous problems with @ref step_27 "step-27" in 3d. But an
623  // explanation for this is still required, and what we do here
624  // is not what we describe in the paper!.
625  if ((unique_sets_of_dofs == 2) && (dim == 2))
626  {
627  const std::set<types::fe_index> fe_indices =
628  line->get_active_fe_indices();
629 
630  // find out which is the most dominating finite element of
631  // the ones that are used on this line
632  // TODO: Change set to types::fe_index
633  const types::fe_index most_dominating_fe_index =
635  {fe_indices.begin(), fe_indices.end()},
636  /*codim=*/dim - 1);
637 
638  // if we found the most dominating element, then use this
639  // to eliminate some of the degrees of freedom by
640  // identification. otherwise, the code that computes
641  // hanging node constraints will have to deal with it by
642  // computing appropriate constraints along this face/edge
643  if (most_dominating_fe_index != numbers::invalid_fe_index)
644  {
645  // loop over the indices of all the finite elements
646  // that are not dominating, and identify their dofs to
647  // the most dominating one
648  for (const auto &other_fe_index : fe_indices)
649  if (other_fe_index != most_dominating_fe_index)
650  {
651  const auto &identities =
652  *ensure_existence_and_return_dof_identities<
653  1>(dof_handler.get_fe_collection(),
654  most_dominating_fe_index,
655  other_fe_index,
656  line_dof_identities
657  [most_dominating_fe_index]
658  [other_fe_index]);
659 
660  for (const auto &identity : identities)
661  {
663  primary_dof_index = line->dof_index(
664  identity.first,
665  most_dominating_fe_index);
667  dependent_dof_index =
668  line->dof_index(identity.second,
669  other_fe_index);
670 
671  // on subdomain boundaries, we will
672  // encounter invalid DoFs on ghost cells,
673  // for which we have not yet distributed
674  // valid indices. depending on which finte
675  // element is dominating the other on this
676  // interface, we either have to constrain
677  // the valid to the invalid indices, or vice
678  // versa.
679  //
680  // we only store an identity if we are about
681  // to overwrite a valid DoF. we will skip
682  // constraining invalid DoFs for now, and
683  // consider them later in Phase 5.
684  if (dependent_dof_index !=
686  {
687  // if the DoF indices of both elements
688  // are already distributed, i.e., both
689  // of these 'fe_indices' are associated
690  // with a locally owned cell, then we
691  // should either not have a dof_identity
692  // yet, or it must come out here to be
693  // exactly as we had computed before
694  if (primary_dof_index !=
696  Assert((dof_identities.find(
697  primary_dof_index) ==
698  dof_identities.end()) ||
699  (dof_identities
700  [dependent_dof_index] ==
701  primary_dof_index),
702  ExcInternalError());
703 
704  dof_identities[dependent_dof_index] =
705  primary_dof_index;
706  }
707  }
708  }
709  }
710  }
711  }
712 
713  return dof_identities;
714  }
715 
716 
717 
722  template <int dim, int spacedim>
723  static std::map<types::global_dof_index, types::global_dof_index>
725  const DoFHandler<dim, spacedim> &dof_handler)
726  {
727  (void)dof_handler;
728  Assert(
729  dof_handler.hp_capability_enabled == true,
731 
732  // this function should only be called for dim<3 where there are
733  // no quad dof identities. for dim==3, the specialization below should
734  // take care of it
735  Assert(dim < 3, ExcInternalError());
736 
737  return std::map<types::global_dof_index, types::global_dof_index>();
738  }
739 
740 
741  template <int spacedim>
742  static std::map<types::global_dof_index, types::global_dof_index>
744  {
745  Assert(dof_handler.hp_capability_enabled == true,
747 
748  const int dim = 3;
749 
750  std::map<types::global_dof_index, types::global_dof_index>
751  dof_identities;
752 
753  // An implementation of the algorithm described in the hp-
754  // paper, including the modification mentioned later in the
755  // "complications in 3-d" subsections
756  //
757  // as explained there, we do something only if there are
758  // exactly 2 finite elements associated with an object. if
759  // there is only one, then there is nothing to do anyway,
760  // and if there are 3 or more, then we can get into
761  // trouble. note that this only happens for lines in 3d and
762  // higher, and for quads only in 4d and higher, so this
763  // isn't a particularly frequent case
764  ::Table<3, std::unique_ptr<DoFIdentities>> quad_dof_identities(
765  dof_handler.fe_collection.size(),
766  dof_handler.fe_collection.size(),
767  2 /*triangle (0) or quadrilateral (1)*/);
768 
769  std::vector<bool> quad_touched(
770  dof_handler.get_triangulation().n_raw_quads());
771  for (const auto &cell : dof_handler.active_cell_iterators())
772  for (const auto q : cell->face_indices())
773  if (!quad_touched[cell->quad(q)->index()] &&
774  (cell->quad(q)->n_active_fe_indices() == 2))
775  {
776  const auto quad = cell->quad(q);
777  quad_touched[quad->index()] = true;
778 
779  const std::set<types::fe_index> fe_indices =
780  quad->get_active_fe_indices();
781 
782  // find out which is the most dominating finite
783  // element of the ones that are used on this quad
784  // TODO: Change set to types::fe_index
785  const types::fe_index most_dominating_fe_index =
787  {fe_indices.begin(), fe_indices.end()},
788  /*codim=*/dim - 2);
789 
790  const unsigned int most_dominating_fe_index_face_no =
791  cell->active_fe_index() == most_dominating_fe_index ?
792  q :
793  cell->neighbor_face_no(q);
794 
795  // if we found the most dominating element, then use
796  // this to eliminate some of the degrees of freedom
797  // by identification. otherwise, the code that
798  // computes hanging node constraints will have to
799  // deal with it by computing appropriate constraints
800  // along this face/edge
801  if (most_dominating_fe_index != numbers::invalid_fe_index)
802  {
803  // loop over the indices of all the finite
804  // elements that are not dominating, and
805  // identify their dofs to the most dominating
806  // one
807  for (const auto &other_fe_index : fe_indices)
808  if (other_fe_index != most_dominating_fe_index)
809  {
810  const auto &identities =
811  *ensure_existence_and_return_dof_identities<2>(
812  dof_handler.get_fe_collection(),
813  most_dominating_fe_index,
814  other_fe_index,
815  quad_dof_identities
816  [most_dominating_fe_index][other_fe_index]
817  [cell->quad(q)->reference_cell() ==
819  most_dominating_fe_index_face_no);
820 
821  for (const auto &identity : identities)
822  {
824  primary_dof_index =
825  quad->dof_index(identity.first,
826  most_dominating_fe_index);
828  dependent_dof_index =
829  quad->dof_index(identity.second,
830  other_fe_index);
831 
832  // we only store an identity if we are about to
833  // overwrite a valid degree of freedom. we will
834  // skip invalid degrees of freedom (that are
835  // associated with ghost cells) for now, and
836  // consider them later in phase 5.
837  if (dependent_dof_index !=
839  {
840  // if the DoF indices of both elements are
841  // already distributed, i.e., both of these
842  // 'fe_indices' are associated with a
843  // locally owned cell, then we should either
844  // not have a dof_identity yet, or it must
845  // come out here to be exactly as we had
846  // computed before
847  if (primary_dof_index !=
849  Assert((dof_identities.find(
850  primary_dof_index) ==
851  dof_identities.end()) ||
852  (dof_identities
853  [dependent_dof_index] ==
854  primary_dof_index),
855  ExcInternalError());
856 
857  dof_identities[dependent_dof_index] =
858  primary_dof_index;
859  }
860  }
861  }
862  }
863  }
864 
865  return dof_identities;
866  }
867 
868 
869 
874  template <int dim, int spacedim>
875  static void
878  &all_constrained_indices,
879  const DoFHandler<dim, spacedim> &dof_handler)
880  {
881  if (dof_handler.hp_capability_enabled == false)
882  return;
883 
884  Assert(all_constrained_indices.size() == dim, ExcInternalError());
885 
886  Threads::TaskGroup<> tasks;
887 
888  unsigned int i = 0;
889  tasks += Threads::new_task([&, i]() {
890  all_constrained_indices[i] =
891  compute_vertex_dof_identities(dof_handler);
892  });
893 
894  if (dim > 1)
895  {
896  ++i;
897  tasks += Threads::new_task([&, i]() {
898  all_constrained_indices[i] =
899  compute_line_dof_identities(dof_handler);
900  });
901  }
902 
903  if (dim > 2)
904  {
905  ++i;
906  tasks += Threads::new_task([&, i]() {
907  all_constrained_indices[i] =
908  compute_quad_dof_identities(dof_handler);
909  });
910  }
911 
912  tasks.join_all();
913  }
914 
915 
916 
938  std::vector<types::global_dof_index> &new_dof_indices,
939  const std::vector<
940  std::map<types::global_dof_index, types::global_dof_index>>
941  &all_constrained_indices,
942  const types::global_dof_index start_dof_index)
943  {
944  // first preset the new DoF indices that are identities
945  for (const auto &constrained_dof_indices : all_constrained_indices)
946  for (const auto &p : constrained_dof_indices)
947  if (new_dof_indices[p.first] != numbers::invalid_dof_index)
948  {
949  Assert(new_dof_indices[p.first] == enumeration_dof_index,
950  ExcInternalError());
951 
952  new_dof_indices[p.first] = p.second;
953  }
954 
955  // then enumerate the rest
956  types::global_dof_index next_free_dof = start_dof_index;
957  for (auto &new_dof_index : new_dof_indices)
958  if (new_dof_index == enumeration_dof_index)
959  new_dof_index = next_free_dof++;
960 
961  // then loop over all those that are constrained and record the
962  // new dof number for those
963  for (const auto &constrained_dof_indices : all_constrained_indices)
964  for (const auto &p : constrained_dof_indices)
965  if (new_dof_indices[p.first] != numbers::invalid_dof_index)
966  {
967  Assert(new_dof_indices[p.first] != enumeration_dof_index,
968  ExcInternalError());
969 
970  if (p.second != numbers::invalid_dof_index)
971  new_dof_indices[p.first] = new_dof_indices[p.second];
972  }
973 
974  for (const types::global_dof_index new_dof_index : new_dof_indices)
975  {
976  (void)new_dof_index;
977  Assert(new_dof_index != enumeration_dof_index,
978  ExcInternalError());
979  Assert(new_dof_index < next_free_dof ||
980  new_dof_index == numbers::invalid_dof_index,
981  ExcInternalError());
982  }
983 
984  return next_free_dof;
985  }
986 
987 
988 
997  template <int dim, int spacedim>
1000  const unsigned int n_dofs_before_identification,
1001  const bool check_validity)
1002  {
1003  if (dof_handler.hp_capability_enabled == false)
1004  return n_dofs_before_identification;
1005 
1006  std::vector<
1007  std::map<types::global_dof_index, types::global_dof_index>>
1008  all_constrained_indices(dim);
1009  compute_dof_identities(all_constrained_indices, dof_handler);
1010 
1011  std::vector<::types::global_dof_index> renumbering(
1012  n_dofs_before_identification, enumeration_dof_index);
1013  const types::global_dof_index n_dofs =
1015  all_constrained_indices,
1016  0);
1017 
1018  renumber_dofs(renumbering, IndexSet(0), dof_handler, check_validity);
1019 
1020  return n_dofs;
1021  }
1022 
1023 
1024 
1029  template <int dim, int spacedim>
1030  static void
1032  DoFHandler<dim, spacedim> &dof_handler)
1033  {
1034  Assert(
1035  dof_handler.hp_capability_enabled == true,
1037 
1038  // Note: we may wish to have something here similar to what
1039  // we do for lines and quads, namely that we only identify
1040  // dofs for any FE towards the most dominating one. however,
1041  // it is not clear whether this is actually necessary for
1042  // vertices at all, I can't think of a finite element that
1043  // would make that necessary...
1045  vertex_dof_identities(dof_handler.get_fe_collection().size(),
1046  dof_handler.get_fe_collection().size());
1047 
1048  // mark all vertices on ghost cells to identify those cells that we
1049  // have already treated
1050  std::vector<bool> include_vertex(
1051  dof_handler.get_triangulation().n_vertices(), false);
1052  if (dynamic_cast<const ::parallel::
1053  DistributedTriangulationBase<dim, spacedim> *>(
1054  &dof_handler.get_triangulation()) != nullptr)
1055  for (const auto &cell : dof_handler.active_cell_iterators())
1056  if (cell->is_ghost())
1057  for (const unsigned int v : cell->vertex_indices())
1058  include_vertex[cell->vertex_index(v)] = true;
1059 
1060  // loop over all vertices and see which one we need to work on
1061  for (unsigned int vertex_index = 0;
1062  vertex_index < dof_handler.get_triangulation().n_vertices();
1063  ++vertex_index)
1064  if ((dof_handler.get_triangulation()
1065  .get_used_vertices()[vertex_index] == true) &&
1066  (include_vertex[vertex_index] == true))
1067  {
1068  const unsigned int n_active_fe_indices =
1069  ::internal::DoFAccessorImplementation::Implementation::
1070  n_active_fe_indices(dof_handler,
1071  0,
1072  vertex_index,
1073  std::integral_constant<int, 0>());
1074 
1075  if (n_active_fe_indices > 1)
1076  {
1077  const std::set<types::fe_index> fe_indices =
1080  dof_handler,
1081  0,
1082  vertex_index,
1083  std::integral_constant<int, 0>());
1084 
1085  // find out which is the most dominating finite
1086  // element of the ones that are used on this vertex
1087  // TODO: Change set to types::fe_index
1088  types::fe_index most_dominating_fe_index =
1089  dof_handler.get_fe_collection().find_dominating_fe(
1090  {fe_indices.begin(), fe_indices.end()},
1091  /*codim=*/dim);
1092 
1093  // if we haven't found a dominating finite element,
1094  // choose the very first one to be dominant similar
1095  // to compute_vertex_dof_identities()
1096  if (most_dominating_fe_index == numbers::invalid_fe_index)
1097  most_dominating_fe_index =
1098  ::internal::DoFAccessorImplementation::
1099  Implementation::nth_active_fe_index(
1100  dof_handler,
1101  0,
1102  vertex_index,
1103  0,
1104  std::integral_constant<int, 0>());
1105 
1106  // loop over the indices of all the finite
1107  // elements that are not dominating, and
1108  // identify their dofs to the most dominating
1109  // one
1110  for (const auto &other_fe_index : fe_indices)
1111  if (other_fe_index != most_dominating_fe_index)
1112  {
1113  // make sure the entry in the equivalence
1114  // table exists
1115  const auto &identities =
1116  *ensure_existence_and_return_dof_identities<0>(
1117  dof_handler.get_fe_collection(),
1118  most_dominating_fe_index,
1119  other_fe_index,
1120  vertex_dof_identities[most_dominating_fe_index]
1121  [other_fe_index]);
1122 
1123  // then loop through the identities we
1124  // have. first get the global numbers of the
1125  // dofs we want to identify and make sure they
1126  // are not yet constrained to anything else,
1127  // except for to each other. use the rule that
1128  // we will always constrain the dof with the
1129  // higher FE index to the one with the lower,
1130  // to avoid circular reasoning.
1131  for (const auto &identity : identities)
1132  {
1133  const types::global_dof_index primary_dof_index =
1134  ::internal::DoFAccessorImplementation::
1135  Implementation::get_dof_index(
1136  dof_handler,
1137  0,
1138  vertex_index,
1139  most_dominating_fe_index,
1140  identity.first,
1141  std::integral_constant<int, 0>());
1143  dependent_dof_index =
1144  ::internal::DoFAccessorImplementation::
1145  Implementation::get_dof_index(
1146  dof_handler,
1147  0,
1148  vertex_index,
1149  other_fe_index,
1150  identity.second,
1151  std::integral_constant<int, 0>());
1152 
1153  // check if we are on an interface between
1154  // a locally owned and a ghost cell on which
1155  // we need to work on.
1156  //
1157  // all degrees of freedom belonging to
1158  // dominating FE indices or to a processor
1159  // with a higher rank have been set at this
1160  // point (either in Phase 2, or after the
1161  // first ghost exchange in Phase 5). thus,
1162  // we only have to set the indices of
1163  // degrees of freedom that have been
1164  // previously flagged invalid.
1165  if ((dependent_dof_index ==
1167  (primary_dof_index !=
1169  ::internal::DoFAccessorImplementation::
1170  Implementation::set_dof_index(
1171  dof_handler,
1172  0,
1173  vertex_index,
1174  other_fe_index,
1175  identity.second,
1176  std::integral_constant<int, 0>(),
1177  primary_dof_index);
1178  }
1179  }
1180  }
1181  }
1182  }
1183 
1184 
1185 
1190  template <int spacedim>
1191  static void
1193  DoFHandler<1, spacedim> &dof_handler)
1194  {
1195  (void)dof_handler;
1196  Assert(dof_handler.hp_capability_enabled == true,
1198  }
1199 
1200 
1201  template <int dim, int spacedim>
1202  static void
1204  DoFHandler<dim, spacedim> &dof_handler)
1205  {
1206  Assert(
1207  dof_handler.hp_capability_enabled == true,
1209 
1210  // mark all lines on ghost cells
1211  std::vector<bool> line_marked(
1212  dof_handler.get_triangulation().n_raw_lines());
1213  for (const auto &cell : dof_handler.active_cell_iterators())
1214  if (cell->is_ghost())
1215  for (const auto l : cell->line_indices())
1216  line_marked[cell->line(l)->index()] = true;
1217 
1218  // An implementation of the algorithm described in the hp-paper,
1219  // including the modification mentioned later in the "complications in
1220  // 3-d" subsections
1221  //
1222  // as explained there, we do something only if there are exactly 2
1223  // finite elements associated with an object. if there is only one,
1224  // then there is nothing to do anyway, and if there are 3 or more,
1225  // then we can get into trouble. note that this only happens for lines
1226  // in 3d and higher, and for quads only in 4d and higher, so this
1227  // isn't a particularly frequent case
1228  //
1229  // there is one case, however, that we would like to handle (see, for
1230  // example, the hp/crash_15 testcase): if we have
1231  // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
1232  // then we should be able to handle this because we can simply unify
1233  // *all* dofs, not only a some. so what we do is to first treat all
1234  // pairs of finite elements that have *identical* dofs, and then only
1235  // deal with those that are not identical of which we can handle at
1236  // most 2
1237  ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
1238  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
1239 
1240  for (const auto &cell : dof_handler.active_cell_iterators())
1241  for (const auto l : cell->line_indices())
1242  if ((cell->is_locally_owned()) &&
1243  line_marked[cell->line(l)->index()])
1244  {
1245  const auto line = cell->line(l);
1246  line_marked[line->index()] = false;
1247 
1248  unsigned int unique_sets_of_dofs =
1249  line->n_active_fe_indices();
1250 
1251  // do a first loop over all sets of dofs and do identity
1252  // uniquification
1253  const unsigned int n_active_fe_indices =
1254  line->n_active_fe_indices();
1255  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1256  for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
1257  {
1258  const types::fe_index fe_index_1 =
1259  line->nth_active_fe_index(f),
1260  fe_index_2 =
1261  line->nth_active_fe_index(g);
1262 
1263  if ((dof_handler.get_fe(fe_index_1).n_dofs_per_line() ==
1264  dof_handler.get_fe(fe_index_2)
1265  .n_dofs_per_line()) &&
1266  (dof_handler.get_fe(fe_index_1).n_dofs_per_line() >
1267  0))
1268  {
1269  // the number of dofs per line is identical
1270  const unsigned int dofs_per_line =
1271  dof_handler.get_fe(fe_index_1).n_dofs_per_line();
1272 
1273  const auto &identities =
1274  *ensure_existence_and_return_dof_identities<1>(
1275  dof_handler.get_fe_collection(),
1276  fe_index_1,
1277  fe_index_2,
1278  line_dof_identities[fe_index_1][fe_index_2]);
1279  // see if these sets of dofs are identical. the
1280  // first condition for this is that indeed there are
1281  // n identities
1282  if (identities.size() == dofs_per_line)
1283  {
1284  unsigned int i = 0;
1285  for (; i < dofs_per_line; ++i)
1286  if ((identities[i].first != i) &&
1287  (identities[i].second != i))
1288  // not an identity
1289  break;
1290 
1291  if (i == dofs_per_line)
1292  {
1293  // The line dofs (i.e., the ones interior to
1294  // a line) of these two finite elements are
1295  // identical. Note that there could be
1296  // situations when one element still
1297  // dominates another, e.g.: FE_Q(2) x
1298  // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
1299 
1300  --unique_sets_of_dofs;
1301 
1302  // determine which one of both finite
1303  // elements is the dominating one.
1304  const std::set<types::fe_index> fe_indices{
1305  fe_index_1, fe_index_2};
1306 
1307  // TODO: Change set to types::fe_index
1308  types::fe_index dominating_fe_index =
1309  dof_handler.get_fe_collection()
1310  .find_dominating_fe({fe_indices.begin(),
1311  fe_indices.end()},
1312  /*codim*/ dim - 1);
1313  types::fe_index other_fe_index =
1315 
1316  if (dominating_fe_index !=
1318  other_fe_index =
1319  (dominating_fe_index == fe_index_1) ?
1320  fe_index_2 :
1321  fe_index_1;
1322  else
1323  {
1324  // if we haven't found a dominating
1325  // finite element, choose the one with
1326  // the lower index to be dominating
1327  dominating_fe_index = fe_index_1;
1328  other_fe_index = fe_index_2;
1329  }
1330 
1331  for (unsigned int j = 0; j < dofs_per_line;
1332  ++j)
1333  {
1335  primary_dof_index = line->dof_index(
1336  j, dominating_fe_index);
1338  dependent_dof_index =
1339  line->dof_index(j, other_fe_index);
1340 
1341  // check if we are on an interface
1342  // between a locally owned and a ghost
1343  // cell on which we need to work on.
1344  //
1345  // all degrees of freedom belonging to
1346  // dominating fe_indices or to a
1347  // processor with a higher rank have
1348  // been set at this point (either in
1349  // Phase 2, or after the first ghost
1350  // exchange in Phase 5). thus, we only
1351  // have to set the indices of degrees
1352  // of freedom that have been previously
1353  // flagged invalid.
1354  if ((dependent_dof_index ==
1356  (primary_dof_index !=
1358  line->set_dof_index(j,
1359  primary_dof_index,
1360  other_fe_index);
1361  }
1362  }
1363  }
1364  }
1365  }
1366 
1367  // if at this point, there is only one unique set of dofs
1368  // left, then we have taken care of everything above. if there
1369  // are two, then we need to deal with them here. if there are
1370  // more, then we punt, as described in the paper (and
1371  // mentioned above)
1372  // TODO: The check for 'dim==2' was inserted by intuition. It
1373  // fixes
1374  // the previous problems with @ref step_27 "step-27" in 3d. But an
1375  // explanation for this is still required, and what we do here
1376  // is not what we describe in the paper!.
1377  if ((unique_sets_of_dofs == 2) && (dim == 2))
1378  {
1379  const std::set<types::fe_index> fe_indices =
1380  line->get_active_fe_indices();
1381 
1382  // find out which is the most dominating finite element of
1383  // the ones that are used on this line
1384  // TODO: Change set to types::fe_index
1385  const types::fe_index most_dominating_fe_index =
1386  dof_handler.get_fe_collection().find_dominating_fe(
1387  {fe_indices.begin(), fe_indices.end()},
1388  /*codim=*/dim - 1);
1389 
1390  // if we found the most dominating element, then use this
1391  // to eliminate some of the degrees of freedom by
1392  // identification. otherwise, the code that computes
1393  // hanging node constraints will have to deal with it by
1394  // computing appropriate constraints along this face/edge
1395  if (most_dominating_fe_index != numbers::invalid_fe_index)
1396  {
1397  // loop over the indices of all the finite elements
1398  // that are not dominating, and identify their dofs to
1399  // the most dominating one
1400  for (const auto &other_fe_index : fe_indices)
1401  if (other_fe_index != most_dominating_fe_index)
1402  {
1403  const auto &identities =
1404  *ensure_existence_and_return_dof_identities<
1405  1>(dof_handler.get_fe_collection(),
1406  most_dominating_fe_index,
1407  other_fe_index,
1408  line_dof_identities
1409  [most_dominating_fe_index]
1410  [other_fe_index]);
1411 
1412  for (const auto &identity : identities)
1413  {
1415  primary_dof_index = line->dof_index(
1416  identity.first,
1417  most_dominating_fe_index);
1419  dependent_dof_index =
1420  line->dof_index(identity.second,
1421  other_fe_index);
1422 
1423  // check if we are on an interface between
1424  // a locally owned and a ghost cell on which
1425  // we need to work on.
1426  //
1427  // all degrees of freedom belonging to
1428  // dominating FE indices or to a processor
1429  // with a higher rank have been set at this
1430  // point (either in Phase 2, or after the
1431  // first ghost exchange in Phase 5). thus,
1432  // we only have to set the indices of
1433  // degrees of freedom that have been
1434  // previously flagged invalid.
1435  if ((dependent_dof_index ==
1437  (primary_dof_index !=
1439  line->set_dof_index(identity.second,
1440  primary_dof_index,
1441  other_fe_index);
1442  }
1443  }
1444  }
1445  }
1446  }
1447  }
1448 
1449 
1450 
1455  template <int dim, int spacedim>
1456  static void
1458  DoFHandler<dim, spacedim> &dof_handler)
1459  {
1460  (void)dof_handler;
1461  Assert(
1462  dof_handler.hp_capability_enabled == true,
1464 
1465  // this function should only be called for dim<3 where there are
1466  // no quad dof identities. for dim>=3, the specialization below should
1467  // take care of it
1468  Assert(dim < 3, ExcInternalError());
1469  }
1470 
1471 
1472  template <int spacedim>
1473  static void
1475  DoFHandler<3, spacedim> &dof_handler)
1476  {
1477  Assert(dof_handler.hp_capability_enabled == true,
1479 
1480  const int dim = 3;
1481 
1482  // mark all quads on ghost cells
1483  std::vector<bool> quad_marked(
1484  dof_handler.get_triangulation().n_raw_quads());
1485  for (const auto &cell : dof_handler.active_cell_iterators())
1486  if (cell->is_ghost())
1487  for (const auto q : cell->face_indices())
1488  quad_marked[cell->quad(q)->index()] = true;
1489 
1490  // An implementation of the algorithm described in the hp-
1491  // paper, including the modification mentioned later in the
1492  // "complications in 3-d" subsections
1493  //
1494  // as explained there, we do something only if there are
1495  // exactly 2 finite elements associated with an object. if
1496  // there is only one, then there is nothing to do anyway,
1497  // and if there are 3 or more, then we can get into
1498  // trouble. note that this only happens for lines in 3d and
1499  // higher, and for quads only in 4d and higher, so this
1500  // isn't a particularly frequent case
1501  ::Table<3, std::unique_ptr<DoFIdentities>> quad_dof_identities(
1502  dof_handler.fe_collection.size(),
1503  dof_handler.fe_collection.size(),
1504  2 /*triangle (0) or quadrilateral (1)*/);
1505 
1506  for (const auto &cell : dof_handler.active_cell_iterators())
1507  for (const auto q : cell->face_indices())
1508  if ((cell->is_locally_owned()) &&
1509  quad_marked[cell->quad(q)->index()] &&
1510  (cell->quad(q)->n_active_fe_indices() == 2))
1511  {
1512  const auto quad = cell->quad(q);
1513  quad_marked[quad->index()] = false;
1514 
1515  const std::set<types::fe_index> fe_indices =
1516  quad->get_active_fe_indices();
1517 
1518  // find out which is the most dominating finite
1519  // element of the ones that are used on this quad
1520  // TODO: Change set to types::fe_index
1521  const types::fe_index most_dominating_fe_index =
1522  dof_handler.get_fe_collection().find_dominating_fe(
1523  {fe_indices.begin(), fe_indices.end()},
1524  /*codim=*/dim - 2);
1525 
1526  const types::fe_index most_dominating_fe_index_face_no =
1527  cell->active_fe_index() == most_dominating_fe_index ?
1528  q :
1529  cell->neighbor_face_no(q);
1530 
1531  // if we found the most dominating element, then use
1532  // this to eliminate some of the degrees of freedom
1533  // by identification. otherwise, the code that
1534  // computes hanging node constraints will have to
1535  // deal with it by computing appropriate constraints
1536  // along this face/edge
1537  if (most_dominating_fe_index != numbers::invalid_fe_index)
1538  {
1539  // loop over the indices of all the finite
1540  // elements that are not dominating, and
1541  // identify their dofs to the most dominating
1542  // one
1543  for (const auto &other_fe_index : fe_indices)
1544  if (other_fe_index != most_dominating_fe_index)
1545  {
1546  const auto &identities =
1547  *ensure_existence_and_return_dof_identities<2>(
1548  dof_handler.get_fe_collection(),
1549  most_dominating_fe_index,
1550  other_fe_index,
1551  quad_dof_identities
1552  [most_dominating_fe_index][other_fe_index]
1553  [cell->quad(q)->reference_cell() ==
1555  most_dominating_fe_index_face_no);
1556 
1557  for (const auto &identity : identities)
1558  {
1560  primary_dof_index =
1561  quad->dof_index(identity.first,
1562  most_dominating_fe_index);
1564  dependent_dof_index =
1565  quad->dof_index(identity.second,
1566  other_fe_index);
1567 
1568  // check if we are on an interface between
1569  // a locally owned and a ghost cell on which
1570  // we need to work on.
1571  //
1572  // all degrees of freedom belonging to
1573  // dominating FE indices or to a processor with
1574  // a higher rank have been set at this point
1575  // (either in Phase 2, or after the first ghost
1576  // exchange in Phase 5). thus, we only have to
1577  // set the indices of degrees of freedom that
1578  // have been previously flagged invalid.
1579  if ((dependent_dof_index ==
1581  (primary_dof_index !=
1583  quad->set_dof_index(identity.second,
1584  primary_dof_index,
1585  other_fe_index);
1586  }
1587  }
1588  }
1589  }
1590  }
1591 
1592 
1593 
1606  template <int dim, int spacedim>
1607  static void
1609  DoFHandler<dim, spacedim> &dof_handler)
1610  {
1611  if (dof_handler.hp_capability_enabled == false)
1612  return;
1613 
1614  {
1615  Threads::TaskGroup<> tasks;
1616 
1617  tasks += Threads::new_task([&]() {
1619  });
1620 
1621  if (dim > 1)
1622  {
1623  tasks += Threads::new_task([&]() {
1625  });
1626  }
1627 
1628  if (dim > 2)
1629  {
1630  tasks += Threads::new_task([&]() {
1632  });
1633  }
1634 
1635  tasks.join_all();
1636  }
1637  }
1638 
1639 
1640 
1647  template <int dim, int spacedim>
1650  DoFHandler<dim, spacedim> &dof_handler)
1651  {
1652  Assert(dof_handler.get_triangulation().n_levels() > 0,
1653  ExcMessage("Empty triangulation"));
1654 
1655  // distribute dofs on all cells excluding artificial ones
1656  types::global_dof_index next_free_dof = 0;
1657 
1658  for (auto cell : dof_handler.active_cell_iterators())
1659  if (!cell->is_artificial() &&
1661  (cell->subdomain_id() == subdomain_id)))
1662  {
1663  // feed the process_dof_indices function with an empty type
1664  // `std::tuple<>`, as we do not want to retrieve any DoF
1665  // indices here and rather modify the stored ones
1666  DoFAccessorImplementation::Implementation::process_dof_indices(
1667  *cell,
1668  std::make_tuple(),
1669  cell->active_fe_index(),
1670  DoFAccessorImplementation::Implementation::
1671  DoFIndexProcessor<dim, spacedim>(),
1672  [&next_free_dof](auto &stored_index, auto) {
1673  if (stored_index == numbers::invalid_dof_index)
1674  {
1675  stored_index = next_free_dof;
1676  Assert(
1677  next_free_dof !=
1678  std::numeric_limits<types::global_dof_index>::max(),
1679  ExcMessage(
1680  "You have reached the maximal number of degrees of "
1681  "freedom that can be stored in the chosen data "
1682  "type. In practice, this can only happen if you "
1683  "are using 32-bit data types. You will have to "
1684  "re-compile deal.II with the "
1685  "`DEAL_II_WITH_64BIT_INDICES' flag set to `ON'."));
1686  ++next_free_dof;
1687  }
1688  },
1689  false);
1690  }
1691 
1692  return next_free_dof;
1693  }
1694 
1695 
1696 
1710  template <int dim, int spacedim>
1711  static void
1713  std::vector<types::global_dof_index> &renumbering,
1715  const DoFHandler<dim, spacedim> &dof_handler)
1716  {
1717  std::vector<types::global_dof_index> local_dof_indices;
1718 
1719  for (const auto &cell : dof_handler.active_cell_iterators())
1720  if (cell->is_ghost() && (cell->subdomain_id() < subdomain_id))
1721  {
1722  // we found a neighboring ghost cell whose subdomain
1723  // is "stronger" than our own subdomain
1724 
1725  // delete all dofs that live there and that we have
1726  // previously assigned a number to (i.e. the ones on
1727  // the interface); make sure to not use the cache
1728  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1729  internal::DoFAccessorImplementation::Implementation::
1730  get_dof_indices(*cell,
1731  local_dof_indices,
1732  cell->active_fe_index());
1733  for (const auto &local_dof_index : local_dof_indices)
1734  if (local_dof_index != numbers::invalid_dof_index)
1735  renumbering[local_dof_index] = numbers::invalid_dof_index;
1736  }
1737  }
1738 
1739 
1740 
1741  /* -------------- distribute_mg_dofs functionality ------------- */
1742 
1743 
1744 
1745  template <int dim, int spacedim>
1748  DoFHandler<dim, spacedim> &dof_handler,
1749  const unsigned int level)
1750  {
1751  Assert(dof_handler.hp_capability_enabled == false,
1752  ExcInternalError());
1753 
1754  const ::Triangulation<dim, spacedim> &tria =
1755  dof_handler.get_triangulation();
1756  Assert(tria.n_levels() > 0, ExcMessage("Empty triangulation"));
1757  if (level >= tria.n_levels())
1758  return 0; // this is allowed for multigrid
1759 
1760  types::global_dof_index next_free_dof = 0;
1761 
1762  for (auto cell : dof_handler.cell_iterators_on_level(level))
1763  if ((level_subdomain_id == numbers::invalid_subdomain_id) ||
1764  (cell->level_subdomain_id() == level_subdomain_id))
1765  {
1766  DoFAccessorImplementation::Implementation::process_dof_indices(
1767  *cell,
1768  std::make_tuple(),
1769  0,
1770  DoFAccessorImplementation::Implementation::
1771  MGDoFIndexProcessor<dim, spacedim>(level),
1772  [&next_free_dof](auto &stored_index, auto) {
1773  if (stored_index == numbers::invalid_dof_index)
1774  {
1775  stored_index = next_free_dof;
1776  Assert(
1777  next_free_dof !=
1779  ExcMessage(
1780  "You have reached the maximal number of degrees of "
1781  "freedom that can be stored in the chosen data "
1782  "type. In practice, this can only happen if you "
1783  "are using 32-bit data types. You will have to "
1784  "re-compile deal.II with the "
1785  "`DEAL_II_WITH_64BIT_INDICES' flag set to `ON'."));
1786  ++next_free_dof;
1787  }
1788  },
1789  true);
1790  }
1791 
1792  return next_free_dof;
1793  }
1794 
1795 
1796 
1797  /* --------------------- renumber_dofs functionality ---------------- */
1798 
1799 
1807  template <int dim, int spacedim>
1808  static void
1810  const std::vector<types::global_dof_index> &new_numbers,
1811  const IndexSet &indices_we_care_about,
1812  DoFHandler<dim, spacedim> &dof_handler)
1813  {
1814  for (unsigned int d = 1; d < dim; ++d)
1815  for (auto &i : dof_handler.object_dof_indices[0][d])
1816  if (i != numbers::invalid_dof_index)
1817  i = ((indices_we_care_about.size() == 0) ?
1818  new_numbers[i] :
1819  new_numbers[indices_we_care_about.index_within_set(i)]);
1820  }
1821 
1822 
1823 
1824  template <int dim, int spacedim>
1825  static void
1827  const std::vector<types::global_dof_index> &new_numbers,
1828  const IndexSet &indices_we_care_about,
1829  DoFHandler<dim, spacedim> &dof_handler,
1830  const bool check_validity)
1831  {
1832  if (dof_handler.hp_capability_enabled == false)
1833  {
1834  // we can not use cell iterators in this function since then
1835  // we would renumber the dofs on the interface of two cells
1836  // more than once. Anyway, this way it's not only more
1837  // correct but also faster; note, however, that dof numbers
1838  // may be invalid_dof_index, namely when the appropriate
1839  // vertex/line/etc is unused
1840  for (std::vector<types::global_dof_index>::iterator i =
1841  dof_handler.object_dof_indices[0][0].begin();
1842  i != dof_handler.object_dof_indices[0][0].end();
1843  ++i)
1844  if (*i != numbers::invalid_dof_index)
1845  *i =
1846  (indices_we_care_about.size() == 0) ?
1847  (new_numbers[*i]) :
1848  (new_numbers[indices_we_care_about.index_within_set(*i)]);
1849  else if (check_validity)
1850  // if index is invalid_dof_index: check if this one
1851  // really is unused
1852  Assert(dof_handler.get_triangulation().vertex_used(
1853  (i - dof_handler.object_dof_indices[0][0].begin()) /
1854  dof_handler.get_fe().n_dofs_per_vertex()) == false,
1855  ExcInternalError());
1856  return;
1857  }
1858 
1859 
1860  for (unsigned int vertex_index = 0;
1861  vertex_index < dof_handler.get_triangulation().n_vertices();
1862  ++vertex_index)
1863  {
1864  const unsigned int n_active_fe_indices =
1865  ::internal::DoFAccessorImplementation::Implementation::
1866  n_active_fe_indices(dof_handler,
1867  0,
1868  vertex_index,
1869  std::integral_constant<int, 0>());
1870 
1871  // if this vertex is unused, then we really ought not to have
1872  // allocated any space for it, i.e., n_active_fe_indices should be
1873  // zero, and there is no space to actually store dof indices for
1874  // this vertex
1875  if (dof_handler.get_triangulation().vertex_used(vertex_index) ==
1876  false)
1877  Assert(n_active_fe_indices == 0, ExcInternalError());
1878 
1879  // otherwise the vertex is used; it may still not hold any dof
1880  // indices if it is located on an artificial cell and not adjacent
1881  // to a ghost cell, but in that case there is simply nothing for
1882  // us to do
1883  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1884  {
1885  const types::fe_index fe_index =
1886  ::internal::DoFAccessorImplementation::
1887  Implementation::nth_active_fe_index(
1888  dof_handler,
1889  0,
1890  vertex_index,
1891  f,
1892  std::integral_constant<int, 0>());
1893 
1894  for (unsigned int d = 0;
1895  d < dof_handler.get_fe(fe_index).n_dofs_per_vertex();
1896  ++d)
1897  {
1898  const types::global_dof_index old_dof_index =
1899  ::internal::DoFAccessorImplementation::
1900  Implementation::get_dof_index(
1901  dof_handler,
1902  0,
1903  vertex_index,
1904  fe_index,
1905  d,
1906  std::integral_constant<int, 0>());
1907 
1908  // if check_validity was set, then we are to verify that
1909  // the previous indices were all valid. this really should
1910  // be the case: we allocated space for these vertex dofs,
1911  // i.e., at least one adjacent cell has a valid
1912  // active FE index, so there are DoFs that really live
1913  // on this vertex. if check_validity is set, then we
1914  // must make sure that they have been set to something
1915  // useful
1916  if (check_validity)
1917  Assert(old_dof_index != numbers::invalid_dof_index,
1918  ExcInternalError());
1919 
1920  if (old_dof_index != numbers::invalid_dof_index)
1921  {
1922  // In the following blocks, we first check whether
1923  // we were given an IndexSet of DoFs to touch. If not
1924  // (the first 'if' case here), then we are in the
1925  // sequential case and are allowed to touch all DoFs.
1926  //
1927  // If yes (the 'else' case), then we need to
1928  // distinguish whether the DoF whose number we want to
1929  // touch is in fact locally owned (i.e., is in the
1930  // index set) and then we can actually assign it a new
1931  // number; otherwise, we have encountered a
1932  // non-locally owned DoF for which we don't know the
1933  // new number yet and so set it to an invalid index.
1934  // This will later be fixed up after the first ghost
1935  // exchange phase when we unify hp-DoFs on neighboring
1936  // cells.
1937  if (indices_we_care_about.size() == 0)
1938  ::internal::DoFAccessorImplementation::
1939  Implementation::set_dof_index(
1940  dof_handler,
1941  0,
1942  vertex_index,
1943  fe_index,
1944  d,
1945  std::integral_constant<int, 0>(),
1946  new_numbers[old_dof_index]);
1947  else
1948  {
1949  if (indices_we_care_about.is_element(
1950  old_dof_index))
1951  ::internal::DoFAccessorImplementation::
1952  Implementation::set_dof_index(
1953  dof_handler,
1954  0,
1955  vertex_index,
1956  fe_index,
1957  d,
1958  std::integral_constant<int, 0>(),
1959  new_numbers[indices_we_care_about
1960  .index_within_set(
1961  old_dof_index)]);
1962  else
1963  ::internal::DoFAccessorImplementation::
1964  Implementation::set_dof_index(
1965  dof_handler,
1966  0,
1967  vertex_index,
1968  fe_index,
1969  d,
1970  std::integral_constant<int, 0>(),
1972  }
1973  }
1974  }
1975  }
1976  }
1977  }
1978 
1979 
1980 
1981  template <int dim, int spacedim>
1982  static void
1984  const std::vector<types::global_dof_index> &new_numbers,
1985  const IndexSet &indices_we_care_about,
1986  DoFHandler<dim, spacedim> &dof_handler)
1987  {
1988  if (dof_handler.hp_capability_enabled == false)
1989  {
1990  for (unsigned int level = 0;
1991  level < dof_handler.object_dof_indices.size();
1992  ++level)
1993  for (auto &i : dof_handler.object_dof_indices[level][dim])
1994  if (i != numbers::invalid_dof_index)
1995  i = ((indices_we_care_about.size() == 0) ?
1996  new_numbers[i] :
1997  new_numbers[indices_we_care_about.index_within_set(
1998  i)]);
1999  return;
2000  }
2001 
2002  for (const auto &cell : dof_handler.active_cell_iterators())
2003  if (!cell->is_artificial())
2004  {
2005  const types::fe_index fe_index = cell->active_fe_index();
2006 
2007  for (unsigned int d = 0;
2008  d < dof_handler.get_fe(fe_index)
2009  .template n_dofs_per_object<dim>();
2010  ++d)
2011  {
2012  const types::global_dof_index old_dof_index =
2013  cell->dof_index(d, fe_index);
2014  if (old_dof_index != numbers::invalid_dof_index)
2015  {
2016  // In the following blocks, we first check whether
2017  // we were given an IndexSet of DoFs to touch. If not
2018  // (the first 'if' case here), then we are in the
2019  // sequential case and are allowed to touch all DoFs.
2020  //
2021  // If yes (the 'else' case), then we need to distinguish
2022  // whether the DoF whose number we want to touch is in
2023  // fact locally owned (i.e., is in the index set) and
2024  // then we can actually assign it a new number;
2025  // otherwise, we have encountered a non-locally owned
2026  // DoF for which we don't know the new number yet and so
2027  // set it to an invalid index. This will later be fixed
2028  // up after the first ghost exchange phase when we unify
2029  // hp-DoFs on neighboring cells.
2030  if (indices_we_care_about.size() == 0)
2031  cell->set_dof_index(d,
2032  new_numbers[old_dof_index],
2033  fe_index);
2034  else
2035  {
2036  if (indices_we_care_about.is_element(old_dof_index))
2037  cell->set_dof_index(
2038  d,
2039  new_numbers[indices_we_care_about
2040  .index_within_set(old_dof_index)],
2041  fe_index);
2042  else
2043  cell->set_dof_index(d,
2045  fe_index);
2046  }
2047  }
2048  }
2049  }
2050  }
2051 
2052 
2053 
2054  template <int spacedim>
2055  static void
2057  const std::vector<types::global_dof_index> & /*new_numbers*/,
2058  const IndexSet & /*indices_we_care_about*/,
2059  DoFHandler<1, spacedim> & /*dof_handler*/)
2060  {
2061  // nothing to do in 1d since there are no separate faces -- we've
2062  // already taken care of this when dealing with the vertices
2063  }
2064 
2065 
2066 
2067  template <int spacedim>
2068  static void
2070  const std::vector<types::global_dof_index> &new_numbers,
2071  const IndexSet &indices_we_care_about,
2072  DoFHandler<2, spacedim> &dof_handler)
2073  {
2074  const unsigned int dim = 2;
2075 
2076  if (dof_handler.hp_capability_enabled == false)
2077  {
2078  for (unsigned int d = 1; d < dim; ++d)
2079  for (auto &i : dof_handler.object_dof_indices[0][d])
2080  if (i != numbers::invalid_dof_index)
2081  i = ((indices_we_care_about.size() == 0) ?
2082  new_numbers[i] :
2083  new_numbers[indices_we_care_about.index_within_set(
2084  i)]);
2085  return;
2086  }
2087 
2088  // deal with DoFs on lines
2089  {
2090  std::vector<bool> line_touched(
2091  dof_handler.get_triangulation().n_raw_lines());
2092  for (const auto &cell : dof_handler.active_cell_iterators())
2093  if (!cell->is_artificial())
2094  for (const auto l : cell->line_indices())
2095  if (!line_touched[cell->line(l)->index()])
2096  {
2097  const auto line = cell->line(l);
2098  line_touched[line->index()] = true;
2099 
2100  const unsigned int n_active_fe_indices =
2101  line->n_active_fe_indices();
2102 
2103  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2104  {
2105  const types::fe_index fe_index =
2106  line->nth_active_fe_index(f);
2107 
2108  for (unsigned int d = 0;
2109  d <
2110  dof_handler.get_fe(fe_index).n_dofs_per_line();
2111  ++d)
2112  {
2113  const types::global_dof_index old_dof_index =
2114  line->dof_index(d, fe_index);
2115  if (old_dof_index != numbers::invalid_dof_index)
2116  {
2117  // In the following blocks, we first check
2118  // whether we were given an IndexSet of DoFs
2119  // to touch. If not (the first 'if' case
2120  // here), then we are in the sequential case
2121  // and are allowed to touch all DoFs.
2122  //
2123  // If yes (the 'else' case), then we need to
2124  // distinguish whether the DoF whose number we
2125  // want to touch is in fact locally owned
2126  // (i.e., is in the index set) and then we can
2127  // actually assign it a new number; otherwise,
2128  // we have encountered a non-locally owned DoF
2129  // for which we don't know the new number yet
2130  // and so set it to an invalid index. This
2131  // will later be fixed up after the first
2132  // ghost exchange phase when we unify hp-DoFs
2133  // on neighboring cells.
2134  if (indices_we_care_about.size() == 0)
2135  line->set_dof_index(
2136  d, new_numbers[old_dof_index], fe_index);
2137  else
2138  {
2139  if (indices_we_care_about.is_element(
2140  old_dof_index))
2141  line->set_dof_index(
2142  d,
2143  new_numbers[indices_we_care_about
2144  .index_within_set(
2145  old_dof_index)],
2146  fe_index);
2147  else
2148  line->set_dof_index(
2149  d,
2151  fe_index);
2152  }
2153  }
2154  }
2155  }
2156  }
2157  }
2158  }
2159 
2160 
2161 
2162  template <int spacedim>
2163  static void
2165  const std::vector<types::global_dof_index> &new_numbers,
2166  const IndexSet &indices_we_care_about,
2167  DoFHandler<3, spacedim> &dof_handler)
2168  {
2169  const unsigned int dim = 3;
2170 
2171  if (dof_handler.hp_capability_enabled == false)
2172  {
2173  for (unsigned int d = 1; d < dim; ++d)
2174  for (auto &i : dof_handler.object_dof_indices[0][d])
2175  if (i != numbers::invalid_dof_index)
2176  i = ((indices_we_care_about.size() == 0) ?
2177  new_numbers[i] :
2178  new_numbers[indices_we_care_about.index_within_set(
2179  i)]);
2180  return;
2181  }
2182 
2183  // deal with DoFs on lines
2184  {
2185  std::vector<bool> line_touched(
2186  dof_handler.get_triangulation().n_raw_lines());
2187  for (const auto &cell : dof_handler.active_cell_iterators())
2188  if (!cell->is_artificial())
2189  for (const auto l : cell->line_indices())
2190  if (!line_touched[cell->line(l)->index()])
2191  {
2192  const auto line = cell->line(l);
2193  line_touched[line->index()] = true;
2194 
2195  const unsigned int n_active_fe_indices =
2196  line->n_active_fe_indices();
2197 
2198  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2199  {
2200  const types::fe_index fe_index =
2201  line->nth_active_fe_index(f);
2202 
2203  for (unsigned int d = 0;
2204  d <
2205  dof_handler.get_fe(fe_index).n_dofs_per_line();
2206  ++d)
2207  {
2208  const types::global_dof_index old_dof_index =
2209  line->dof_index(d, fe_index);
2210  if (old_dof_index != numbers::invalid_dof_index)
2211  {
2212  // In the following blocks, we first check
2213  // whether we were given an IndexSet of DoFs
2214  // to touch. If not (the first 'if' case
2215  // here), then we are in the sequential case
2216  // and are allowed to touch all DoFs.
2217  //
2218  // If yes (the 'else' case), then we need to
2219  // distinguish whether the DoF whose number we
2220  // want to touch is in fact locally owned
2221  // (i.e., is in the index set) and then we can
2222  // actually assign it a new number; otherwise,
2223  // we have encountered a non-locally owned DoF
2224  // for which we don't know the new number yet
2225  // and so set it to an invalid index. This
2226  // will later be fixed up after the first
2227  // ghost exchange phase when we unify hp-DoFs
2228  // on neighboring cells.
2229  if (indices_we_care_about.size() == 0)
2230  line->set_dof_index(
2231  d, new_numbers[old_dof_index], fe_index);
2232  else if (indices_we_care_about.is_element(
2233  old_dof_index))
2234  line->set_dof_index(
2235  d,
2236  new_numbers[indices_we_care_about
2237  .index_within_set(
2238  old_dof_index)],
2239  fe_index);
2240  else
2241  line->set_dof_index(
2243  }
2244  }
2245  }
2246  }
2247  }
2248 
2249  // then deal with dofs on quads
2250  {
2251  std::vector<bool> quad_touched(
2252  dof_handler.get_triangulation().n_raw_quads());
2253  for (const auto &cell : dof_handler.active_cell_iterators())
2254  if (!cell->is_artificial())
2255  for (const auto q : cell->face_indices())
2256  if (!quad_touched[cell->quad(q)->index()])
2257  {
2258  const auto quad = cell->quad(q);
2259  quad_touched[quad->index()] = true;
2260 
2261  const unsigned int n_active_fe_indices =
2262  quad->n_active_fe_indices();
2263 
2264  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2265  {
2266  const types::fe_index fe_index =
2267  quad->nth_active_fe_index(f);
2268 
2269  for (unsigned int d = 0;
2270  d <
2271  dof_handler.get_fe(fe_index).n_dofs_per_quad(q);
2272  ++d)
2273  {
2274  const types::global_dof_index old_dof_index =
2275  quad->dof_index(d, fe_index);
2276  if (old_dof_index != numbers::invalid_dof_index)
2277  {
2278  // In the following blocks, we first check
2279  // whether we were given an IndexSet of DoFs
2280  // to touch. If not (the first 'if' case
2281  // here), then we are in the sequential case
2282  // and are allowed to touch all DoFs.
2283  //
2284  // If yes (the 'else' case), then we need to
2285  // distinguish whether the DoF whose number we
2286  // want to touch is in fact locally owned
2287  // (i.e., is in the index set) and then we can
2288  // actually assign it a new number; otherwise,
2289  // we have encountered a non-locally owned DoF
2290  // for which we don't know the new number yet
2291  // and so set it to an invalid index. This
2292  // will later be fixed up after the first
2293  // ghost exchange phase when we unify hp-DoFs
2294  // on neighboring cells.
2295  if (indices_we_care_about.size() == 0)
2296  quad->set_dof_index(
2297  d, new_numbers[old_dof_index], fe_index);
2298  else
2299  {
2300  if (indices_we_care_about.is_element(
2301  old_dof_index))
2302  quad->set_dof_index(
2303  d,
2304  new_numbers[indices_we_care_about
2305  .index_within_set(
2306  old_dof_index)],
2307  fe_index);
2308  else
2309  quad->set_dof_index(
2310  d,
2312  fe_index);
2313  }
2314  }
2315  }
2316  }
2317  }
2318  }
2319  }
2320 
2321 
2322 
2334  template <int dim, int space_dim>
2335  static void
2336  renumber_dofs(const std::vector<types::global_dof_index> &new_numbers,
2337  const IndexSet &indices_we_care_about,
2338  const DoFHandler<dim, space_dim> &dof_handler,
2339  const bool check_validity)
2340  {
2341  if (dim == 1)
2342  Assert(indices_we_care_about == IndexSet(0), ExcNotImplemented());
2343 
2344  // renumber DoF indices on vertices, cells, and faces. this
2345  // can be done in parallel because the respective functions
2346  // work on separate data structures
2347  Threads::TaskGroup<> tasks;
2348  tasks += Threads::new_task([&]() {
2349  renumber_vertex_dofs(new_numbers,
2350  indices_we_care_about,
2351  const_cast<DoFHandler<dim, space_dim> &>(
2352  dof_handler),
2353  check_validity);
2354  });
2355  tasks += Threads::new_task([&]() {
2356  renumber_face_dofs(new_numbers,
2357  indices_we_care_about,
2358  const_cast<DoFHandler<dim, space_dim> &>(
2359  dof_handler));
2360  });
2361  tasks += Threads::new_task([&]() {
2362  renumber_cell_dofs(new_numbers,
2363  indices_we_care_about,
2364  const_cast<DoFHandler<dim, space_dim> &>(
2365  dof_handler));
2366  });
2367  tasks.join_all();
2368  }
2369 
2370 
2371 
2372  /* --------------------- renumber_mg_dofs functionality ----------------
2373  */
2374 
2382  template <int dim, int spacedim>
2383  static void
2385  const std::vector<::types::global_dof_index> &new_numbers,
2386  const IndexSet &indices_we_care_about,
2387  DoFHandler<dim, spacedim> &dof_handler,
2388  const unsigned int level)
2389  {
2390  Assert(level < dof_handler.get_triangulation().n_levels(),
2391  ExcInternalError());
2392 
2393  for (auto i = dof_handler.mg_vertex_dofs.begin();
2394  i != dof_handler.mg_vertex_dofs.end();
2395  ++i)
2396  // if the present vertex lives on the current level
2397  if ((i->get_coarsest_level() <= level) &&
2398  (i->get_finest_level() >= level))
2399  for (unsigned int d = 0;
2400  d < dof_handler.get_fe().n_dofs_per_vertex();
2401  ++d)
2402  {
2404  i->access_index(level,
2405  d,
2406  dof_handler.get_fe().n_dofs_per_vertex());
2407 
2408  if (idx != numbers::invalid_dof_index)
2409  {
2410  Assert(indices_we_care_about.size() > 0 ?
2411  indices_we_care_about.is_element(idx) :
2412  (idx < new_numbers.size()),
2413  ExcInternalError());
2414  i->access_index(
2415  level, d, dof_handler.get_fe().n_dofs_per_vertex()) =
2416  (indices_we_care_about.size() == 0) ?
2417  new_numbers[idx] :
2418  new_numbers[indices_we_care_about.index_within_set(
2419  idx)];
2420  }
2421  }
2422  }
2423 
2424 
2425 
2433  template <int dim, int spacedim>
2434  static void
2436  const std::vector<::types::global_dof_index> &new_numbers,
2437  const IndexSet &indices_we_care_about,
2438  DoFHandler<dim, spacedim> &dof_handler,
2439  const unsigned int level)
2440  {
2441  for (std::vector<types::global_dof_index>::iterator i =
2442  dof_handler.mg_levels[level]->dof_object.dofs.begin();
2443  i != dof_handler.mg_levels[level]->dof_object.dofs.end();
2444  ++i)
2445  {
2446  if (*i != numbers::invalid_dof_index)
2447  {
2448  Assert((indices_we_care_about.size() > 0 ?
2449  indices_we_care_about.is_element(*i) :
2450  (*i < new_numbers.size())),
2451  ExcInternalError());
2452  *i =
2453  (indices_we_care_about.size() == 0) ?
2454  (new_numbers[*i]) :
2455  (new_numbers[indices_we_care_about.index_within_set(*i)]);
2456  }
2457  }
2458  }
2459 
2460 
2461 
2469  template <int spacedim>
2470  static void
2472  const std::vector<types::global_dof_index> & /*new_numbers*/,
2473  const IndexSet & /*indices_we_care_about*/,
2474  DoFHandler<1, spacedim> & /*dof_handler*/,
2475  const unsigned int /*level*/,
2476  const bool /*check_validity*/)
2477  {
2478  // nothing to do in 1d because there are no separate faces
2479  }
2480 
2481 
2482 
2483  template <int dim, int spacedim>
2484  static void
2486  const std::vector<::types::global_dof_index> &new_numbers,
2487  const IndexSet &indices_we_care_about,
2488  DoFHandler<dim, spacedim> &dof_handler,
2489  const unsigned int level,
2490  const bool check_validity)
2491  {
2492  const unsigned int dofs_per_line =
2493  dof_handler.get_fe().n_dofs_per_line();
2494  if (dofs_per_line > 0 ||
2495  (dim > 2 && dof_handler.get_fe().max_dofs_per_quad() > 0))
2496  {
2497  // visit all lines/quads adjacent to cells of the current level
2498  // exactly once, as those lines/quads logically belong to the same
2499  // level as the cell, at least for isotropic refinement
2500  std::vector<bool> line_touched(
2501  dof_handler.get_triangulation().n_raw_lines());
2502  std::vector<bool> quad_touched(
2503  dim > 2 ? dof_handler.get_triangulation().n_raw_quads() : 0);
2504  for (const auto &cell :
2505  dof_handler.cell_iterators_on_level(level))
2506  if (cell->level_subdomain_id() !=
2508  {
2509  // lines
2510  if (dofs_per_line > 0)
2511  {
2512  const auto line_indices =
2513  internal::TriaAccessorImplementation::Implementation::
2514  get_line_indices_of_cell(*cell);
2515  for (const auto line : cell->line_indices())
2516  {
2517  if (!line_touched[line_indices[line]])
2518  {
2519  line_touched[line_indices[line]] = true;
2520  ::types::global_dof_index *indices =
2521  &internal::DoFAccessorImplementation::
2522  Implementation::get_mg_dof_index(
2523  dof_handler,
2524  dof_handler.mg_levels[level],
2525  dof_handler.mg_faces,
2526  line_indices[line],
2527  0,
2528  0,
2529  std::integral_constant<int, 1>());
2530  for (unsigned int d = 0; d < dofs_per_line; ++d)
2531  {
2532  if (check_validity)
2533  Assert(indices[d] !=
2535  ExcInternalError());
2536 
2537  if (indices[d] !=
2539  indices[d] =
2540  (indices_we_care_about.size() == 0) ?
2541  new_numbers[indices[d]] :
2542  new_numbers[indices_we_care_about
2543  .index_within_set(
2544  indices[d])];
2545  }
2546  }
2547  }
2548  }
2549 
2550  // quads
2551  if (dim > 2)
2552  for (const auto quad : cell->face_indices())
2553  if (!quad_touched[cell->quad(quad)->index()])
2554  {
2555  quad_touched[cell->quad(quad)->index()] = true;
2556  const unsigned int dofs_per_quad =
2557  dof_handler.get_fe().n_dofs_per_quad(quad);
2558  if (dofs_per_quad > 0)
2559  {
2560  ::types::global_dof_index *indices =
2561  &internal::DoFAccessorImplementation::
2562  Implementation::get_mg_dof_index(
2563  dof_handler,
2564  dof_handler.mg_levels[level],
2565  dof_handler.mg_faces,
2566  cell->quad(quad)->index(),
2567  0,
2568  0,
2569  std::integral_constant<int, 2>());
2570  for (unsigned int d = 0; d < dofs_per_quad; ++d)
2571  {
2572  if (check_validity)
2573  Assert(indices[d] !=
2575  ExcInternalError());
2576 
2577  if (indices[d] !=
2579  indices[d] =
2580  (indices_we_care_about.size() == 0) ?
2581  new_numbers[indices[d]] :
2582  new_numbers[indices_we_care_about
2583  .index_within_set(
2584  indices[d])];
2585  }
2586  }
2587  }
2588  }
2589  }
2590  }
2591 
2592 
2593 
2594  template <int dim, int spacedim>
2595  static void
2597  const std::vector<::types::global_dof_index> &new_numbers,
2598  const IndexSet &indices_we_care_about,
2599  DoFHandler<dim, spacedim> &dof_handler,
2600  const unsigned int level,
2601  const bool check_validity)
2602  {
2603  Assert(
2604  dof_handler.hp_capability_enabled == false,
2606 
2607  Assert(level < dof_handler.get_triangulation().n_global_levels(),
2608  ExcInternalError());
2609 
2610  // renumber DoF indices on vertices, cells, and faces. this
2611  // can be done in parallel because the respective functions
2612  // work on separate data structures
2613  Threads::TaskGroup<> tasks;
2614  tasks += Threads::new_task([&]() {
2615  renumber_vertex_mg_dofs(new_numbers,
2616  indices_we_care_about,
2617  dof_handler,
2618  level);
2619  });
2620  tasks += Threads::new_task([&]() {
2621  renumber_face_mg_dofs(new_numbers,
2622  indices_we_care_about,
2623  dof_handler,
2624  level,
2625  check_validity);
2626  });
2627  tasks += Threads::new_task([&]() {
2628  renumber_cell_mg_dofs(new_numbers,
2629  indices_we_care_about,
2630  dof_handler,
2631  level);
2632  });
2633  tasks.join_all();
2634  }
2635  };
2636 
2637 
2638 
2639  /* --------------------- class Sequential ---------------- */
2640 
2641 
2642 
2643  template <int dim, int spacedim>
2645  DoFHandler<dim, spacedim> &dof_handler)
2646  : dof_handler(&dof_handler)
2647  {}
2648 
2649 
2650 
2651  template <int dim, int spacedim>
2652  NumberCache
2654  {
2655  const types::global_dof_index n_initial_dofs =
2657  *dof_handler);
2658 
2659  const types::global_dof_index n_dofs =
2660  Implementation::unify_dof_indices(*dof_handler,
2661  n_initial_dofs,
2662  /*check_validity=*/true);
2663 
2664  // return a sequential, complete index set
2665  return NumberCache(n_dofs);
2666  }
2667 
2668 
2669 
2670  template <int dim, int spacedim>
2671  std::vector<NumberCache>
2673  {
2674  std::vector<NumberCache> number_caches;
2675  number_caches.reserve(dof_handler->get_triangulation().n_levels());
2676  for (unsigned int level = 0;
2677  level < dof_handler->get_triangulation().n_levels();
2678  ++level)
2679  {
2680  // first distribute dofs on this level
2681  const types::global_dof_index n_level_dofs =
2683  numbers::invalid_subdomain_id, *dof_handler, level);
2684 
2685  // then add a complete, sequential index set
2686  number_caches.emplace_back(n_level_dofs);
2687  }
2688 
2689  return number_caches;
2690  }
2691 
2692 
2693 
2694  template <int dim, int spacedim>
2695  NumberCache
2697  const std::vector<types::global_dof_index> &new_numbers) const
2698  {
2699  Implementation::renumber_dofs(new_numbers,
2700  IndexSet(0),
2701  *dof_handler,
2702  /*check_validity=*/true);
2703 
2704  // return a sequential, complete index set. take into account that the
2705  // number of DoF indices may in fact be smaller than there were before
2706  // if some previously separately numbered dofs have been identified.
2707  // this is, for example, what we do when the DoFHandler has hp-
2708  // capabilities enabled: it first enumerates all DoFs on cells
2709  // independently, and then unifies some located at vertices or faces;
2710  // this leaves us with fewer DoFs than there were before, so use the
2711  // largest index as the one to determine the size of the index space
2712  if (new_numbers.empty())
2713  return NumberCache();
2714  else
2715  return NumberCache(
2716  *std::max_element(new_numbers.begin(), new_numbers.end()) + 1);
2717  }
2718 
2719 
2720 
2721  template <int dim, int spacedim>
2722  NumberCache
2724  const unsigned int level,
2725  const std::vector<types::global_dof_index> &new_numbers) const
2726  {
2728  new_numbers, IndexSet(0), *dof_handler, level, true);
2729 
2730  // return a sequential, complete index set
2731  return NumberCache(new_numbers.size());
2732  }
2733 
2734 
2735  /* --------------------- class ParallelShared ---------------- */
2736 
2737 
2738  template <int dim, int spacedim>
2740  DoFHandler<dim, spacedim> &dof_handler)
2741  : dof_handler(&dof_handler)
2742  {}
2743 
2744 
2745 
2746  namespace
2747  {
2756  template <int dim, int spacedim>
2757  std::vector<types::subdomain_id>
2758  get_dof_subdomain_association(
2759  const DoFHandler<dim, spacedim> &dof_handler,
2760  const types::global_dof_index n_dofs,
2761  const unsigned int n_procs)
2762  {
2763  (void)n_procs;
2764  std::vector<types::subdomain_id> subdomain_association(
2766  std::vector<types::global_dof_index> local_dof_indices;
2767  local_dof_indices.reserve(
2768  dof_handler.get_fe_collection().max_dofs_per_cell());
2769 
2770  // loop over all cells and record which subdomain a DoF belongs to.
2771  // give to the smaller subdomain_id in case it is on an interface
2772  for (const auto &cell : dof_handler.active_cell_iterators())
2773  {
2774  // get the owner of the cell; note that we have made sure above
2775  // that all cells are either locally owned or ghosts (not
2776  // artificial), so this call will always yield the true owner;
2777  // note that the cache is not assigned yet, so we must bypass it
2778  const types::subdomain_id subdomain_id = cell->subdomain_id();
2779  const unsigned int dofs_per_cell =
2780  cell->get_fe().n_dofs_per_cell();
2781  local_dof_indices.resize(dofs_per_cell);
2782  internal::DoFAccessorImplementation::Implementation::
2783  get_dof_indices(*cell,
2784  local_dof_indices,
2785  cell->active_fe_index());
2786 
2787  // set subdomain ids. if dofs already have their values set then
2788  // they must be on partition interfaces. In that case assign them
2789  // to the processor with the smaller subdomain id.
2790  for (unsigned int i = 0; i < dofs_per_cell; ++i)
2791  if (subdomain_association[local_dof_indices[i]] ==
2793  subdomain_association[local_dof_indices[i]] = subdomain_id;
2794  else if (subdomain_association[local_dof_indices[i]] >
2795  subdomain_id)
2796  {
2797  subdomain_association[local_dof_indices[i]] = subdomain_id;
2798  }
2799  }
2800 
2801  Assert(std::find(subdomain_association.begin(),
2802  subdomain_association.end(),
2804  subdomain_association.end(),
2805  ExcInternalError());
2806 
2807  Assert(*std::max_element(subdomain_association.begin(),
2808  subdomain_association.end()) < n_procs,
2809  ExcInternalError());
2810 
2811  return subdomain_association;
2812  }
2813 
2814 
2821  template <int dim, int spacedim>
2822  std::vector<types::subdomain_id>
2823  get_dof_level_subdomain_association(
2824  const DoFHandler<dim, spacedim> &dof_handler,
2825  const types::global_dof_index n_dofs_on_level,
2826  const unsigned int n_procs,
2827  const unsigned int level)
2828  {
2829  (void)n_procs;
2830  std::vector<types::subdomain_id> level_subdomain_association(
2831  n_dofs_on_level, numbers::invalid_subdomain_id);
2832  std::vector<types::global_dof_index> local_dof_indices;
2833  local_dof_indices.reserve(
2834  dof_handler.get_fe_collection().max_dofs_per_cell());
2835 
2836  // loop over all cells and record which subdomain a DoF belongs to.
2837  // interface goes to processor with smaller subdomain id
2838  for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2839  {
2840  // get the owner of the cell; note that we have made sure above
2841  // that all cells are either locally owned or ghosts (not
2842  // artificial), so this call will always yield the true owner
2843  const types::subdomain_id level_subdomain_id =
2844  cell->level_subdomain_id();
2845  const unsigned int dofs_per_cell =
2846  cell->get_fe().n_dofs_per_cell();
2847  local_dof_indices.resize(dofs_per_cell);
2848  cell->get_mg_dof_indices(local_dof_indices);
2849 
2850  // set level subdomain ids. if dofs already have their values set
2851  // then they must be on partition interfaces. In that case assign
2852  // them to the processor with the smaller subdomain id.
2853  for (unsigned int i = 0; i < dofs_per_cell; ++i)
2854  if (level_subdomain_association[local_dof_indices[i]] ==
2856  level_subdomain_association[local_dof_indices[i]] =
2857  level_subdomain_id;
2858  else if (level_subdomain_association[local_dof_indices[i]] >
2859  level_subdomain_id)
2860  {
2861  level_subdomain_association[local_dof_indices[i]] =
2862  level_subdomain_id;
2863  }
2864  }
2865 
2866  Assert(std::find(level_subdomain_association.begin(),
2867  level_subdomain_association.end(),
2869  level_subdomain_association.end(),
2870  ExcInternalError());
2871 
2872  Assert(*std::max_element(level_subdomain_association.begin(),
2873  level_subdomain_association.end()) < n_procs,
2874  ExcInternalError());
2875 
2876  return level_subdomain_association;
2877  }
2878  } // namespace
2879 
2880 
2881 
2882  template <int dim, int spacedim>
2883  NumberCache
2885  {
2886  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
2887  (dynamic_cast<
2888  const ::parallel::shared::Triangulation<dim, spacedim> *>(
2889  &this->dof_handler->get_triangulation()));
2890  Assert(tr != nullptr, ExcInternalError());
2891 
2892  const unsigned int n_procs =
2893  Utilities::MPI::n_mpi_processes(tr->get_communicator());
2894 
2895  // If an underlying shared::Tria allows artificial cells, we need to
2896  // restore the true cell owners temporarily.
2897  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we
2898  // save the current set of subdomain ids, set subdomain ids to the
2899  // "true" owner of each cell upon construction of the
2900  // TemporarilyRestoreSubdomainIds object, and later restore these flags
2901  // when it is destroyed.
2902  const internal::parallel::shared::
2903  TemporarilyRestoreSubdomainIds<dim, spacedim>
2904  subdomain_modifier(*tr);
2905 
2906  // first let the sequential algorithm do its magic. it is going to
2907  // enumerate DoFs on all cells, regardless of owner
2908  const types::global_dof_index n_initial_dofs =
2910  *this->dof_handler);
2911 
2912  const types::global_dof_index n_dofs =
2913  Implementation::unify_dof_indices(*this->dof_handler,
2914  n_initial_dofs,
2915  /*check_validity=*/true);
2916 
2917  // then re-enumerate them based on their subdomain association.
2918  // for this, we first have to identify for each current DoF
2919  // index which subdomain they belong to. ideally, we would
2920  // like to call DoFRenumbering::subdomain_wise(), but
2921  // because the NumberCache of the current DoFHandler is not
2922  // fully set up yet, we can't quite do that. also, that
2923  // function has to deal with other kinds of triangulations as
2924  // well, whereas we here know what kind of triangulation
2925  // we have and can simplify the code accordingly
2926  std::vector<types::global_dof_index> new_dof_indices(
2927  n_dofs, enumeration_dof_index);
2928  {
2929  // first get the association of each dof with a subdomain and
2930  // determine the total number of subdomain ids used
2931  const std::vector<types::subdomain_id> subdomain_association =
2932  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
2933 
2934  // then renumber the subdomains by first looking at those belonging
2935  // to subdomain 0, then those of subdomain 1, etc. note that the
2936  // algorithm is stable, i.e. if two dofs i,j have i<j and belong to
2937  // the same subdomain, then they will be in this order also after
2938  // reordering
2939  types::global_dof_index next_free_index = 0;
2940  for (types::subdomain_id subdomain = 0; subdomain < n_procs;
2941  ++subdomain)
2942  for (types::global_dof_index i = 0; i < n_dofs; ++i)
2943  if (subdomain_association[i] == subdomain)
2944  {
2945  Assert(new_dof_indices[i] == enumeration_dof_index,
2946  ExcInternalError());
2947  new_dof_indices[i] = next_free_index;
2948  ++next_free_index;
2949  }
2950 
2951  // we should have numbered all dofs
2952  Assert(next_free_index == n_dofs, ExcInternalError());
2953  Assert(std::find(new_dof_indices.begin(),
2954  new_dof_indices.end(),
2955  enumeration_dof_index) == new_dof_indices.end(),
2956  ExcInternalError());
2957  }
2958  // finally do the renumbering. we can use the sequential
2959  // version of the function because we do things on all
2960  // cells and all cells have their subdomain ids and DoFs
2961  // correctly set
2962  Implementation::renumber_dofs(new_dof_indices,
2963  IndexSet(0),
2964  *this->dof_handler,
2965  /*check_validity=*/true);
2966 
2967  // update the number cache. for this, we first have to find the
2968  // subdomain association for each DoF again following renumbering, from
2969  // which we can then compute the IndexSets of locally owned DoFs for all
2970  // processors. all other fields then follow from this
2971  //
2972  // given the way we enumerate degrees of freedom, the locally owned
2973  // ranges must all be contiguous and consecutive. this makes filling
2974  // the IndexSets cheap. an assertion at the top verifies that this
2975  // assumption is true
2976  const std::vector<types::subdomain_id> subdomain_association =
2977  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
2978 
2979  for (unsigned int i = 1; i < n_dofs; ++i)
2980  Assert(subdomain_association[i] >= subdomain_association[i - 1],
2981  ExcInternalError());
2982 
2983  std::vector<IndexSet> locally_owned_dofs_per_processor(
2984  n_procs, IndexSet(n_dofs));
2985  {
2986  // we know that the set of subdomain indices is contiguous from
2987  // the assertion above; find the start and end index for each
2988  // processor, taking into account that sometimes a processor
2989  // may not in fact have any DoFs at all. we do the latter
2990  // by just identifying contiguous ranges of subdomain_ids
2991  // and filling IndexSets for those subdomains; subdomains
2992  // that don't appear will lead to IndexSets that are simply
2993  // never touched and remain empty as initialized above.
2994  unsigned int start_index = 0;
2995  unsigned int end_index = 0;
2996  while (start_index < n_dofs)
2997  {
2998  while ((end_index < n_dofs) &&
2999  (subdomain_association[end_index] ==
3000  subdomain_association[start_index]))
3001  ++end_index;
3002 
3003  // we've now identified a range of same indices. set that
3004  // range in the corresponding IndexSet
3005  if (end_index > start_index)
3006  {
3007  const unsigned int subdomain_owner =
3008  subdomain_association[start_index];
3009  locally_owned_dofs_per_processor[subdomain_owner].add_range(
3010  start_index, end_index);
3011  }
3012 
3013  // then move on to thinking about the next range
3014  start_index = end_index;
3015  }
3016  }
3017 
3018  // return a NumberCache object made up from the sets of locally
3019  // owned DoFs
3020  return NumberCache(
3021  locally_owned_dofs_per_processor,
3022  this->dof_handler->get_triangulation().locally_owned_subdomain());
3023  }
3024 
3025 
3026 
3027  template <int dim, int spacedim>
3028  std::vector<NumberCache>
3030  {
3031  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3032  (dynamic_cast<
3033  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3034  &this->dof_handler->get_triangulation()));
3035  Assert(tr != nullptr, ExcInternalError());
3036 
3037  AssertThrow((tr->is_multilevel_hierarchy_constructed()),
3038  ExcMessage(
3039  "Multigrid DoFs can only be distributed on a parallel "
3040  "Triangulation if the flag construct_multigrid_hierarchy "
3041  "is set in the constructor."));
3042 
3043  const unsigned int n_procs =
3044  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3045  const unsigned int n_levels = tr->n_global_levels();
3046 
3047  std::vector<NumberCache> number_caches;
3048  number_caches.reserve(n_levels);
3049 
3050  // We create an index set for each level
3051  for (unsigned int lvl = 0; lvl < n_levels; ++lvl)
3052  {
3053  // If the underlying shared::Tria allows artificial cells,
3054  // then save the current set of level subdomain ids, and set
3055  // subdomain ids to the "true" owner of each cell. we later
3056  // restore these flags
3057  // Note: "allows_artificial_cells" is currently enforced for
3058  // MG computations.
3059  std::vector<types::subdomain_id> saved_level_subdomain_ids;
3060  saved_level_subdomain_ids.resize(tr->n_cells(lvl));
3061  {
3062  typename ::parallel::shared::Triangulation<dim, spacedim>::
3063  cell_iterator cell =
3064  this->dof_handler->get_triangulation().begin(
3065  lvl),
3066  endc =
3067  this->dof_handler->get_triangulation().end(lvl);
3068 
3069  const std::vector<types::subdomain_id> &true_level_subdomain_ids =
3070  tr->get_true_level_subdomain_ids_of_cells(lvl);
3071 
3072  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3073  {
3074  saved_level_subdomain_ids[index] = cell->level_subdomain_id();
3075  cell->set_level_subdomain_id(true_level_subdomain_ids[index]);
3076  }
3077  }
3078 
3079  // Next let the sequential algorithm do its magic. it is going to
3080  // enumerate DoFs on all cells on the given level, regardless of
3081  // owner
3082  const types::global_dof_index n_dofs_on_level =
3084  numbers::invalid_subdomain_id, *this->dof_handler, lvl);
3085 
3086  // then re-enumerate them based on their level subdomain
3087  // association. for this, we first have to identify for each current
3088  // DoF index which subdomain they belong to. ideally, we would like
3089  // to call DoFRenumbering::subdomain_wise(), but because the
3090  // NumberCache of the current DoFHandler is not fully set up yet, we
3091  // can't quite do that. also, that function has to deal with other
3092  // kinds of triangulations as well, whereas we here know what kind
3093  // of triangulation we have and can simplify the code accordingly
3094  std::vector<types::global_dof_index> new_dof_indices(
3095  n_dofs_on_level, numbers::invalid_dof_index);
3096  {
3097  // first get the association of each dof with a subdomain and
3098  // determine the total number of subdomain ids used
3099  const std::vector<types::subdomain_id>
3100  level_subdomain_association =
3101  get_dof_level_subdomain_association(*this->dof_handler,
3102  n_dofs_on_level,
3103  n_procs,
3104  lvl);
3105 
3106  // then renumber the subdomains by first looking at those
3107  // belonging to subdomain 0, then those of subdomain 1, etc. note
3108  // that the algorithm is stable, i.e. if two dofs i,j have i<j and
3109  // belong to the same subdomain, then they will be in this order
3110  // also after reordering
3111  types::global_dof_index next_free_index = 0;
3112  for (types::subdomain_id level_subdomain = 0;
3113  level_subdomain < n_procs;
3114  ++level_subdomain)
3115  for (types::global_dof_index i = 0; i < n_dofs_on_level; ++i)
3116  if (level_subdomain_association[i] == level_subdomain)
3117  {
3118  Assert(new_dof_indices[i] == numbers::invalid_dof_index,
3119  ExcInternalError());
3120  new_dof_indices[i] = next_free_index;
3121  ++next_free_index;
3122  }
3123 
3124  // we should have numbered all dofs
3125  Assert(next_free_index == n_dofs_on_level, ExcInternalError());
3126  Assert(std::find(new_dof_indices.begin(),
3127  new_dof_indices.end(),
3129  new_dof_indices.end(),
3130  ExcInternalError());
3131  }
3132 
3133  // finally do the renumbering. we can use the sequential
3134  // version of the function because we do things on all
3135  // cells and all cells have their subdomain ids and DoFs
3136  // correctly set
3138  new_dof_indices, IndexSet(0), *this->dof_handler, lvl, true);
3139 
3140  // update the number cache. for this, we first have to find the
3141  // level subdomain association for each DoF again following
3142  // renumbering, from which we can then compute the IndexSets of
3143  // locally owned DoFs for all processors. all other fields then
3144  // follow from this
3145  //
3146  // given the way we enumerate degrees of freedom, the locally owned
3147  // ranges must all be contiguous and consecutive. this makes filling
3148  // the IndexSets cheap. an assertion at the top verifies that this
3149  // assumption is true
3150  const std::vector<types::subdomain_id> level_subdomain_association =
3151  get_dof_level_subdomain_association(*this->dof_handler,
3152  n_dofs_on_level,
3153  n_procs,
3154  lvl);
3155 
3156  for (unsigned int i = 1; i < n_dofs_on_level; ++i)
3157  Assert(level_subdomain_association[i] >=
3158  level_subdomain_association[i - 1],
3159  ExcInternalError());
3160 
3161  std::vector<IndexSet> locally_owned_dofs_per_processor(
3162  n_procs, IndexSet(n_dofs_on_level));
3163  {
3164  // we know that the set of subdomain indices is contiguous from
3165  // the assertion above; find the start and end index for each
3166  // processor, taking into account that sometimes a processor
3167  // may not in fact have any DoFs at all. we do the latter
3168  // by just identifying contiguous ranges of level_subdomain_ids
3169  // and filling IndexSets for those subdomains; subdomains
3170  // that don't appear will lead to IndexSets that are simply
3171  // never touched and remain empty as initialized above.
3172  unsigned int start_index = 0;
3173  unsigned int end_index = 0;
3174  while (start_index < n_dofs_on_level)
3175  {
3176  while ((end_index) < n_dofs_on_level &&
3177  (level_subdomain_association[end_index] ==
3178  level_subdomain_association[start_index]))
3179  ++end_index;
3180 
3181  // we've now identified a range of same indices. set that
3182  // range in the corresponding IndexSet
3183  if (end_index > start_index)
3184  {
3185  const unsigned int level_subdomain_owner =
3186  level_subdomain_association[start_index];
3187  locally_owned_dofs_per_processor[level_subdomain_owner]
3188  .add_range(start_index, end_index);
3189  }
3190 
3191  // then move on to thinking about the next range
3192  start_index = end_index;
3193  }
3194  }
3195 
3196  // finally, restore current level subdomain ids
3197  {
3198  typename ::parallel::shared::Triangulation<dim, spacedim>::
3199  cell_iterator cell =
3200  this->dof_handler->get_triangulation().begin(
3201  lvl),
3202  endc =
3203  this->dof_handler->get_triangulation().end(lvl);
3204 
3205  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3206  cell->set_level_subdomain_id(saved_level_subdomain_ids[index]);
3207 
3208  // add NumberCache for current level
3209  number_caches.emplace_back(
3210  NumberCache(locally_owned_dofs_per_processor,
3211  this->dof_handler->get_triangulation()
3213  }
3214  }
3215 
3216  return number_caches;
3217  }
3218 
3219 
3220 
3221  template <int dim, int spacedim>
3222  NumberCache
3224  const std::vector<types::global_dof_index> &new_numbers) const
3225  {
3226 #ifndef DEAL_II_WITH_MPI
3227  (void)new_numbers;
3228  Assert(false, ExcNotImplemented());
3229  return NumberCache();
3230 #else
3231  // Similar to distribute_dofs() we need to have a special treatment in
3232  // case artificial cells are present.
3233  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3234  (dynamic_cast<
3235  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3236  &this->dof_handler->get_triangulation()));
3237  Assert(tr != nullptr, ExcInternalError());
3238 
3239  // Set subdomain IDs to the "true" owner of each cell.
3240  const internal::parallel::shared::
3241  TemporarilyRestoreSubdomainIds<dim, spacedim>
3242  subdomain_modifier(*tr);
3243 
3244  std::vector<types::global_dof_index> global_gathered_numbers(
3245  this->dof_handler->n_dofs(), 0);
3246  // as we call DoFRenumbering::subdomain_wise(*dof_handler) from
3247  // distribute_dofs(), we need to support sequential-like input.
3248  // Distributed-like input from, for example, component_wise renumbering
3249  // is also supported.
3250  const bool uses_sequential_numbering =
3251  new_numbers.size() == this->dof_handler->n_dofs();
3252  bool all_use_sequential_numbering = false;
3253  Utilities::MPI::internal::all_reduce<bool>(
3254  MPI_LAND,
3255  ArrayView<const bool>(&uses_sequential_numbering, 1),
3256  tr->get_communicator(),
3257  ArrayView<bool>(&all_use_sequential_numbering, 1));
3258  if (all_use_sequential_numbering)
3259  {
3260  global_gathered_numbers = new_numbers;
3261  }
3262  else
3263  {
3264  Assert(new_numbers.size() ==
3265  this->dof_handler->locally_owned_dofs().n_elements(),
3266  ExcInternalError());
3267  const unsigned int n_cpu =
3268  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3269  std::vector<types::global_dof_index> gathered_new_numbers(
3270  this->dof_handler->n_dofs(), 0);
3271  Assert(Utilities::MPI::this_mpi_process(tr->get_communicator()) ==
3272  this->dof_handler->get_triangulation()
3273  .locally_owned_subdomain(),
3274  ExcInternalError());
3275 
3276  // gather new numbers among processors into one vector
3277  {
3278  std::vector<types::global_dof_index> new_numbers_copy(
3279  new_numbers);
3280 
3281  // store the number of elements that are to be received from each
3282  // process
3283  std::vector<int> rcounts(n_cpu);
3284 
3286  // set rcounts based on new_numbers:
3287  int cur_count = new_numbers_copy.size();
3288  int ierr = MPI_Allgather(&cur_count,
3289  1,
3290  MPI_INT,
3291  rcounts.data(),
3292  1,
3293  MPI_INT,
3294  tr->get_communicator());
3295  AssertThrowMPI(ierr);
3296 
3297  // compute the displacements (relative to recvbuf)
3298  // at which to place the incoming data from process i
3299  std::vector<int> displacements(n_cpu);
3300  for (unsigned int i = 0; i < n_cpu; ++i)
3301  {
3302  displacements[i] = shift;
3303  shift += rcounts[i];
3304  }
3305  Assert(new_numbers_copy.size() ==
3306  static_cast<unsigned int>(
3308  tr->get_communicator())]),
3309  ExcInternalError());
3310  ierr = MPI_Allgatherv(new_numbers_copy.data(),
3311  new_numbers_copy.size(),
3313  gathered_new_numbers.data(),
3314  rcounts.data(),
3315  displacements.data(),
3317  tr->get_communicator());
3318  AssertThrowMPI(ierr);
3319  }
3320 
3321  // put new numbers according to the current
3322  // locally_owned_dofs_per_processor IndexSets
3324  // flag_1 and flag_2 are
3325  // used to control that there is a
3326  // one-to-one relation between old and new DoFs.
3327  std::vector<unsigned int> flag_1(this->dof_handler->n_dofs(), 0);
3328  std::vector<unsigned int> flag_2(this->dof_handler->n_dofs(), 0);
3329  std::vector<IndexSet> locally_owned_dofs_per_processor =
3331  tr->get_communicator(),
3332  this->dof_handler->locally_owned_dofs());
3333  for (unsigned int i = 0; i < n_cpu; ++i)
3334  {
3335  const IndexSet iset = locally_owned_dofs_per_processor[i];
3336  for (types::global_dof_index ind = 0; ind < iset.n_elements();
3337  ind++)
3338  {
3339  const types::global_dof_index target =
3340  iset.nth_index_in_set(ind);
3342  gathered_new_numbers[shift + ind];
3343  Assert(target < this->dof_handler->n_dofs(),
3344  ExcInternalError());
3345  Assert(value < this->dof_handler->n_dofs(),
3346  ExcInternalError());
3347  global_gathered_numbers[target] = value;
3348  flag_1[target]++;
3349  flag_2[value]++;
3350  }
3351  shift += iset.n_elements();
3352  }
3353 
3354  Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1,
3355  ExcInternalError());
3356  Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1,
3357  ExcInternalError());
3358  Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1,
3359  ExcInternalError());
3360  Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1,
3361  ExcInternalError());
3362  }
3363 
3364  // let the sequential algorithm do its magic; ignore the
3365  // return type, but reconstruct the number cache based on
3366  // which DoFs each process owns
3367  Implementation::renumber_dofs(global_gathered_numbers,
3368  IndexSet(0),
3369  *this->dof_handler,
3370  /*check_validity=*/true);
3371 
3372  const NumberCache number_cache(
3373  DoFTools::locally_owned_dofs_per_subdomain(*this->dof_handler),
3374  this->dof_handler->get_triangulation().locally_owned_subdomain());
3375 
3376  return number_cache;
3377 #endif
3378  }
3379 
3380 
3381 
3382  template <int dim, int spacedim>
3383  NumberCache
3385  const unsigned int /*level*/,
3386  const std::vector<types::global_dof_index> & /*new_numbers*/) const
3387  {
3388  // multigrid is not currently implemented for shared triangulations
3389  Assert(false, ExcNotImplemented());
3390 
3391  return {};
3392  }
3393 
3394 
3395 
3396  /* --------------------- class ParallelDistributed ---------------- */
3397 
3398 #ifdef DEAL_II_WITH_MPI
3399 
3400  namespace
3401  {
3402  template <int dim, int spacedim>
3403  void
3404  communicate_mg_ghost_cells(DoFHandler<dim, spacedim> &dof_handler,
3405  std::vector<std::vector<bool>> &cell_marked)
3406  {
3407  const auto pack = [](const auto &cell) {
3408  // why would somebody request a cell that is not ours?
3409  Assert(cell->is_locally_owned_on_level(), ExcInternalError());
3410 
3411  std::vector<::types::global_dof_index> data(
3412  cell->get_fe().n_dofs_per_cell());
3413  cell->get_mg_dof_indices(data);
3414 
3415  return data;
3416  };
3417 
3418  const auto unpack = [&cell_marked](const auto &cell,
3419  const auto &dofs) {
3420  Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3421  ExcInternalError());
3422 
3423  Assert(cell->level_subdomain_id() !=
3425  ExcInternalError());
3426 
3427  bool complete = true;
3428  DoFAccessorImplementation::Implementation::process_dof_indices(
3429  *cell,
3430  dofs,
3431  0,
3432  DoFAccessorImplementation::Implementation::
3433  MGDoFIndexProcessor<dim, spacedim>(cell->level()),
3434 
3435  // Intel ICC 18 and earlier for some reason believe that
3436  // numbers::invalid_dof_index is not a valid object
3437  // inside the lambda function. Fix this by creating a
3438  // local variable initialized by the global one.
3439  //
3440  // Intel ICC 19 and earlier have trouble with our Assert
3441  // macros inside the lambda function. We disable the macro
3442  // for these compilers.
3444  auto &stored_index, auto received_index) {
3445  if (*received_index != invalid_dof_index)
3446  {
3447 # if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
3448  Assert((stored_index == invalid_dof_index) ||
3449  (stored_index == *received_index),
3450  ExcInternalError());
3451 # endif
3452  stored_index = *received_index;
3453  }
3454  else
3455  complete = false;
3456  },
3457  true);
3458 
3459  if (!complete)
3460  {
3461  // We should have the cell already marked
3462  Assert(cell_marked[cell->level()][cell->index()],
3463  ExcInternalError());
3464  }
3465  else
3466  cell_marked[cell->level()][cell->index()] = false;
3467  };
3468 
3469  const auto filter = [&cell_marked](const auto &cell) {
3470  return cell_marked[cell->level()][cell->index()];
3471  };
3472 
3474  std::vector<types::global_dof_index>,
3475  DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter);
3476  }
3477 
3478 
3479 
3498  template <int dim, int spacedim>
3499  void
3500  communicate_dof_indices_on_marked_cells(
3501  const DoFHandler<dim, spacedim> &dof_handler,
3502  std::vector<bool> &cell_marked)
3503  {
3504 # ifndef DEAL_II_WITH_MPI
3505  (void)dof_handler;
3506  Assert(false, ExcNotImplemented());
3507 # else
3508 
3509  // define functions that pack data on cells that are ghost cells
3510  // somewhere else, and unpack data on cells where we get information
3511  // from elsewhere
3512  const auto pack = [](const auto &cell) {
3513  Assert(cell->is_locally_owned(), ExcInternalError());
3514 
3515  std::vector<::types::global_dof_index> data(
3516  cell->get_fe().n_dofs_per_cell());
3517 
3518  // bypass the cache which is not filled yet
3519  internal::DoFAccessorImplementation::Implementation::
3520  get_dof_indices(*cell, data, cell->active_fe_index());
3521 
3522  return data;
3523  };
3524 
3525  const auto unpack = [&cell_marked](const auto &cell,
3526  const auto &dofs) {
3527  Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3528  ExcInternalError());
3529 
3530  Assert(cell->is_ghost(), ExcInternalError());
3531 
3532  // Use a combined read/set function on the entities of the dof
3533  // indices to speed things up against get_dof_indices +
3534  // set_dof_indices
3535  bool complete = true;
3536  DoFAccessorImplementation::Implementation::process_dof_indices(
3537  *cell,
3538  dofs,
3539  cell->active_fe_index(),
3540  DoFAccessorImplementation::Implementation::
3541  DoFIndexProcessor<dim, spacedim>(),
3542 
3543  // Intel ICC 18 and earlier for some reason believe that
3544  // numbers::invalid_dof_index is not a valid object
3545  // inside the lambda function. Fix this by creating a
3546  // local variable initialized by the global one.
3547  //
3548  // Intel ICC 19 and earlier have trouble with our Assert
3549  // macros inside the lambda function. We disable the macro
3550  // for these compilers.
3552  auto &stored_index, const auto received_index) {
3553  if (*received_index != invalid_dof_index)
3554  {
3555 # if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
3556  Assert((stored_index == invalid_dof_index) ||
3557  (stored_index == *received_index),
3558  ExcInternalError());
3559 # endif
3560  stored_index = *received_index;
3561  }
3562  else
3563  complete = false;
3564  },
3565  false);
3566 
3567  if (!complete)
3568  {
3569  // We should have the cell already marked
3570  Assert(cell_marked[cell->active_cell_index()],
3571  ExcInternalError());
3572  }
3573  else
3574  cell_marked[cell->active_cell_index()] = false;
3575  };
3576 
3577  const auto filter = [&cell_marked](const auto &cell) {
3578  return cell_marked[cell->active_cell_index()];
3579  };
3580 
3582  std::vector<types::global_dof_index>,
3583  DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter);
3584 # endif
3585  }
3586 
3587 
3588 
3589  } // namespace
3590 
3591 #endif // DEAL_II_WITH_MPI
3592 
3593 
3594 
3595  template <int dim, int spacedim>
3597  DoFHandler<dim, spacedim> &dof_handler)
3598  : dof_handler(&dof_handler)
3599  {}
3600 
3601 
3602 
3603  template <int dim, int spacedim>
3604  NumberCache
3606  {
3607 #ifndef DEAL_II_WITH_MPI
3608  Assert(false, ExcNotImplemented());
3609  return NumberCache();
3610 #else
3611 
3613  *triangulation =
3614  (dynamic_cast<
3616  const_cast<::Triangulation<dim, spacedim> *>(
3617  &dof_handler->get_triangulation())));
3618  Assert(triangulation != nullptr, ExcInternalError());
3619 
3621  triangulation->locally_owned_subdomain();
3622 
3623 
3624  /*
3625  The following algorithm has a number of stages that are all
3626  documented in the paper that describes the parallel::distributed
3627  functionality:
3628 
3629  1/ locally enumerate dofs on locally owned cells
3630  2/ eliminate dof duplicates on all cells.
3631  un-numerate those that are on interfaces with ghost
3632  cells and that we don't own based on the tie-breaking
3633  criterion. unify dofs afterwards.
3634  3/ unify dofs and re-enumerate the remaining valid ones.
3635  the end result is that we only enumerate locally owned
3636  DoFs
3637  4/ shift indices so that each processor has a unique
3638  range of indices
3639  5/ for all locally owned cells that are ghost
3640  cells somewhere else, send our own DoF indices
3641  to the appropriate set of other processors.
3642  overwrite invalid DoF indices on ghost interfaces
3643  with the corresponding valid ones that we now know.
3644  6/ send DoF indices again to get the correct indices
3645  on ghost cells that we may not have known earlier
3646  */
3647 
3648  // --------- Phase 1: enumerate dofs on locally owned cells
3649  const types::global_dof_index n_initial_local_dofs =
3651 
3652  // --------- Phase 2: eliminate dof duplicates on all cells:
3653  // - un-numerate dofs on interfaces to ghost cells
3654  // that we don't own
3655  // - in case of hp-support, unify dofs
3656  std::vector<::types::global_dof_index> renumbering(
3657  n_initial_local_dofs, enumeration_dof_index);
3658 
3659  // first, we invalidate degrees of freedom that belong to processors
3660  // of a lower rank, from which we will receive the final (and lower)
3661  // degrees of freedom later.
3664  renumbering, subdomain_id, *dof_handler);
3665 
3666  // then, we identify DoF duplicates if the DoFHandler has hp-
3667  // capabilities
3668  std::vector<std::map<types::global_dof_index, types::global_dof_index>>
3669  all_constrained_indices(dim);
3670  Implementation::compute_dof_identities(all_constrained_indices,
3671  *dof_handler);
3672 
3673  // --------- Phase 3: re-enumerate the valid degrees of freedom
3674  // consecutively. thus, we finally receive the
3675  // correct number of locally owned DoFs after
3676  // this step.
3677  //
3678  // the order in which we handle Phases 2 and 3 is important,
3679  // since we want to clarify ownership of degrees of freedom before
3680  // we actually unify and enumerate their indices. otherwise, we could
3681  // end up having a degree of freedom to which only invalid indices will
3682  // be assigned.
3683  types::global_dof_index n_identity_constrained_indices = 0;
3684  for (const auto &constrained_indices : all_constrained_indices)
3685  for (const auto &index : constrained_indices)
3686  if (renumbering[index.first] != numbers::invalid_dof_index)
3687  ++n_identity_constrained_indices;
3688 
3689  const types::global_dof_index n_locally_owned_dofs =
3690  std::count(renumbering.begin(),
3691  renumbering.end(),
3692  enumeration_dof_index) -
3693  n_identity_constrained_indices;
3694 
3695  // --------- Phase 4: shift indices so that each processor has a unique
3696  // range of indices
3697  ::types::global_dof_index my_shift = 0;
3698  const int ierr = MPI_Exscan(&n_locally_owned_dofs,
3699  &my_shift,
3700  1,
3702  MPI_SUM,
3703  triangulation->get_communicator());
3704  AssertThrowMPI(ierr);
3705 
3706  // make dof indices globally consecutive
3708  renumbering, all_constrained_indices, my_shift);
3709 
3710  // now re-enumerate all dofs to this shifted and condensed
3711  // numbering form. we renumber some dofs as invalid, so
3712  // choose the nocheck-version.
3713  Implementation::renumber_dofs(renumbering,
3714  IndexSet(0),
3715  *dof_handler,
3716  /*check_validity=*/false);
3717 
3718  // now a little bit of housekeeping
3719  const ::types::global_dof_index n_global_dofs =
3720  Utilities::MPI::sum(n_locally_owned_dofs,
3721  triangulation->get_communicator());
3722 
3723  NumberCache number_cache;
3724  number_cache.n_global_dofs = n_global_dofs;
3725  number_cache.n_locally_owned_dofs = n_locally_owned_dofs;
3726  number_cache.locally_owned_dofs = IndexSet(n_global_dofs);
3727  number_cache.locally_owned_dofs.add_range(my_shift,
3728  my_shift +
3729  n_locally_owned_dofs);
3730  number_cache.locally_owned_dofs.compress();
3731 
3732  // this ends the phase where we enumerate degrees of freedom on
3733  // each processor. what is missing is communicating DoF indices
3734  // on ghost cells
3735 
3736  // --------- Phase 5: for all locally owned cells that are ghost
3737  // cells somewhere else, send our own DoF indices
3738  // to the appropriate set of other processors
3739  {
3740  // mark all cells that either have to send data (locally
3741  // owned cells that are adjacent to ghost neighbors in some
3742  // way) or receive data (all ghost cells) via the user flags
3743  std::vector<bool> cell_marked(triangulation->n_active_cells());
3744  for (const auto &cell : dof_handler->active_cell_iterators())
3745  if (cell->is_ghost())
3746  cell_marked[cell->active_cell_index()] = true;
3747 
3748  // Send and receive cells. After this, only the local cells
3749  // are marked, that received new data. This has to be
3750  // communicated in a second communication step.
3751  //
3752  // as explained in the 'distributed' paper, this has to be
3753  // done twice
3754  communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
3755 
3756  // If the DoFHandler has hp-capabilities enabled, then we may have
3757  // received valid indices of degrees of freedom that are dominated
3758  // by a FE object adjacent to a ghost interface. thus, we overwrite
3759  // the remaining invalid indices with the valid ones in this step.
3761  *dof_handler);
3762 
3763  // --------- Phase 6: all locally owned cells have their correct
3764  // DoF indices set. however, some ghost cells
3765  // may still have invalid ones. thus, exchange
3766  // one more time.
3767  communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
3768 
3769  // at this point, we must have taken care of the data transfer
3770  // on all cells we had previously marked. verify this
3771 # ifdef DEBUG
3772  for (const auto &cell : dof_handler->active_cell_iterators())
3773  Assert(cell_marked[cell->active_cell_index()] == false,
3774  ExcInternalError());
3775 # endif
3776  }
3777 
3778 # ifdef DEBUG
3779  // check that we are really done
3780  {
3781  std::vector<::types::global_dof_index> local_dof_indices;
3782 
3783  for (const auto &cell : dof_handler->active_cell_iterators())
3784  if (!cell->is_artificial())
3785  {
3786  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
3787  cell->get_dof_indices(local_dof_indices);
3788  if (local_dof_indices.end() !=
3789  std::find(local_dof_indices.begin(),
3790  local_dof_indices.end(),
3792  {
3793  if (cell->is_ghost())
3794  {
3795  Assert(false,
3796  ExcMessage(
3797  "A ghost cell ended up with incomplete "
3798  "DoF index information. This should not "
3799  "have happened!"));
3800  }
3801  else
3802  {
3803  Assert(
3804  false,
3805  ExcMessage(
3806  "A locally owned cell ended up with incomplete "
3807  "DoF index information. This should not "
3808  "have happened!"));
3809  }
3810  }
3811  }
3812  }
3813 # endif // DEBUG
3814  return number_cache;
3815 #endif // DEAL_II_WITH_MPI
3816  }
3817 
3818 
3819 
3820  template <int dim, int spacedim>
3821  std::vector<NumberCache>
3823  {
3824 #ifndef DEAL_II_WITH_MPI
3825  Assert(false, ExcNotImplemented());
3826  return std::vector<NumberCache>();
3827 #else
3828 
3830  *triangulation =
3831  (dynamic_cast<
3833  const_cast<::Triangulation<dim, spacedim> *>(
3834  &dof_handler->get_triangulation())));
3835  Assert(triangulation != nullptr, ExcInternalError());
3836 
3837  AssertThrow((triangulation->is_multilevel_hierarchy_constructed()),
3838  ExcMessage(
3839  "Multigrid DoFs can only be distributed on a parallel "
3840  "Triangulation if the flag construct_multigrid_hierarchy "
3841  "is set in the constructor."));
3842 
3843  // loop over all levels that exist globally (across all
3844  // processors), even if the current processor does not in fact
3845  // have any cells on that level or if the local part of the
3846  // Triangulation has fewer levels. we need to do this because
3847  // we need to communicate across all processors on all levels
3848  const unsigned int n_levels = triangulation->n_global_levels();
3849  std::vector<NumberCache> number_caches;
3850  number_caches.reserve(n_levels);
3851  for (unsigned int level = 0; level < n_levels; ++level)
3852  {
3853  NumberCache level_number_cache;
3854 
3855  //* 1. distribute on own subdomain
3856  const unsigned int n_initial_local_dofs =
3858  triangulation->locally_owned_subdomain(), *dof_handler, level);
3859 
3860  //* 2. iterate over ghostcells and kill dofs that are not
3861  // owned by us, which we mark by invalid_dof_index
3862  std::vector<::types::global_dof_index> renumbering(
3863  n_initial_local_dofs, enumeration_dof_index);
3864 
3865  if (level < triangulation->n_levels())
3866  {
3867  std::vector<::types::global_dof_index> local_dof_indices;
3868 
3869  for (const auto &cell :
3870  dof_handler->cell_iterators_on_level(level))
3871  if (cell->level_subdomain_id() !=
3873  (cell->level_subdomain_id() <
3874  triangulation->locally_owned_subdomain()))
3875  {
3876  // we found a neighboring ghost cell whose
3877  // subdomain is "stronger" than our own
3878  // subdomain
3879 
3880  // delete all dofs that live there and that we
3881  // have previously assigned a number to
3882  // (i.e. the ones on the interface)
3883  local_dof_indices.resize(
3884  cell->get_fe().n_dofs_per_cell());
3885  cell->get_mg_dof_indices(local_dof_indices);
3886  for (unsigned int i = 0;
3887  i < cell->get_fe().n_dofs_per_cell();
3888  ++i)
3889  if (local_dof_indices[i] != numbers::invalid_dof_index)
3890  renumbering[local_dof_indices[i]] =
3892  }
3893  }
3894 
3895  level_number_cache.n_locally_owned_dofs =
3896  std::count(renumbering.begin(),
3897  renumbering.end(),
3898  enumeration_dof_index);
3899 
3900  //* 3. communicate local dofcount and shift ids to make
3901  // them unique
3902  ::types::global_dof_index my_shift = 0;
3903  int ierr = MPI_Exscan(&level_number_cache.n_locally_owned_dofs,
3904  &my_shift,
3905  1,
3907  MPI_SUM,
3908  triangulation->get_communicator());
3909  AssertThrowMPI(ierr);
3910 
3911  // The last processor knows about the total number of dofs, so we
3912  // can use a cheaper broadcast rather than an MPI_Allreduce via
3913  // MPI::sum().
3914  level_number_cache.n_global_dofs =
3915  my_shift + level_number_cache.n_locally_owned_dofs;
3916  ierr = MPI_Bcast(&level_number_cache.n_global_dofs,
3917  1,
3920  triangulation->get_communicator()) -
3921  1,
3922  triangulation->get_communicator());
3923  AssertThrowMPI(ierr);
3924 
3925  // assign appropriate indices
3926  for (types::global_dof_index &index : renumbering)
3927  if (index == enumeration_dof_index)
3928  index = my_shift++;
3929 
3930  // now re-enumerate all dofs to this shifted and condensed
3931  // numbering form. we renumber some dofs as invalid, so
3932  // choose the nocheck-version of the function
3933  //
3934  // of course there is nothing for us to renumber if the
3935  // level we are currently dealing with doesn't even exist
3936  // within the current triangulation, so skip renumbering
3937  // in that case
3938  if (level < triangulation->n_levels())
3940  renumbering, IndexSet(0), *dof_handler, level, false);
3941 
3942  // now a little bit of housekeeping
3943  level_number_cache.locally_owned_dofs =
3944  IndexSet(level_number_cache.n_global_dofs);
3945  level_number_cache.locally_owned_dofs.add_range(
3946  my_shift - level_number_cache.n_locally_owned_dofs, my_shift);
3947  level_number_cache.locally_owned_dofs.compress();
3948 
3949  number_caches.emplace_back(level_number_cache);
3950  }
3951 
3952 
3953  //* communicate ghost DoFs
3954  // We mark all ghost cells by setting the user_flag and then request
3955  // these cells from the corresponding owners. As this information
3956  // can be incomplete,
3957  {
3958  std::vector<std::vector<bool>> cell_marked(triangulation->n_levels());
3959  for (unsigned int l = 0; l < triangulation->n_levels(); ++l)
3960  cell_marked[l].resize(triangulation->n_raw_cells(l));
3961  for (const auto &cell : dof_handler->cell_iterators())
3962  if (cell->is_ghost_on_level())
3963  cell_marked[cell->level()][cell->index()] = true;
3964 
3965  // Phase 1. Request all marked cells from corresponding owners. If we
3966  // managed to get every DoF, remove the user_flag, otherwise we
3967  // will request them again in the step below.
3968  communicate_mg_ghost_cells(*dof_handler, cell_marked);
3969 
3970  // Phase 2, only request the cells that were not completed
3971  // in Phase 1.
3972  communicate_mg_ghost_cells(*dof_handler, cell_marked);
3973 
3974 # ifdef DEBUG
3975  // make sure we have finished all cells:
3976  for (const auto &cell : dof_handler->cell_iterators())
3977  Assert(cell_marked[cell->level()][cell->index()] == false,
3978  ExcInternalError());
3979 # endif
3980  }
3981 
3982 
3983 
3984 # ifdef DEBUG
3985  // check that we are really done
3986  {
3987  std::vector<::types::global_dof_index> local_dof_indices;
3988  for (const auto &cell : dof_handler->cell_iterators())
3989  if (cell->level_subdomain_id() !=
3991  {
3992  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
3993  cell->get_mg_dof_indices(local_dof_indices);
3994  if (local_dof_indices.end() !=
3995  std::find(local_dof_indices.begin(),
3996  local_dof_indices.end(),
3998  {
3999  Assert(false, ExcMessage("not all DoFs got distributed!"));
4000  }
4001  }
4002  }
4003 # endif // DEBUG
4004 
4005  return number_caches;
4006 
4007 #endif // DEAL_II_WITH_MPI
4008  }
4009 
4010 
4011  template <int dim, int spacedim>
4012  NumberCache
4014  const std::vector<::types::global_dof_index> &new_numbers) const
4015  {
4016  (void)new_numbers;
4017 
4018  Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
4019  ExcInternalError());
4020 
4021 #ifndef DEAL_II_WITH_MPI
4022  Assert(false, ExcNotImplemented());
4023  return NumberCache();
4024 #else
4025 
4027  *triangulation =
4028  (dynamic_cast<
4030  const_cast<::Triangulation<dim, spacedim> *>(
4031  &dof_handler->get_triangulation())));
4032  Assert(triangulation != nullptr, ExcInternalError());
4033 
4034 
4035  // We start by checking whether only the numbering within the MPI
4036  // ranks changed, in which case we do not need to find a new index
4037  // set.
4038  const IndexSet &owned_dofs = dof_handler->locally_owned_dofs();
4039  const bool locally_owned_set_changes =
4040  std::any_of(new_numbers.cbegin(),
4041  new_numbers.cend(),
4042  [&owned_dofs](const types::global_dof_index i) {
4043  return owned_dofs.is_element(i) == false;
4044  });
4045 
4046  IndexSet my_locally_owned_new_dof_indices = owned_dofs;
4047  if (locally_owned_set_changes && owned_dofs.n_elements() > 0)
4048  {
4049  std::vector<::types::global_dof_index> new_numbers_sorted =
4050  new_numbers;
4051  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4052 
4053  my_locally_owned_new_dof_indices = IndexSet(dof_handler->n_dofs());
4054  my_locally_owned_new_dof_indices.add_indices(
4055  new_numbers_sorted.begin(), new_numbers_sorted.end());
4056  my_locally_owned_new_dof_indices.compress();
4057 
4058  Assert(my_locally_owned_new_dof_indices.n_elements() ==
4059  new_numbers.size(),
4060  ExcInternalError());
4061  }
4062 
4063  // delete all knowledge of DoF indices that are not locally
4064  // owned. we do so by getting DoF indices on cells, checking
4065  // whether they are locally owned, if not, setting them to
4066  // an invalid value, and then setting them again on the current
4067  // cell
4068  //
4069  // DoFs we (i) know about, and (ii) don't own locally must be
4070  // located either on ghost cells, or on the interface between a
4071  // locally owned cell and a ghost cell. In any case, it is
4072  // sufficient to kill them only from the ghost side cell, so loop
4073  // only over ghost cells
4074  for (auto cell : dof_handler->active_cell_iterators())
4075  if (cell->is_ghost())
4076  {
4077  DoFAccessorImplementation::Implementation::process_dof_indices(
4078  *cell,
4079  std::make_tuple(),
4080  cell->active_fe_index(),
4081  DoFAccessorImplementation::Implementation::
4082  DoFIndexProcessor<dim, spacedim>(),
4083  [&owned_dofs](auto &stored_index, auto) {
4084  // delete a DoF index if it has not already been
4085  // deleted (e.g., by visiting a neighboring cell, if
4086  // it is on the boundary), and if we don't own it
4087  if (stored_index != numbers::invalid_dof_index &&
4088  (!owned_dofs.is_element(stored_index)))
4089  stored_index = numbers::invalid_dof_index;
4090  },
4091  false);
4092  }
4093 
4094 
4095  // renumber. Skip when there is nothing to do because we own no DoF.
4096  if (owned_dofs.n_elements() > 0)
4097  Implementation::renumber_dofs(new_numbers,
4098  owned_dofs,
4099  *dof_handler,
4100  /*check_validity=*/false);
4101 
4102  // Communicate newly assigned DoF indices to other processors
4103  // and get the same information for our own ghost cells.
4104  //
4105  // This is the same as phase 5+6 in the distribute_dofs() algorithm,
4106  // taking into account that we have to unify a few DoFs in between
4107  // then communication phases if we do hp-numbering
4108  {
4109  // mark all ghost cells for transfer
4110  std::vector<bool> cell_marked(triangulation->n_active_cells());
4111  for (const auto &cell : dof_handler->active_cell_iterators())
4112  if (cell->is_ghost())
4113  cell_marked[cell->active_cell_index()] = true;
4114 
4115  // Send and receive cells. After this, only the local cells
4116  // are marked, that received new data. This has to be
4117  // communicated in a second communication step.
4118  //
4119  // as explained in the 'distributed' paper, this has to be
4120  // done twice
4121  communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
4122 
4123  // if the DoFHandler has hp-capabilities then we may have
4124  // received valid indices of degrees of freedom that are
4125  // dominated by a FE object adjacent to a ghost interface.
4126  // thus, we overwrite the remaining invalid indices with the
4127  // valid ones in this step.
4129  *dof_handler);
4130 
4131  communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
4132  }
4133 
4134  NumberCache number_cache;
4135  number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
4136  number_cache.n_global_dofs = dof_handler->n_dofs();
4137  number_cache.n_locally_owned_dofs =
4138  number_cache.locally_owned_dofs.n_elements();
4139  return number_cache;
4140 #endif
4141  }
4142 
4143 
4144 
4145  template <int dim, int spacedim>
4146  NumberCache
4148  const unsigned int level,
4149  const std::vector<types::global_dof_index> &new_numbers) const
4150  {
4151 #ifndef DEAL_II_WITH_MPI
4152 
4153  (void)level;
4154  (void)new_numbers;
4155 
4156  Assert(false, ExcNotImplemented());
4157  return NumberCache();
4158 #else
4159 
4161  *triangulation =
4162  (dynamic_cast<
4164  const_cast<::Triangulation<dim, spacedim> *>(
4165  &dof_handler->get_triangulation())));
4166  Assert(triangulation != nullptr, ExcInternalError());
4167 
4168  // This code is very close to the respective code in renumber_dofs,
4169  // with the difference that we work on different entities with
4170  // different objects.
4171  const IndexSet &owned_dofs = dof_handler->locally_owned_mg_dofs(level);
4172  AssertDimension(new_numbers.size(), owned_dofs.n_elements());
4173 
4174  const bool locally_owned_set_changes =
4175  std::any_of(new_numbers.cbegin(),
4176  new_numbers.cend(),
4177  [&owned_dofs](const types::global_dof_index i) {
4178  return owned_dofs.is_element(i) == false;
4179  });
4180 
4181  IndexSet my_locally_owned_new_dof_indices = owned_dofs;
4182  if (locally_owned_set_changes && owned_dofs.n_elements() > 0)
4183  {
4184  std::vector<::types::global_dof_index> new_numbers_sorted =
4185  new_numbers;
4186  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4187 
4188  my_locally_owned_new_dof_indices =
4189  IndexSet(dof_handler->n_dofs(level));
4190  my_locally_owned_new_dof_indices.add_indices(
4191  new_numbers_sorted.begin(), new_numbers_sorted.end());
4192  my_locally_owned_new_dof_indices.compress();
4193 
4194  Assert(my_locally_owned_new_dof_indices.n_elements() ==
4195  new_numbers.size(),
4196  ExcInternalError());
4197  }
4198 
4199  // delete all knowledge of DoF indices that are not locally
4200  // owned
4201  for (auto cell : dof_handler->cell_iterators_on_level(level))
4202  if (cell->is_ghost_on_level())
4203  {
4204  DoFAccessorImplementation::Implementation::process_dof_indices(
4205  *cell,
4206  std::make_tuple(),
4207  0,
4208  DoFAccessorImplementation::Implementation::
4209  MGDoFIndexProcessor<dim, spacedim>(cell->level()),
4210  [&owned_dofs](auto &stored_index, auto) {
4211  if ((stored_index != numbers::invalid_dof_index) &&
4212  (!owned_dofs.is_element(stored_index)))
4213  stored_index = numbers::invalid_dof_index;
4214  },
4215  true);
4216  }
4217 
4218  // renumber. Skip when there is nothing to do because we own no DoF.
4219  if (level < triangulation->n_levels() && owned_dofs.n_elements() > 0)
4221  new_numbers, owned_dofs, *dof_handler, level, false);
4222 
4223  // communicate newly assigned DoF indices with other processors
4224  {
4225  std::vector<std::vector<bool>> cell_marked(triangulation->n_levels());
4226  for (unsigned int l = 0; l < triangulation->n_levels(); ++l)
4227  cell_marked[l].resize(triangulation->n_raw_cells(l));
4228  for (const auto &cell : dof_handler->cell_iterators_on_level(level))
4229  if (cell->is_ghost_on_level())
4230  cell_marked[cell->level()][cell->index()] = true;
4231 
4232  communicate_mg_ghost_cells(*dof_handler, cell_marked);
4233 
4234  communicate_mg_ghost_cells(*dof_handler, cell_marked);
4235  }
4236 
4237  NumberCache number_cache;
4238  number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
4239  number_cache.n_global_dofs = dof_handler->n_dofs(level);
4240  number_cache.n_locally_owned_dofs =
4241  number_cache.locally_owned_dofs.n_elements();
4242  return number_cache;
4243 #endif
4244  }
4245  } // namespace Policy
4246  } // namespace DoFHandlerImplementation
4247 } // namespace internal
4248 
4249 
4250 
4251 /*-------------- Explicit Instantiations -------------------------------*/
4252 #include "dof_handler_policy.inst"
4253 
4254 
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition: dof_handler.h:1597
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1517
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1590
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition: dof_handler.h:1548
const Triangulation< dim, spacedim > & get_triangulation() const
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
Definition: dof_handler.h:1603
bool hp_capability_enabled
Definition: dof_handler.h:1504
types::global_dof_index n_dofs() const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const IndexSet & locally_owned_dofs() const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
types::global_dof_index n_locally_owned_dofs() const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_line() const
unsigned int max_dofs_per_quad() const
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
size_type size() const
Definition: index_set.h:1761
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1986
size_type n_elements() const
Definition: index_set.h:1919
bool is_element(const size_type index) const
Definition: index_set.h:1879
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1788
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1967
void compress() const
Definition: index_set.h:1769
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1816
const std::vector< bool > & get_used_vertices() const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_raw_lines() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_levels() const
cell_iterator end() const
bool vertex_used(const unsigned int index) const
virtual unsigned int n_global_levels() const
unsigned int n_raw_quads() const
unsigned int n_vertices() const
unsigned int size() const
Definition: collection.h:265
unsigned int find_dominating_fe(const std::set< unsigned int > &fes, const unsigned int codim=0) const
unsigned int max_dofs_per_cell() const
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
virtual std::vector< NumberCache > distribute_mg_dofs() const override
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
virtual std::vector< NumberCache > distribute_mg_dofs() const override
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
ParallelShared(DoFHandler< dim, spacedim > &dof_handler)
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
virtual std::vector< NumberCache > distribute_mg_dofs() const override
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 2 > second
Definition: grid_out.cc:4615
Point< 2 > first
Definition: grid_out.cc:4614
unsigned int level
Definition: grid_out.cc:4617
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1947
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
Task< RT > new_task(const std::function< RT()> &function)
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1425
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1412
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:191
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std::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})
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std::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})
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr const ReferenceCell Quadrilateral
std::vector< T > all_gather(const MPI_Comm comm, const T &object_to_send)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition: mpi.cc:149
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition: mpi.cc:164
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1343
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1500
const types::fe_index invalid_fe_index
Definition: types.h:244
const types::subdomain_id artificial_subdomain_id
Definition: types.h:363
const types::subdomain_id invalid_subdomain_id
Definition: types.h:342
static const unsigned int invalid_unsigned_int
Definition: types.h:221
const types::global_dof_index invalid_dof_index
Definition: types.h:253
unsigned int global_dof_index
Definition: types.h:82
unsigned int subdomain_id
Definition: types.h:44
unsigned short int fe_index
Definition: types.h:60
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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)
static void merge_invalid_quad_dofs_on_ghost_interfaces(DoFHandler< 3, spacedim > &dof_handler)
static void renumber_face_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 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)
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const DoFHandler< dim, spacedim > &dof_handler)
static void merge_invalid_line_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
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)
static void merge_invalid_dof_indices_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
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 merge_invalid_vertex_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
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_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 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)
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)
static void merge_invalid_line_dofs_on_ghost_interfaces(DoFHandler< 1, spacedim > &dof_handler)
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)
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)
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)
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)
static types::global_dof_index distribute_dofs(const types::subdomain_id subdomain_id, 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 merge_invalid_quad_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
static void renumber_face_mg_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &, const unsigned int, const bool)
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const DoFHandler< 3, spacedim > &dof_handler)
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)
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const DoFHandler< 1, spacedim > &dof_handler)
static std::map< types::global_dof_index, types::global_dof_index > compute_line_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 > &)
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 types::global_dof_index start_dof_index)
const ::Triangulation< dim, spacedim > & tria
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition: types.h:103