Reference documentation for deal.II version GIT 5983d193e2 2023-05-27 16: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 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
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.size() == 0)
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  [&complete](auto &stored_index, auto received_index) {
3435  if (*received_index != numbers::invalid_dof_index)
3436  {
3437 # if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
3438  Assert((stored_index == (numbers::invalid_dof_index)) ||
3439  (stored_index == *received_index),
3440  ExcInternalError());
3441 # endif
3442  stored_index = *received_index;
3443  }
3444  else
3445  complete = false;
3446  },
3447  true);
3448 
3449  if (!complete)
3450  {
3451  // We should have the cell already marked
3452  Assert(cell_marked[cell->level()][cell->index()],
3453  ExcInternalError());
3454  }
3455  else
3456  cell_marked[cell->level()][cell->index()] = false;
3457  };
3458 
3459  const auto filter = [&cell_marked](const auto &cell) {
3460  return cell_marked[cell->level()][cell->index()];
3461  };
3462 
3464  std::vector<types::global_dof_index>,
3465  DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter);
3466  }
3467 
3468 
3469 
3488  template <int dim, int spacedim>
3489  void
3490  communicate_dof_indices_on_marked_cells(
3491  const DoFHandler<dim, spacedim> &dof_handler,
3492  std::vector<bool> & cell_marked)
3493  {
3494 # ifndef DEAL_II_WITH_MPI
3495  (void)dof_handler;
3496  Assert(false, ExcNotImplemented());
3497 # else
3498 
3499  // define functions that pack data on cells that are ghost cells
3500  // somewhere else, and unpack data on cells where we get information
3501  // from elsewhere
3502  const auto pack = [](const auto &cell) {
3503  Assert(cell->is_locally_owned(), ExcInternalError());
3504 
3505  std::vector<::types::global_dof_index> data(
3506  cell->get_fe().n_dofs_per_cell());
3507 
3508  // bypass the cache which is not filled yet
3509  internal::DoFAccessorImplementation::Implementation::
3510  get_dof_indices(*cell, data, cell->active_fe_index());
3511 
3512  return data;
3513  };
3514 
3515  const auto unpack = [&cell_marked](const auto &cell,
3516  const auto &dofs) {
3517  Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3518  ExcInternalError());
3519 
3520  Assert(cell->is_ghost(), ExcInternalError());
3521 
3522  // Use a combined read/set function on the entities of the dof
3523  // indices to speed things up against get_dof_indices +
3524  // set_dof_indices
3525  bool complete = true;
3526  DoFAccessorImplementation::Implementation::process_dof_indices(
3527  *cell,
3528  dofs,
3529  cell->active_fe_index(),
3530  DoFAccessorImplementation::Implementation::
3531  DoFIndexProcessor<dim, spacedim>(),
3532  [&complete](auto &stored_index, auto received_index) {
3533  if (*received_index != numbers::invalid_dof_index)
3534  {
3535 # if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
3536  Assert((stored_index == (numbers::invalid_dof_index)) ||
3537  (stored_index == *received_index),
3538  ExcInternalError());
3539 # endif
3540  stored_index = *received_index;
3541  }
3542  else
3543  complete = false;
3544  },
3545  false);
3546 
3547  if (!complete)
3548  {
3549  // We should have the cell already marked
3550  Assert(cell_marked[cell->active_cell_index()],
3551  ExcInternalError());
3552  }
3553  else
3554  cell_marked[cell->active_cell_index()] = false;
3555  };
3556 
3557  const auto filter = [&cell_marked](const auto &cell) {
3558  return cell_marked[cell->active_cell_index()];
3559  };
3560 
3562  std::vector<types::global_dof_index>,
3563  DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter);
3564 # endif
3565  }
3566 
3567 
3568 
3569  } // namespace
3570 
3571 #endif // DEAL_II_WITH_MPI
3572 
3573 
3574 
3575  template <int dim, int spacedim>
3577  DoFHandler<dim, spacedim> &dof_handler)
3578  : dof_handler(&dof_handler)
3579  {}
3580 
3581 
3582 
3583  template <int dim, int spacedim>
3584  NumberCache
3586  {
3587 #ifndef DEAL_II_WITH_MPI
3588  Assert(false, ExcNotImplemented());
3589  return NumberCache();
3590 #else
3591 
3593  *triangulation =
3594  (dynamic_cast<
3596  const_cast<::Triangulation<dim, spacedim> *>(
3597  &dof_handler->get_triangulation())));
3598  Assert(triangulation != nullptr, ExcInternalError());
3599 
3601  triangulation->locally_owned_subdomain();
3602 
3603 
3604  /*
3605  The following algorithm has a number of stages that are all
3606  documented in the paper that describes the parallel::distributed
3607  functionality:
3608 
3609  1/ locally enumerate dofs on locally owned cells
3610  2/ eliminate dof duplicates on all cells.
3611  un-numerate those that are on interfaces with ghost
3612  cells and that we don't own based on the tie-breaking
3613  criterion. unify dofs afterwards.
3614  3/ unify dofs and re-enumerate the remaining valid ones.
3615  the end result is that we only enumerate locally owned
3616  DoFs
3617  4/ shift indices so that each processor has a unique
3618  range of indices
3619  5/ for all locally owned cells that are ghost
3620  cells somewhere else, send our own DoF indices
3621  to the appropriate set of other processors.
3622  overwrite invalid DoF indices on ghost interfaces
3623  with the corresponding valid ones that we now know.
3624  6/ send DoF indices again to get the correct indices
3625  on ghost cells that we may not have known earlier
3626  */
3627 
3628  // --------- Phase 1: enumerate dofs on locally owned cells
3629  const types::global_dof_index n_initial_local_dofs =
3631 
3632  // --------- Phase 2: eliminate dof duplicates on all cells:
3633  // - un-numerate dofs on interfaces to ghost cells
3634  // that we don't own
3635  // - in case of hp-support, unify dofs
3636  std::vector<::types::global_dof_index> renumbering(
3637  n_initial_local_dofs, enumeration_dof_index);
3638 
3639  // first, we invalidate degrees of freedom that belong to processors
3640  // of a lower rank, from which we will receive the final (and lower)
3641  // degrees of freedom later.
3644  renumbering, subdomain_id, *dof_handler);
3645 
3646  // then, we identify DoF duplicates if the DoFHandler has hp-
3647  // capabilities
3648  std::vector<std::map<types::global_dof_index, types::global_dof_index>>
3649  all_constrained_indices(dim);
3650  Implementation::compute_dof_identities(all_constrained_indices,
3651  *dof_handler);
3652 
3653  // --------- Phase 3: re-enumerate the valid degrees of freedom
3654  // consecutively. thus, we finally receive the
3655  // correct number of locally owned DoFs after
3656  // this step.
3657  //
3658  // the order in which we handle Phases 2 and 3 is important,
3659  // since we want to clarify ownership of degrees of freedom before
3660  // we actually unify and enumerate their indices. otherwise, we could
3661  // end up having a degree of freedom to which only invalid indices will
3662  // be assigned.
3663  types::global_dof_index n_identity_constrained_indices = 0;
3664  for (const auto &constrained_indices : all_constrained_indices)
3665  for (const auto &index : constrained_indices)
3666  if (renumbering[index.first] != numbers::invalid_dof_index)
3667  ++n_identity_constrained_indices;
3668 
3669  const types::global_dof_index n_locally_owned_dofs =
3670  std::count(renumbering.begin(),
3671  renumbering.end(),
3672  enumeration_dof_index) -
3673  n_identity_constrained_indices;
3674 
3675  // --------- Phase 4: shift indices so that each processor has a unique
3676  // range of indices
3677  ::types::global_dof_index my_shift = 0;
3678  const int ierr = MPI_Exscan(&n_locally_owned_dofs,
3679  &my_shift,
3680  1,
3682  MPI_SUM,
3683  triangulation->get_communicator());
3684  AssertThrowMPI(ierr);
3685 
3686  // make dof indices globally consecutive
3688  renumbering, all_constrained_indices, my_shift);
3689 
3690  // now re-enumerate all dofs to this shifted and condensed
3691  // numbering form. we renumber some dofs as invalid, so
3692  // choose the nocheck-version.
3693  Implementation::renumber_dofs(renumbering,
3694  IndexSet(0),
3695  *dof_handler,
3696  /*check_validity=*/false);
3697 
3698  // now a little bit of housekeeping
3699  const ::types::global_dof_index n_global_dofs =
3700  Utilities::MPI::sum(n_locally_owned_dofs,
3701  triangulation->get_communicator());
3702 
3703  NumberCache number_cache;
3704  number_cache.n_global_dofs = n_global_dofs;
3705  number_cache.n_locally_owned_dofs = n_locally_owned_dofs;
3706  number_cache.locally_owned_dofs = IndexSet(n_global_dofs);
3707  number_cache.locally_owned_dofs.add_range(my_shift,
3708  my_shift +
3709  n_locally_owned_dofs);
3710  number_cache.locally_owned_dofs.compress();
3711 
3712  // this ends the phase where we enumerate degrees of freedom on
3713  // each processor. what is missing is communicating DoF indices
3714  // on ghost cells
3715 
3716  // --------- Phase 5: for all locally owned cells that are ghost
3717  // cells somewhere else, send our own DoF indices
3718  // to the appropriate set of other processors
3719  {
3720  // mark all cells that either have to send data (locally
3721  // owned cells that are adjacent to ghost neighbors in some
3722  // way) or receive data (all ghost cells) via the user flags
3723  std::vector<bool> cell_marked(triangulation->n_active_cells());
3724  for (const auto &cell : dof_handler->active_cell_iterators())
3725  if (cell->is_ghost())
3726  cell_marked[cell->active_cell_index()] = true;
3727 
3728  // Send and receive cells. After this, only the local cells
3729  // are marked, that received new data. This has to be
3730  // communicated in a second communication step.
3731  //
3732  // as explained in the 'distributed' paper, this has to be
3733  // done twice
3734  communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
3735 
3736  // If the DoFHandler has hp-capabilities enabled, then we may have
3737  // received valid indices of degrees of freedom that are dominated
3738  // by a FE object adjacent to a ghost interface. thus, we overwrite
3739  // the remaining invalid indices with the valid ones in this step.
3741  *dof_handler);
3742 
3743  // --------- Phase 6: all locally owned cells have their correct
3744  // DoF indices set. however, some ghost cells
3745  // may still have invalid ones. thus, exchange
3746  // one more time.
3747  communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
3748 
3749  // at this point, we must have taken care of the data transfer
3750  // on all cells we had previously marked. verify this
3751 # ifdef DEBUG
3752  for (const auto &cell : dof_handler->active_cell_iterators())
3753  Assert(cell_marked[cell->active_cell_index()] == false,
3754  ExcInternalError());
3755 # endif
3756  }
3757 
3758 # ifdef DEBUG
3759  // check that we are really done
3760  {
3761  std::vector<::types::global_dof_index> local_dof_indices;
3762 
3763  for (const auto &cell : dof_handler->active_cell_iterators())
3764  if (!cell->is_artificial())
3765  {
3766  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
3767  cell->get_dof_indices(local_dof_indices);
3768  if (local_dof_indices.end() !=
3769  std::find(local_dof_indices.begin(),
3770  local_dof_indices.end(),
3772  {
3773  if (cell->is_ghost())
3774  {
3775  Assert(false,
3776  ExcMessage(
3777  "A ghost cell ended up with incomplete "
3778  "DoF index information. This should not "
3779  "have happened!"));
3780  }
3781  else
3782  {
3783  Assert(
3784  false,
3785  ExcMessage(
3786  "A locally owned cell ended up with incomplete "
3787  "DoF index information. This should not "
3788  "have happened!"));
3789  }
3790  }
3791  }
3792  }
3793 # endif // DEBUG
3794  return number_cache;
3795 #endif // DEAL_II_WITH_MPI
3796  }
3797 
3798 
3799 
3800  template <int dim, int spacedim>
3801  std::vector<NumberCache>
3803  {
3804 #ifndef DEAL_II_WITH_MPI
3805  Assert(false, ExcNotImplemented());
3806  return std::vector<NumberCache>();
3807 #else
3808 
3810  *triangulation =
3811  (dynamic_cast<
3813  const_cast<::Triangulation<dim, spacedim> *>(
3814  &dof_handler->get_triangulation())));
3815  Assert(triangulation != nullptr, ExcInternalError());
3816 
3817  AssertThrow((triangulation->is_multilevel_hierarchy_constructed()),
3818  ExcMessage(
3819  "Multigrid DoFs can only be distributed on a parallel "
3820  "Triangulation if the flag construct_multigrid_hierarchy "
3821  "is set in the constructor."));
3822 
3823  // loop over all levels that exist globally (across all
3824  // processors), even if the current processor does not in fact
3825  // have any cells on that level or if the local part of the
3826  // Triangulation has fewer levels. we need to do this because
3827  // we need to communicate across all processors on all levels
3828  const unsigned int n_levels = triangulation->n_global_levels();
3829  std::vector<NumberCache> number_caches;
3830  number_caches.reserve(n_levels);
3831  for (unsigned int level = 0; level < n_levels; ++level)
3832  {
3833  NumberCache level_number_cache;
3834 
3835  //* 1. distribute on own subdomain
3836  const unsigned int n_initial_local_dofs =
3838  triangulation->locally_owned_subdomain(), *dof_handler, level);
3839 
3840  //* 2. iterate over ghostcells and kill dofs that are not
3841  // owned by us, which we mark by invalid_dof_index
3842  std::vector<::types::global_dof_index> renumbering(
3843  n_initial_local_dofs, enumeration_dof_index);
3844 
3845  if (level < triangulation->n_levels())
3846  {
3847  std::vector<::types::global_dof_index> local_dof_indices;
3848 
3849  for (const auto &cell :
3850  dof_handler->cell_iterators_on_level(level))
3851  if (cell->level_subdomain_id() !=
3853  (cell->level_subdomain_id() <
3854  triangulation->locally_owned_subdomain()))
3855  {
3856  // we found a neighboring ghost cell whose
3857  // subdomain is "stronger" than our own
3858  // subdomain
3859 
3860  // delete all dofs that live there and that we
3861  // have previously assigned a number to
3862  // (i.e. the ones on the interface)
3863  local_dof_indices.resize(
3864  cell->get_fe().n_dofs_per_cell());
3865  cell->get_mg_dof_indices(local_dof_indices);
3866  for (unsigned int i = 0;
3867  i < cell->get_fe().n_dofs_per_cell();
3868  ++i)
3869  if (local_dof_indices[i] != numbers::invalid_dof_index)
3870  renumbering[local_dof_indices[i]] =
3872  }
3873  }
3874 
3875  level_number_cache.n_locally_owned_dofs =
3876  std::count(renumbering.begin(),
3877  renumbering.end(),
3878  enumeration_dof_index);
3879 
3880  //* 3. communicate local dofcount and shift ids to make
3881  // them unique
3882  ::types::global_dof_index my_shift = 0;
3883  int ierr = MPI_Exscan(&level_number_cache.n_locally_owned_dofs,
3884  &my_shift,
3885  1,
3887  MPI_SUM,
3888  triangulation->get_communicator());
3889  AssertThrowMPI(ierr);
3890 
3891  // The last processor knows about the total number of dofs, so we
3892  // can use a cheaper broadcast rather than an MPI_Allreduce via
3893  // MPI::sum().
3894  level_number_cache.n_global_dofs =
3895  my_shift + level_number_cache.n_locally_owned_dofs;
3896  ierr = MPI_Bcast(&level_number_cache.n_global_dofs,
3897  1,
3900  triangulation->get_communicator()) -
3901  1,
3902  triangulation->get_communicator());
3903  AssertThrowMPI(ierr);
3904 
3905  // assign appropriate indices
3906  for (types::global_dof_index &index : renumbering)
3907  if (index == enumeration_dof_index)
3908  index = my_shift++;
3909 
3910  // now re-enumerate all dofs to this shifted and condensed
3911  // numbering form. we renumber some dofs as invalid, so
3912  // choose the nocheck-version of the function
3913  //
3914  // of course there is nothing for us to renumber if the
3915  // level we are currently dealing with doesn't even exist
3916  // within the current triangulation, so skip renumbering
3917  // in that case
3918  if (level < triangulation->n_levels())
3920  renumbering, IndexSet(0), *dof_handler, level, false);
3921 
3922  // now a little bit of housekeeping
3923  level_number_cache.locally_owned_dofs =
3924  IndexSet(level_number_cache.n_global_dofs);
3925  level_number_cache.locally_owned_dofs.add_range(
3926  my_shift - level_number_cache.n_locally_owned_dofs, my_shift);
3927  level_number_cache.locally_owned_dofs.compress();
3928 
3929  number_caches.emplace_back(level_number_cache);
3930  }
3931 
3932 
3933  //* communicate ghost DoFs
3934  // We mark all ghost cells by setting the user_flag and then request
3935  // these cells from the corresponding owners. As this information
3936  // can be incomplete,
3937  {
3938  std::vector<std::vector<bool>> cell_marked(triangulation->n_levels());
3939  for (unsigned int l = 0; l < triangulation->n_levels(); ++l)
3940  cell_marked[l].resize(triangulation->n_raw_cells(l));
3941  for (const auto &cell : dof_handler->cell_iterators())
3942  if (cell->is_ghost_on_level())
3943  cell_marked[cell->level()][cell->index()] = true;
3944 
3945  // Phase 1. Request all marked cells from corresponding owners. If we
3946  // managed to get every DoF, remove the user_flag, otherwise we
3947  // will request them again in the step below.
3948  communicate_mg_ghost_cells(*dof_handler, cell_marked);
3949 
3950  // Phase 2, only request the cells that were not completed
3951  // in Phase 1.
3952  communicate_mg_ghost_cells(*dof_handler, cell_marked);
3953 
3954 # ifdef DEBUG
3955  // make sure we have finished all cells:
3956  for (const auto &cell : dof_handler->cell_iterators())
3957  Assert(cell_marked[cell->level()][cell->index()] == false,
3958  ExcInternalError());
3959 # endif
3960  }
3961 
3962 
3963 
3964 # ifdef DEBUG
3965  // check that we are really done
3966  {
3967  std::vector<::types::global_dof_index> local_dof_indices;
3968  for (const auto &cell : dof_handler->cell_iterators())
3969  if (cell->level_subdomain_id() !=
3971  {
3972  local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
3973  cell->get_mg_dof_indices(local_dof_indices);
3974  if (local_dof_indices.end() !=
3975  std::find(local_dof_indices.begin(),
3976  local_dof_indices.end(),
3978  {
3979  Assert(false, ExcMessage("not all DoFs got distributed!"));
3980  }
3981  }
3982  }
3983 # endif // DEBUG
3984 
3985  return number_caches;
3986 
3987 #endif // DEAL_II_WITH_MPI
3988  }
3989 
3990 
3991  template <int dim, int spacedim>
3992  NumberCache
3994  const std::vector<::types::global_dof_index> &new_numbers) const
3995  {
3996  (void)new_numbers;
3997 
3998  Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
3999  ExcInternalError());
4000 
4001 #ifndef DEAL_II_WITH_MPI
4002  Assert(false, ExcNotImplemented());
4003  return NumberCache();
4004 #else
4005 
4007  *triangulation =
4008  (dynamic_cast<
4010  const_cast<::Triangulation<dim, spacedim> *>(
4011  &dof_handler->get_triangulation())));
4012  Assert(triangulation != nullptr, ExcInternalError());
4013 
4014 
4015  // We start by checking whether only the numbering within the MPI
4016  // ranks changed, in which case we do not need to find a new index
4017  // set.
4018  const IndexSet &owned_dofs = dof_handler->locally_owned_dofs();
4019  const bool locally_owned_set_changes =
4020  std::any_of(new_numbers.cbegin(),
4021  new_numbers.cend(),
4022  [&owned_dofs](const types::global_dof_index i) {
4023  return owned_dofs.is_element(i) == false;
4024  });
4025 
4026  IndexSet my_locally_owned_new_dof_indices = owned_dofs;
4027  if (locally_owned_set_changes && owned_dofs.n_elements() > 0)
4028  {
4029  std::vector<::types::global_dof_index> new_numbers_sorted =
4030  new_numbers;
4031  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4032 
4033  my_locally_owned_new_dof_indices = IndexSet(dof_handler->n_dofs());
4034  my_locally_owned_new_dof_indices.add_indices(
4035  new_numbers_sorted.begin(), new_numbers_sorted.end());
4036  my_locally_owned_new_dof_indices.compress();
4037 
4038  Assert(my_locally_owned_new_dof_indices.n_elements() ==
4039  new_numbers.size(),
4040  ExcInternalError());
4041  }
4042 
4043  // delete all knowledge of DoF indices that are not locally
4044  // owned. we do so by getting DoF indices on cells, checking
4045  // whether they are locally owned, if not, setting them to
4046  // an invalid value, and then setting them again on the current
4047  // cell
4048  //
4049  // DoFs we (i) know about, and (ii) don't own locally must be
4050  // located either on ghost cells, or on the interface between a
4051  // locally owned cell and a ghost cell. In any case, it is
4052  // sufficient to kill them only from the ghost side cell, so loop
4053  // only over ghost cells
4054  for (auto cell : dof_handler->active_cell_iterators())
4055  if (cell->is_ghost())
4056  {
4057  DoFAccessorImplementation::Implementation::process_dof_indices(
4058  *cell,
4059  std::make_tuple(),
4060  cell->active_fe_index(),
4061  DoFAccessorImplementation::Implementation::
4062  DoFIndexProcessor<dim, spacedim>(),
4063  [&owned_dofs](auto &stored_index, auto) {
4064  // delete a DoF index if it has not already been
4065  // deleted (e.g., by visiting a neighboring cell, if
4066  // it is on the boundary), and if we don't own it
4067  if (stored_index != numbers::invalid_dof_index &&
4068  (!owned_dofs.is_element(stored_index)))
4069  stored_index = numbers::invalid_dof_index;
4070  },
4071  false);
4072  }
4073 
4074 
4075  // renumber. Skip when there is nothing to do because we own no DoF.
4076  if (owned_dofs.n_elements() > 0)
4077  Implementation::renumber_dofs(new_numbers,
4078  owned_dofs,
4079  *dof_handler,
4080  /*check_validity=*/false);
4081 
4082  // Communicate newly assigned DoF indices to other processors
4083  // and get the same information for our own ghost cells.
4084  //
4085  // This is the same as phase 5+6 in the distribute_dofs() algorithm,
4086  // taking into account that we have to unify a few DoFs in between
4087  // then communication phases if we do hp-numbering
4088  {
4089  // mark all ghost cells for transfer
4090  std::vector<bool> cell_marked(triangulation->n_active_cells());
4091  for (const auto &cell : dof_handler->active_cell_iterators())
4092  if (cell->is_ghost())
4093  cell_marked[cell->active_cell_index()] = true;
4094 
4095  // Send and receive cells. After this, only the local cells
4096  // are marked, that received new data. This has to be
4097  // communicated in a second communication step.
4098  //
4099  // as explained in the 'distributed' paper, this has to be
4100  // done twice
4101  communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
4102 
4103  // if the DoFHandler has hp-capabilities then we may have
4104  // received valid indices of degrees of freedom that are
4105  // dominated by a FE object adjacent to a ghost interface.
4106  // thus, we overwrite the remaining invalid indices with the
4107  // valid ones in this step.
4109  *dof_handler);
4110 
4111  communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
4112  }
4113 
4114  NumberCache number_cache;
4115  number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
4116  number_cache.n_global_dofs = dof_handler->n_dofs();
4117  number_cache.n_locally_owned_dofs =
4118  number_cache.locally_owned_dofs.n_elements();
4119  return number_cache;
4120 #endif
4121  }
4122 
4123 
4124 
4125  template <int dim, int spacedim>
4126  NumberCache
4128  const unsigned int level,
4129  const std::vector<types::global_dof_index> &new_numbers) const
4130  {
4131 #ifndef DEAL_II_WITH_MPI
4132 
4133  (void)level;
4134  (void)new_numbers;
4135 
4136  Assert(false, ExcNotImplemented());
4137  return NumberCache();
4138 #else
4139 
4141  *triangulation =
4142  (dynamic_cast<
4144  const_cast<::Triangulation<dim, spacedim> *>(
4145  &dof_handler->get_triangulation())));
4146  Assert(triangulation != nullptr, ExcInternalError());
4147 
4148  // This code is very close to the respective code in renumber_dofs,
4149  // with the difference that we work on different entities with
4150  // different objects.
4151  const IndexSet &owned_dofs = dof_handler->locally_owned_mg_dofs(level);
4152  AssertDimension(new_numbers.size(), owned_dofs.n_elements());
4153 
4154  const bool locally_owned_set_changes =
4155  std::any_of(new_numbers.cbegin(),
4156  new_numbers.cend(),
4157  [&owned_dofs](const types::global_dof_index i) {
4158  return owned_dofs.is_element(i) == false;
4159  });
4160 
4161  IndexSet my_locally_owned_new_dof_indices = owned_dofs;
4162  if (locally_owned_set_changes && owned_dofs.n_elements() > 0)
4163  {
4164  std::vector<::types::global_dof_index> new_numbers_sorted =
4165  new_numbers;
4166  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4167 
4168  my_locally_owned_new_dof_indices =
4169  IndexSet(dof_handler->n_dofs(level));
4170  my_locally_owned_new_dof_indices.add_indices(
4171  new_numbers_sorted.begin(), new_numbers_sorted.end());
4172  my_locally_owned_new_dof_indices.compress();
4173 
4174  Assert(my_locally_owned_new_dof_indices.n_elements() ==
4175  new_numbers.size(),
4176  ExcInternalError());
4177  }
4178 
4179  // delete all knowledge of DoF indices that are not locally
4180  // owned
4181  for (auto cell : dof_handler->cell_iterators_on_level(level))
4182  if (cell->is_ghost_on_level())
4183  {
4184  DoFAccessorImplementation::Implementation::process_dof_indices(
4185  *cell,
4186  std::make_tuple(),
4187  0,
4188  DoFAccessorImplementation::Implementation::
4189  MGDoFIndexProcessor<dim, spacedim>(cell->level()),
4190  [&owned_dofs](auto &stored_index, auto) {
4191  if ((stored_index != numbers::invalid_dof_index) &&
4192  (!owned_dofs.is_element(stored_index)))
4193  stored_index = numbers::invalid_dof_index;
4194  },
4195  true);
4196  }
4197 
4198  // renumber. Skip when there is nothing to do because we own no DoF.
4199  if (level < triangulation->n_levels() && owned_dofs.n_elements() > 0)
4201  new_numbers, owned_dofs, *dof_handler, level, false);
4202 
4203  // communicate newly assigned DoF indices with other processors
4204  {
4205  std::vector<std::vector<bool>> cell_marked(triangulation->n_levels());
4206  for (unsigned int l = 0; l < triangulation->n_levels(); ++l)
4207  cell_marked[l].resize(triangulation->n_raw_cells(l));
4208  for (const auto &cell : dof_handler->cell_iterators_on_level(level))
4209  if (cell->is_ghost_on_level())
4210  cell_marked[cell->level()][cell->index()] = true;
4211 
4212  communicate_mg_ghost_cells(*dof_handler, cell_marked);
4213 
4214  communicate_mg_ghost_cells(*dof_handler, cell_marked);
4215  }
4216 
4217  NumberCache number_cache;
4218  number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
4219  number_cache.n_global_dofs = dof_handler->n_dofs(level);
4220  number_cache.n_locally_owned_dofs =
4221  number_cache.locally_owned_dofs.n_elements();
4222  return number_cache;
4223 #endif
4224  }
4225  } // namespace Policy
4226  } // namespace DoFHandlerImplementation
4227 } // namespace internal
4228 
4229 
4230 
4231 /*-------------- Explicit Instantiations -------------------------------*/
4232 #include "dof_handler_policy.inst"
4233 
4234 
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition: dof_handler.h:1599
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1519
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1592
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition: dof_handler.h:1550
const Triangulation< dim, spacedim > & get_triangulation() const
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
Definition: dof_handler.h:1605
bool hp_capability_enabled
Definition: dof_handler.h:1506
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:1661
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1883
size_type n_elements() const
Definition: index_set.h:1816
bool is_element(const size_type index) const
Definition: index_set.h:1776
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1697
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1864
void compress() const
Definition: index_set.h:1669
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1727
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:264
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:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
Point< 2 > second
Definition: grid_out.cc:4616
Point< 2 > first
Definition: grid_out.cc:4615
unsigned int level
Definition: grid_out.cc:4618
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:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1914
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
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:1376
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1363
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2131
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
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:150
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition: mpi.cc:161
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1351
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1510
const types::fe_index invalid_fe_index
Definition: types.h:228
const types::subdomain_id artificial_subdomain_id
Definition: types.h:315
const types::subdomain_id invalid_subdomain_id
Definition: types.h:298
static const unsigned int invalid_unsigned_int
Definition: types.h:213
const types::global_dof_index invalid_dof_index
Definition: types.h:233
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:99