Reference documentation for deal.II version Git 92c1895bf8 2020-06-04 17:48:28 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_handler_policy.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
21 #include <deal.II/base/utilities.h>
23 
26 
30 
31 #include <deal.II/fe/fe.h>
32 
34 #include <deal.II/grid/tria.h>
36 
37 #include <algorithm>
38 #include <memory>
39 #include <numeric>
40 #include <set>
41 
43 
44 
45 namespace internal
46 {
47  namespace DoFHandlerImplementation
48  {
49  namespace Policy
50  {
51  // use class ::DoFHandler instead
52  // of namespace internal::DoFHandler in
53  // the following
54  using ::DoFHandler;
55 
56  namespace hp
57  {
58  using ::hp::DoFHandler;
59  }
60 
61 
62  namespace
63  {
70  const types::global_dof_index enumeration_dof_index =
72 
77  template <class DoFHandlerType>
78  void
79  update_all_active_cell_dof_indices_caches(
80  const DoFHandlerType &dof_handler)
81  {
82  typename DoFHandlerType::active_cell_iterator
83  beginc = dof_handler.begin_active(),
84  endc = dof_handler.end();
85 
86  auto worker =
87  [](const typename DoFHandlerType::active_cell_iterator &cell,
88  void *,
89  void *) {
90  if (!cell->is_artificial())
91  cell->update_cell_dof_indices_cache();
92  };
93 
94  // parallelize filling all of the cell caches. by using
95  // WorkStream, we make sure that we only run through the
96  // range of iterators once, whereas a parallel_for loop
97  // for example has to split the range multiple times,
98  // which is expensive because cell iterators are not
99  // random access iterators with a cheap operator-
100  WorkStream::run(beginc,
101  endc,
102  worker,
103  /* copier */ std::function<void(void *)>(),
104  /* scratch_data */ nullptr,
105  /* copy_data */ nullptr,
107  /* chunk_size = */ 32);
108  }
109 
110 
115  template <class DoFHandlerType>
116  void
117  update_all_level_cell_dof_indices_caches(
118  const DoFHandlerType &dof_handler)
119  {
120  typename DoFHandlerType::level_cell_iterator beginc =
121  dof_handler.begin(),
122  endc = dof_handler.end();
123 
124  auto worker =
125  [](const typename DoFHandlerType::level_cell_iterator &cell,
126  void *,
127  void *) {
128  if (cell->has_children() || !cell->is_artificial())
129  cell->update_cell_dof_indices_cache();
130  };
131 
132  // parallelize filling all of the cell caches. by using
133  // WorkStream, we make sure that we only run through the
134  // range of iterators once, whereas a parallel_for loop
135  // for example has to split the range multiple times,
136  // which is expensive because cell iterators are not
137  // random access iterators with a cheap operator-
138  WorkStream::run(beginc,
139  endc,
140  worker,
141  /* copier */ std::function<void(void *)>(),
142  /* scratch_data */ nullptr,
143  /* copy_data */ nullptr,
145  /* chunk_size = */ 32);
146  }
147 
148 
149  using DoFIdentities =
150  std::vector<std::pair<unsigned int, unsigned int>>;
151 
152 
163  template <int structdim, int dim, int spacedim>
164  void
165  ensure_existence_of_dof_identities(
166  const FiniteElement<dim, spacedim> &fe1,
167  const FiniteElement<dim, spacedim> &fe2,
168  std::unique_ptr<DoFIdentities> & identities)
169  {
170  // see if we need to fill this entry, or whether it already
171  // exists
172  if (identities.get() == nullptr)
173  {
174  switch (structdim)
175  {
176  case 0:
177  {
178  identities = std::make_unique<DoFIdentities>(
179  fe1.hp_vertex_dof_identities(fe2));
180  break;
181  }
182 
183  case 1:
184  {
185  identities = std::make_unique<DoFIdentities>(
186  fe1.hp_line_dof_identities(fe2));
187  break;
188  }
189 
190  case 2:
191  {
192  identities = std::make_unique<DoFIdentities>(
193  fe1.hp_quad_dof_identities(fe2));
194  break;
195  }
196 
197  default:
198  Assert(false, ExcNotImplemented());
199  }
200 
201  // double check whether the newly created entries make
202  // any sense at all
203  for (unsigned int i = 0; i < identities->size(); ++i)
204  {
205  Assert((*identities)[i].first <
206  fe1.template n_dofs_per_object<structdim>(),
207  ExcInternalError());
208  Assert((*identities)[i].second <
209  fe2.template n_dofs_per_object<structdim>(),
210  ExcInternalError());
211  }
212  }
213  }
214  } // namespace
215 
216 
217 
219  {
220  /* -------------- distribute_dofs functionality ------------- */
221 
226  template <int dim, int spacedim>
227  static std::map<types::global_dof_index, types::global_dof_index>
229  const hp::DoFHandler<dim, spacedim> &dof_handler)
230  {
231  std::map<types::global_dof_index, types::global_dof_index>
232  dof_identities;
233 
234  // Note: we may wish to have something here similar to what
235  // we do for lines and quads, namely that we only identify
236  // dofs for any fe towards the most dominating one. however,
237  // it is not clear whether this is actually necessary for
238  // vertices at all, I can't think of a finite element that
239  // would make that necessary...
241  vertex_dof_identities(dof_handler.get_fe_collection().size(),
242  dof_handler.get_fe_collection().size());
243 
244  // loop over all vertices and see which one we need to work on
245  for (unsigned int vertex_index = 0;
246  vertex_index < dof_handler.get_triangulation().n_vertices();
247  ++vertex_index)
248  if (dof_handler.get_triangulation()
249  .get_used_vertices()[vertex_index] == true)
250  {
251  const unsigned int n_active_fe_indices =
252  ::internal::DoFAccessorImplementation::Implementation::
253  n_active_vertex_fe_indices(dof_handler, vertex_index);
254 
255  if (n_active_fe_indices > 1)
256  {
257  const std::set<unsigned int> fe_indices =
258  ::internal::DoFAccessorImplementation::
259  Implementation::get_active_vertex_fe_indices(
260  dof_handler, vertex_index);
261 
262  // find out which is the most dominating finite
263  // element of the ones that are used on this vertex
264  unsigned int most_dominating_fe_index =
265  dof_handler.get_fe_collection().find_dominating_fe(
266  fe_indices,
267  /*codim*/ dim);
268 
269  // if we haven't found a dominating finite element,
270  // choose the very first one to be dominant
271  if (most_dominating_fe_index ==
273  most_dominating_fe_index =
274  ::internal::DoFAccessorImplementation::
275  Implementation::nth_active_vertex_fe_index(
276  dof_handler, vertex_index, 0);
277 
278  // loop over the indices of all the finite
279  // elements that are not dominating, and
280  // identify their dofs to the most dominating
281  // one
282  for (const auto &other_fe_index : fe_indices)
283  if (other_fe_index != most_dominating_fe_index)
284  {
285  // make sure the entry in the equivalence
286  // table exists
287  ensure_existence_of_dof_identities<0>(
288  dof_handler.get_fe(most_dominating_fe_index),
289  dof_handler.get_fe(other_fe_index),
290  vertex_dof_identities[most_dominating_fe_index]
291  [other_fe_index]);
292 
293  // then loop through the identities we
294  // have. first get the global numbers of the
295  // dofs we want to identify and make sure they
296  // are not yet constrained to anything else,
297  // except for to each other. use the rule that
298  // we will always constrain the dof with the
299  // higher fe index to the one with the lower,
300  // to avoid circular reasoning.
301  DoFIdentities &identities =
302  *vertex_dof_identities[most_dominating_fe_index]
303  [other_fe_index];
304  for (const auto &identity : identities)
305  {
306  const types::global_dof_index master_dof_index =
307  ::internal::DoFAccessorImplementation::
308  Implementation::get_vertex_dof_index(
309  dof_handler,
310  vertex_index,
311  most_dominating_fe_index,
312  identity.first);
313  const types::global_dof_index slave_dof_index =
314  ::internal::DoFAccessorImplementation::
315  Implementation::get_vertex_dof_index(
316  dof_handler,
317  vertex_index,
318  other_fe_index,
319  identity.second);
320 
321  // on subdomain boundaries, we will
322  // encounter invalid DoFs on ghost cells,
323  // for which we have not yet distributed
324  // valid indices. depending on which finte
325  // element is dominating the other on this
326  // interface, we either have to constrain
327  // the valid to the invalid indices, or vice
328  // versa.
329  //
330  // we only store an identity if we are about
331  // to overwrite a valid DoF. we will skip
332  // constraining invalid DoFs for now, and
333  // consider them later in Phase 5.
334  if (slave_dof_index != numbers::invalid_dof_index)
335  {
336  // if the DoF indices of both elements
337  // are already distributed, i.e., both
338  // of these 'fe_indices' are associated
339  // with a locally owned cell, then we
340  // should either not have a dof_identity
341  // yet, or it must come out here to be
342  // exactly as we had computed before
343  if (master_dof_index !=
345  Assert((dof_identities.find(
346  master_dof_index) ==
347  dof_identities.end()) ||
348  (dof_identities[slave_dof_index] ==
349  master_dof_index),
350  ExcInternalError());
351 
352  dof_identities[slave_dof_index] =
353  master_dof_index;
354  }
355  }
356  }
357  }
358  }
359 
360  return dof_identities;
361  }
362 
363 
368  template <int spacedim>
369  static std::map<types::global_dof_index, types::global_dof_index>
371  {
372  return std::map<types::global_dof_index, types::global_dof_index>();
373  }
374 
375 
376  template <int dim, int spacedim>
377  static std::map<types::global_dof_index, types::global_dof_index>
379  const hp::DoFHandler<dim, spacedim> &dof_handler)
380  {
381  std::map<types::global_dof_index, types::global_dof_index>
382  dof_identities;
383 
384  // we will mark lines that we have already treated, so first save and
385  // clear the user flags on lines and later restore them
386  std::vector<bool> user_flags;
387  dof_handler.get_triangulation().save_user_flags_line(user_flags);
388  const_cast<::Triangulation<dim, spacedim> &>(
389  dof_handler.get_triangulation())
390  .clear_user_flags_line();
391 
392  // An implementation of the algorithm described in the hp paper,
393  // including the modification mentioned later in the "complications in
394  // 3-d" subsections
395  //
396  // as explained there, we do something only if there are exactly 2
397  // finite elements associated with an object. if there is only one,
398  // then there is nothing to do anyway, and if there are 3 or more,
399  // then we can get into trouble. note that this only happens for lines
400  // in 3d and higher, and for quads only in 4d and higher, so this
401  // isn't a particularly frequent case
402  //
403  // there is one case, however, that we would like to handle (see, for
404  // example, the hp/crash_15 testcase): if we have
405  // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
406  // then we should be able to handle this because we can simply unify
407  // *all* dofs, not only a some. so what we do is to first treat all
408  // pairs of finite elements that have *identical* dofs, and then only
409  // deal with those that are not identical of which we can handle at
410  // most 2
411  ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
412  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
413 
414  for (const auto &cell : dof_handler.active_cell_iterators())
415  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
416  if (cell->line(l)->user_flag_set() == false)
417  {
419  line = cell->line(l);
420  line->set_user_flag();
421 
422  unsigned int unique_sets_of_dofs =
423  line->n_active_fe_indices();
424 
425  // do a first loop over all sets of dofs and do identity
426  // uniquification
427  const unsigned int n_active_fe_indices =
428  line->n_active_fe_indices();
429  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
430  for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
431  {
432  const unsigned int fe_index_1 =
433  line->nth_active_fe_index(f),
434  fe_index_2 =
435  line->nth_active_fe_index(g);
436 
437  // as described in the hp paper, we only unify on lines
438  // when there are at most two different FE objects
439  // assigned on it.
440  // however, more than two 'active_fe_indices' can be
441  // attached that still fulfill the above criterion,
442  // i.e. when two different FiniteElement objects are
443  // assigned to neighboring cells that map their degrees
444  // of freedom one-to-one.
445  // we cannot verify with certainty if two dofs each of
446  // separate FiniteElement objects actually map
447  // one-to-one. however, checking for the number of
448  // 'dofs_per_line' turned out to be a reasonable
449  // approach, that also works for e.g. two different
450  // FE_Q objects of the same order, from which one is
451  // enhanced by a bubble function that is zero on the
452  // boundary.
453  if ((dof_handler.get_fe(fe_index_1).dofs_per_line ==
454  dof_handler.get_fe(fe_index_2).dofs_per_line) &&
455  (dof_handler.get_fe(fe_index_1).dofs_per_line > 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).dofs_per_line;
460 
461  ensure_existence_of_dof_identities<1>(
462  dof_handler.get_fe(fe_index_1),
463  dof_handler.get_fe(fe_index_2),
464  line_dof_identities[fe_index_1][fe_index_2]);
465  // see if these sets of dofs are identical. the
466  // first condition for this is that indeed there are
467  // n identities
468  if (line_dof_identities[fe_index_1][fe_index_2]
469  ->size() == dofs_per_line)
470  {
471  unsigned int i = 0;
472  for (; i < dofs_per_line; ++i)
473  if (((*(line_dof_identities[fe_index_1]
474  [fe_index_2]))[i]
475  .first != i) &&
476  ((*(line_dof_identities[fe_index_1]
477  [fe_index_2]))[i]
478  .second != i))
479  // not an identity
480  break;
481 
482  if (i == dofs_per_line)
483  {
484  // The line dofs (i.e., the ones interior to
485  // a line) of these two finite elements are
486  // identical. Note that there could be
487  // situations when one element still
488  // dominates another, e.g.: FE_Q(2) x
489  // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
490 
491  --unique_sets_of_dofs;
492 
493  // determine which one of both finite
494  // elements is the dominating one.
495  const std::set<unsigned int> fe_indices{
496  fe_index_1, fe_index_2};
497 
498  unsigned int dominating_fe_index =
499  dof_handler.get_fe_collection()
500  .find_dominating_fe(fe_indices,
501  /*codim=*/dim - 1);
502  unsigned int other_fe_index =
504 
505  if (dominating_fe_index !=
507  other_fe_index =
508  (dominating_fe_index == fe_index_1) ?
509  fe_index_2 :
510  fe_index_1;
511  else
512  {
513  // if we haven't found a dominating
514  // finite element, choose the one with
515  // the lower index to be dominating
516  dominating_fe_index = fe_index_1;
517  other_fe_index = fe_index_2;
518  }
519 
520  for (unsigned int j = 0; j < dofs_per_line;
521  ++j)
522  {
524  master_dof_index = line->dof_index(
525  j, dominating_fe_index);
527  slave_dof_index =
528  line->dof_index(j, other_fe_index);
529 
530  // on subdomain boundaries, we will
531  // encounter invalid DoFs on ghost
532  // cells, for which we have not yet
533  // distributed valid indices. depending
534  // on which finte element is dominating
535  // the other on this interface, we
536  // either have to constrain the valid to
537  // the invalid indices, or vice versa.
538  //
539  // we only store an identity if we are
540  // about to overwrite a valid DoF. we
541  // will skip constraining invalid DoFs
542  // for now, and consider them later in
543  // Phase 5.
544  if (slave_dof_index !=
546  {
547  if (master_dof_index !=
549  {
550  // if master dof was already
551  // constrained, constrain to
552  // that one, otherwise constrain
553  // slave to master
554  if (dof_identities.find(
555  master_dof_index) !=
556  dof_identities.end())
557  {
558  // if the DoF indices of
559  // both elements are already
560  // distributed, i.e., both
561  // of these 'fe_indices' are
562  // associated with a locally
563  // owned cell, then we
564  // should either not have a
565  // dof_identity yet, or it
566  // must come out here to be
567  // exactly as we had
568  // computed before
569  Assert(
570  dof_identities.find(
571  dof_identities
572  [master_dof_index]) ==
573  dof_identities.end(),
574  ExcInternalError());
575 
576  dof_identities
577  [slave_dof_index] =
578  dof_identities
579  [master_dof_index];
580  }
581  else
582  {
583  // see comment above for an
584  // explanation of this
585  // assertion
586  Assert(
587  (dof_identities.find(
588  master_dof_index) ==
589  dof_identities.end()) ||
590  (dof_identities
591  [slave_dof_index] ==
592  master_dof_index),
593  ExcInternalError());
594 
595  dof_identities
596  [slave_dof_index] =
597  master_dof_index;
598  }
599  }
600  else
601  {
602  // set slave_dof to
603  // master_dof_index, which is
604  // invalid
605  dof_identities
606  [slave_dof_index] =
608  }
609  }
610  }
611  }
612  }
613  }
614  }
615 
616  // if at this point, there is only one unique set of dofs
617  // left, then we have taken care of everything above. if there
618  // are two, then we need to deal with them here. if there are
619  // more, then we punt, as described in the paper (and
620  // mentioned above)
621  // TODO: The check for 'dim==2' was inserted by intuition. It
622  // fixes
623  // the previous problems with @ref step_27 "step-27" in 3D. But an
624  // explanation for this is still required, and what we do here
625  // is not what we describe in the paper!.
626  if ((unique_sets_of_dofs == 2) && (dim == 2))
627  {
628  const std::set<unsigned int> fe_indices =
629  line->get_active_fe_indices();
630 
631  // find out which is the most dominating finite element of
632  // the ones that are used on this line
633  const unsigned int most_dominating_fe_index =
634  dof_handler.get_fe_collection().find_dominating_fe(
635  fe_indices,
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 !=
645  {
646  // loop over the indices of all the finite elements
647  // that are not dominating, and identify their dofs to
648  // the most dominating one
649  for (const auto &other_fe_index : fe_indices)
650  if (other_fe_index != most_dominating_fe_index)
651  {
652  ensure_existence_of_dof_identities<1>(
653  dof_handler.get_fe(most_dominating_fe_index),
654  dof_handler.get_fe(other_fe_index),
655  line_dof_identities[most_dominating_fe_index]
656  [other_fe_index]);
657 
658  DoFIdentities &identities =
659  *line_dof_identities[most_dominating_fe_index]
660  [other_fe_index];
661  for (const auto &identity : identities)
662  {
664  master_dof_index = line->dof_index(
665  identity.first,
666  most_dominating_fe_index);
668  slave_dof_index =
669  line->dof_index(identity.second,
670  other_fe_index);
671 
672  // on subdomain boundaries, we will
673  // encounter invalid DoFs on ghost cells,
674  // for which we have not yet distributed
675  // valid indices. depending on which finte
676  // element is dominating the other on this
677  // interface, we either have to constrain
678  // the valid to the invalid indices, or vice
679  // versa.
680  //
681  // we only store an identity if we are about
682  // to overwrite a valid DoF. we will skip
683  // constraining invalid DoFs for now, and
684  // consider them later in Phase 5.
685  if (slave_dof_index !=
687  {
688  // if the DoF indices of both elements
689  // are already distributed, i.e., both
690  // of these 'fe_indices' are associated
691  // with a locally owned cell, then we
692  // should either not have a dof_identity
693  // yet, or it must come out here to be
694  // exactly as we had computed before
695  if (master_dof_index !=
697  Assert((dof_identities.find(
698  master_dof_index) ==
699  dof_identities.end()) ||
700  (dof_identities
701  [slave_dof_index] ==
702  master_dof_index),
703  ExcInternalError());
704 
705  dof_identities[slave_dof_index] =
706  master_dof_index;
707  }
708  }
709  }
710  }
711  }
712  }
713 
714  // finally restore the user flags
715  const_cast<::Triangulation<dim, spacedim> &>(
716  dof_handler.get_triangulation())
717  .load_user_flags_line(user_flags);
718 
719  return dof_identities;
720  }
721 
722 
723 
728  template <int dim, int spacedim>
729  static std::map<types::global_dof_index, types::global_dof_index>
731  {
732  // this function should only be called for dim<3 where there are
733  // no quad dof identies. 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  const hp::DoFHandler<3, spacedim> &dof_handler)
745  {
746  const int dim = 3;
747 
748  std::map<types::global_dof_index, types::global_dof_index>
749  dof_identities;
750 
751 
752  // we will mark quads that we have already treated, so first
753  // save and clear the user flags on quads and later restore
754  // them
755  std::vector<bool> user_flags;
756  dof_handler.get_triangulation().save_user_flags_quad(user_flags);
757  const_cast<::Triangulation<dim, spacedim> &>(
758  dof_handler.get_triangulation())
759  .clear_user_flags_quad();
760 
761  // An implementation of the algorithm described in the hp
762  // paper, including the modification mentioned later in the
763  // "complications in 3-d" subsections
764  //
765  // as explained there, we do something only if there are
766  // exactly 2 finite elements associated with an object. if
767  // there is only one, then there is nothing to do anyway,
768  // and if there are 3 or more, then we can get into
769  // trouble. note that this only happens for lines in 3d and
770  // higher, and for quads only in 4d and higher, so this
771  // isn't a particularly frequent case
772  ::Table<2, std::unique_ptr<DoFIdentities>> quad_dof_identities(
773  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
774 
775  for (const auto &cell : dof_handler.active_cell_iterators())
776  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
777  if ((cell->quad(q)->user_flag_set() == false) &&
778  (cell->quad(q)->n_active_fe_indices() == 2))
779  {
781  quad = cell->quad(q);
782  quad->set_user_flag();
783 
784  const std::set<unsigned int> fe_indices =
785  quad->get_active_fe_indices();
786 
787  // find out which is the most dominating finite
788  // element of the ones that are used on this quad
789  const unsigned int most_dominating_fe_index =
790  dof_handler.get_fe_collection().find_dominating_fe(
791  fe_indices,
792  /*codim=*/dim - 2);
793 
794  // if we found the most dominating element, then use
795  // this to eliminate some of the degrees of freedom
796  // by identification. otherwise, the code that
797  // computes hanging node constraints will have to
798  // deal with it by computing appropriate constraints
799  // along this face/edge
800  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
801  {
802  // loop over the indices of all the finite
803  // elements that are not dominating, and
804  // identify their dofs to the most dominating
805  // one
806  for (const auto &other_fe_index : fe_indices)
807  if (other_fe_index != most_dominating_fe_index)
808  {
809  ensure_existence_of_dof_identities<2>(
810  dof_handler.get_fe(most_dominating_fe_index),
811  dof_handler.get_fe(other_fe_index),
812  quad_dof_identities[most_dominating_fe_index]
813  [other_fe_index]);
814 
815  DoFIdentities &identities =
816  *quad_dof_identities[most_dominating_fe_index]
817  [other_fe_index];
818  for (const auto &identity : identities)
819  {
820  const types::global_dof_index master_dof_index =
821  quad->dof_index(identity.first,
822  most_dominating_fe_index);
823  const types::global_dof_index slave_dof_index =
824  quad->dof_index(identity.second,
825  other_fe_index);
826 
827  // we only store an identity if we are about to
828  // overwrite a valid degree of freedom. we will
829  // skip invalid degrees of freedom (that are
830  // associated with ghost cells) for now, and
831  // consider them later in phase 5.
832  if (slave_dof_index !=
834  {
835  // if the DoF indices of both elements are
836  // already distributed, i.e., both of these
837  // 'fe_indices' are associated with a
838  // locally owned cell, then we should either
839  // not have a dof_identity yet, or it must
840  // come out here to be exactly as we had
841  // computed before
842  if (master_dof_index !=
844  Assert(
845  (dof_identities.find(
846  master_dof_index) ==
847  dof_identities.end()) ||
848  (dof_identities[slave_dof_index] ==
849  master_dof_index),
850  ExcInternalError());
851 
852  dof_identities[slave_dof_index] =
853  master_dof_index;
854  }
855  }
856  }
857  }
858  }
859 
860  // finally restore the user flags
861  const_cast<::Triangulation<dim, spacedim> &>(
862  dof_handler.get_triangulation())
863  .load_user_flags_quad(user_flags);
864 
865  return dof_identities;
866  }
867 
868 
869 
874  template <int dim, int spacedim>
875  static void
877  const std::vector<
878  std::map<types::global_dof_index, types::global_dof_index>> &,
880  {}
881 
882 
883  template <int dim, int spacedim>
884  static void
887  &all_constrained_indices,
888  const hp::DoFHandler<dim, spacedim> &dof_handler)
889  {
890  Assert(all_constrained_indices.size() == dim, ExcInternalError());
891 
892  Threads::TaskGroup<> tasks;
893 
894  unsigned int i = 0;
895  tasks += Threads::new_task([&, i]() {
896  all_constrained_indices[i] =
897  compute_vertex_dof_identities(dof_handler);
898  });
899 
900  if (dim > 1)
901  {
902  ++i;
903  tasks += Threads::new_task([&, i]() {
904  all_constrained_indices[i] =
905  compute_line_dof_identities(dof_handler);
906  });
907  }
908 
909  if (dim > 2)
910  {
911  ++i;
912  tasks += Threads::new_task([&, i]() {
913  all_constrained_indices[i] =
914  compute_quad_dof_identities(dof_handler);
915  });
916  }
917 
918  tasks.join_all();
919  }
920 
921 
922 
942  template <class DoFHandlerType>
945  std::vector<types::global_dof_index> &new_dof_indices,
946  const std::vector<
947  std::map<types::global_dof_index, types::global_dof_index>>
948  &all_constrained_indices,
949  const DoFHandlerType &)
950  {
951  Assert(all_constrained_indices.size() == DoFHandlerType::dimension,
952  ExcInternalError());
953 
954  // first preset the new DoF indices that are identities
955  for (const auto &constrained_dof_indices : all_constrained_indices)
956  for (const auto &p : constrained_dof_indices)
957  if (new_dof_indices[p.first] != numbers::invalid_dof_index)
958  {
959  Assert(new_dof_indices[p.first] == enumeration_dof_index,
960  ExcInternalError());
961 
962  new_dof_indices[p.first] = p.second;
963  }
964 
965  // then enumerate the rest
966  types::global_dof_index next_free_dof = 0;
967  for (auto &new_dof_index : new_dof_indices)
968  if (new_dof_index == enumeration_dof_index)
969  new_dof_index = next_free_dof++;
970 
971  // then loop over all those that are constrained and record the
972  // new dof number for those
973  for (const auto &constrained_dof_indices : all_constrained_indices)
974  for (const auto &p : constrained_dof_indices)
975  if (new_dof_indices[p.first] != numbers::invalid_dof_index)
976  {
977  Assert(new_dof_indices[p.first] != enumeration_dof_index,
978  ExcInternalError());
979 
980  if (p.second != numbers::invalid_dof_index)
981  new_dof_indices[p.first] = new_dof_indices[p.second];
982  }
983 
984  for (const types::global_dof_index new_dof_index : new_dof_indices)
985  {
986  (void)new_dof_index;
987  Assert(new_dof_index != enumeration_dof_index,
988  ExcInternalError());
989  Assert(new_dof_index < next_free_dof ||
990  new_dof_index == numbers::invalid_dof_index,
991  ExcInternalError());
992  }
993 
994  return next_free_dof;
995  }
996 
997 
998 
1008  template <int dim, int spacedim>
1011  const unsigned int n_dofs_before_identification,
1012  const bool)
1013  {
1014  return n_dofs_before_identification;
1015  }
1016 
1017 
1018  template <int dim, int spacedim>
1021  const unsigned int n_dofs_before_identification,
1022  const bool check_validity)
1023  {
1024  std::vector<
1025  std::map<types::global_dof_index, types::global_dof_index>>
1026  all_constrained_indices(dim);
1027  compute_dof_identities(all_constrained_indices, dof_handler);
1028 
1029  std::vector<::types::global_dof_index> renumbering(
1030  n_dofs_before_identification, enumeration_dof_index);
1031  const types::global_dof_index n_dofs =
1032  enumerate_dof_indices_for_renumbering(renumbering,
1033  all_constrained_indices,
1034  dof_handler);
1035 
1036  renumber_dofs(renumbering, IndexSet(0), dof_handler, check_validity);
1037 
1038  update_all_active_cell_dof_indices_caches(dof_handler);
1039 
1040  return n_dofs;
1041  }
1042 
1043 
1044 
1049  template <int dim, int spacedim>
1050  static void
1052  hp::DoFHandler<dim, spacedim> &dof_handler)
1053  {
1054  // Note: we may wish to have something here similar to what
1055  // we do for lines and quads, namely that we only identify
1056  // dofs for any fe towards the most dominating one. however,
1057  // it is not clear whether this is actually necessary for
1058  // vertices at all, I can't think of a finite element that
1059  // would make that necessary...
1061  vertex_dof_identities(dof_handler.get_fe_collection().size(),
1062  dof_handler.get_fe_collection().size());
1063 
1064  // mark all vertices on ghost cells
1065  std::vector<bool> include_vertex(
1066  dof_handler.get_triangulation().n_vertices(), false);
1067  if (dynamic_cast<const ::parallel::
1068  DistributedTriangulationBase<dim, spacedim> *>(
1069  &dof_handler.get_triangulation()) != nullptr)
1070  for (const auto &cell : dof_handler.active_cell_iterators())
1071  if (cell->is_ghost())
1072  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
1073  include_vertex[cell->vertex_index(v)] = true;
1074 
1075  // loop over all vertices and see which one we need to work on
1076  for (unsigned int vertex_index = 0;
1077  vertex_index < dof_handler.get_triangulation().n_vertices();
1078  ++vertex_index)
1079  if ((dof_handler.get_triangulation()
1080  .get_used_vertices()[vertex_index] == true) &&
1081  (include_vertex[vertex_index] == true))
1082  {
1083  const unsigned int n_active_fe_indices =
1084  ::internal::DoFAccessorImplementation::Implementation::
1085  n_active_vertex_fe_indices(dof_handler, vertex_index);
1086 
1087  if (n_active_fe_indices > 1)
1088  {
1089  const std::set<unsigned int> fe_indices =
1090  ::internal::DoFAccessorImplementation::
1091  Implementation::get_active_vertex_fe_indices(
1092  dof_handler, vertex_index);
1093 
1094  // find out which is the most dominating finite
1095  // element of the ones that are used on this vertex
1096  const unsigned int most_dominating_fe_index =
1097  dof_handler.get_fe_collection().find_dominating_fe(
1098  fe_indices,
1099  /*codim=*/dim);
1100 
1101  // if we found the most dominating element, then use
1102  // this to eliminate some of the degrees of freedom
1103  // by identification. otherwise, the code that
1104  // computes hanging node constraints will have to
1105  // deal with it by computing appropriate constraints
1106  // along this face/edge
1107  if (most_dominating_fe_index !=
1109  {
1110  // loop over the indices of all the finite
1111  // elements that are not dominating, and
1112  // identify their dofs to the most dominating
1113  // one
1114  for (const auto &other_fe_index : fe_indices)
1115  if (other_fe_index != most_dominating_fe_index)
1116  {
1117  // make sure the entry in the equivalence
1118  // table exists
1119  ensure_existence_of_dof_identities<0>(
1120  dof_handler.get_fe(most_dominating_fe_index),
1121  dof_handler.get_fe(other_fe_index),
1122  vertex_dof_identities[most_dominating_fe_index]
1123  [other_fe_index]);
1124 
1125  // then loop through the identities we
1126  // have. first get the global numbers of the
1127  // dofs we want to identify and make sure they
1128  // are not yet constrained to anything else,
1129  // except for to each other. use the rule that
1130  // we will always constrain the dof with the
1131  // higher fe index to the one with the lower,
1132  // to avoid circular reasoning.
1133  DoFIdentities &identities =
1134  *vertex_dof_identities[most_dominating_fe_index]
1135  [other_fe_index];
1136  for (const auto &identity : identities)
1137  {
1139  master_dof_index = ::internal::
1140  DoFAccessorImplementation::
1141  Implementation::get_vertex_dof_index(
1142  dof_handler,
1143  vertex_index,
1144  most_dominating_fe_index,
1145  identity.first);
1147  slave_dof_index = ::internal::
1148  DoFAccessorImplementation::
1149  Implementation::get_vertex_dof_index(
1150  dof_handler,
1151  vertex_index,
1152  other_fe_index,
1153  identity.second);
1154 
1155  // check if we are on an interface between
1156  // a locally owned and a ghost cell on which
1157  // we need to work on.
1158  //
1159  // all degrees of freedom belonging to
1160  // dominating fe indices or to a processor
1161  // with a higher rank have been set at this
1162  // point (either in Phase 2, or after the
1163  // first ghost exchange in Phase 5). thus,
1164  // we only have to set the indices of
1165  // degrees of freedom that have been
1166  // previously flagged invalid.
1167  if ((slave_dof_index ==
1169  (master_dof_index !=
1171  ::internal::
1172  DoFAccessorImplementation::
1173  Implementation::set_vertex_dof_index(
1174  dof_handler,
1175  vertex_index,
1176  other_fe_index,
1177  identity.second,
1178  master_dof_index);
1179  }
1180  }
1181  }
1182  }
1183  }
1184  }
1185 
1186 
1187 
1192  template <int spacedim>
1195  {}
1196 
1197 
1198  template <int dim, int spacedim>
1199  static void
1201  hp::DoFHandler<dim, spacedim> &dof_handler)
1202  {
1203  // we will mark lines that we have already treated, so first save and
1204  // clear the user flags on lines and later restore them
1205  std::vector<bool> user_flags;
1206  dof_handler.get_triangulation().save_user_flags_line(user_flags);
1207  const_cast<::Triangulation<dim, spacedim> &>(
1208  dof_handler.get_triangulation())
1209  .clear_user_flags_line();
1210 
1211  // mark all lines on ghost cells
1212  for (const auto &cell : dof_handler.active_cell_iterators())
1213  if (cell->is_ghost())
1214  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
1215  ++l)
1216  cell->line(l)->set_user_flag();
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 (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1242  if ((cell->is_locally_owned()) &&
1243  (cell->line(l)->user_flag_set() == true))
1244  {
1246  line = cell->line(l);
1247  line->clear_user_flag();
1248 
1249  unsigned int unique_sets_of_dofs =
1250  line->n_active_fe_indices();
1251 
1252  // do a first loop over all sets of dofs and do identity
1253  // uniquification
1254  const unsigned int n_active_fe_indices =
1255  line->n_active_fe_indices();
1256  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1257  for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
1258  {
1259  const unsigned int fe_index_1 =
1260  line->nth_active_fe_index(f),
1261  fe_index_2 =
1262  line->nth_active_fe_index(g);
1263 
1264  if ((dof_handler.get_fe(fe_index_1).dofs_per_line ==
1265  dof_handler.get_fe(fe_index_2).dofs_per_line) &&
1266  (dof_handler.get_fe(fe_index_1).dofs_per_line > 0))
1267  {
1268  // the number of dofs per line is identical
1269  const unsigned int dofs_per_line =
1270  dof_handler.get_fe(fe_index_1).dofs_per_line;
1271 
1272  ensure_existence_of_dof_identities<1>(
1273  dof_handler.get_fe(fe_index_1),
1274  dof_handler.get_fe(fe_index_2),
1275  line_dof_identities[fe_index_1][fe_index_2]);
1276  // see if these sets of dofs are identical. the
1277  // first condition for this is that indeed there are
1278  // n identities
1279  if (line_dof_identities[fe_index_1][fe_index_2]
1280  ->size() == dofs_per_line)
1281  {
1282  unsigned int i = 0;
1283  for (; i < dofs_per_line; ++i)
1284  if (((*(line_dof_identities[fe_index_1]
1285  [fe_index_2]))[i]
1286  .first != i) &&
1287  ((*(line_dof_identities[fe_index_1]
1288  [fe_index_2]))[i]
1289  .second != i))
1290  // not an identity
1291  break;
1292 
1293  if (i == dofs_per_line)
1294  {
1295  // The line dofs (i.e., the ones interior to
1296  // a line) of these two finite elements are
1297  // identical. Note that there could be
1298  // situations when one element still
1299  // dominates another, e.g.: FE_Q(2) x
1300  // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
1301 
1302  --unique_sets_of_dofs;
1303 
1304  // determine which one of both finite
1305  // elements is the dominating one.
1306  const std::set<unsigned int> fe_indices{
1307  fe_index_1, fe_index_2};
1308 
1309  unsigned int dominating_fe_index =
1310  dof_handler.get_fe_collection()
1311  .find_dominating_fe(fe_indices,
1312  /*codim*/ dim - 1);
1313  unsigned int 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  master_dof_index = line->dof_index(
1336  j, dominating_fe_index);
1338  slave_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 ((slave_dof_index ==
1356  (master_dof_index !=
1358  line->set_dof_index(j,
1359  master_dof_index,
1360  fe_index_2);
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<unsigned int> 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  const unsigned int most_dominating_fe_index =
1385  dof_handler.get_fe_collection().find_dominating_fe(
1386  fe_indices,
1387  /*codim=*/dim - 1);
1388 
1389  // if we found the most dominating element, then use this
1390  // to eliminate some of the degrees of freedom by
1391  // identification. otherwise, the code that computes
1392  // hanging node constraints will have to deal with it by
1393  // computing appropriate constraints along this face/edge
1394  if (most_dominating_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  ensure_existence_of_dof_identities<1>(
1404  dof_handler.get_fe(most_dominating_fe_index),
1405  dof_handler.get_fe(other_fe_index),
1406  line_dof_identities[most_dominating_fe_index]
1407  [other_fe_index]);
1408 
1409  DoFIdentities &identities =
1410  *line_dof_identities[most_dominating_fe_index]
1411  [other_fe_index];
1412  for (const auto &identity : identities)
1413  {
1415  master_dof_index = line->dof_index(
1416  identity.first,
1417  most_dominating_fe_index);
1419  slave_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 ((slave_dof_index ==
1437  (master_dof_index !=
1439  line->set_dof_index(identity.second,
1440  master_dof_index,
1441  other_fe_index);
1442  }
1443  }
1444  }
1445  }
1446  }
1447 
1448  // finally restore the user flags
1449  const_cast<::Triangulation<dim, spacedim> &>(
1450  dof_handler.get_triangulation())
1451  .load_user_flags_line(user_flags);
1452  }
1453 
1454 
1455 
1460  template <int dim, int spacedim>
1461  static void
1464  {
1465  // this function should only be called for dim<3 where there are
1466  // no quad dof identies. for dim>=3, the specialization below should
1467  // take care of it
1468  Assert(dim < 3, ExcInternalError());
1469  }
1470 
1471 
1472  template <int spacedim>
1474  hp::DoFHandler<3, spacedim> &dof_handler)
1475  {
1476  const int dim = 3;
1477 
1478  // we will mark quads that we have already treated, so first
1479  // save and clear the user flags on quads and later restore
1480  // them
1481  std::vector<bool> user_flags;
1482  dof_handler.get_triangulation().save_user_flags_quad(user_flags);
1483  const_cast<::Triangulation<dim, spacedim> &>(
1484  dof_handler.get_triangulation())
1485  .clear_user_flags_quad();
1486 
1487  // mark all quads on ghost cells
1488  for (const auto &cell : dof_handler.active_cell_iterators())
1489  if (cell->is_ghost())
1490  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell;
1491  ++q)
1492  cell->quad(q)->set_user_flag();
1493 
1494  // An implementation of the algorithm described in the hp
1495  // paper, including the modification mentioned later in the
1496  // "complications in 3-d" subsections
1497  //
1498  // as explained there, we do something only if there are
1499  // exactly 2 finite elements associated with an object. if
1500  // there is only one, then there is nothing to do anyway,
1501  // and if there are 3 or more, then we can get into
1502  // trouble. note that this only happens for lines in 3d and
1503  // higher, and for quads only in 4d and higher, so this
1504  // isn't a particularly frequent case
1505  ::Table<2, std::unique_ptr<DoFIdentities>> quad_dof_identities(
1506  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
1507 
1508  for (const auto &cell : dof_handler.active_cell_iterators())
1509  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
1510  if ((cell->is_locally_owned()) &&
1511  (cell->quad(q)->user_flag_set() == true) &&
1512  (cell->quad(q)->n_active_fe_indices() == 2))
1513  {
1515  quad = cell->quad(q);
1516  quad->clear_user_flag();
1517 
1518  const std::set<unsigned int> fe_indices =
1519  quad->get_active_fe_indices();
1520 
1521  // find out which is the most dominating finite
1522  // element of the ones that are used on this quad
1523  const unsigned int most_dominating_fe_index =
1524  dof_handler.get_fe_collection().find_dominating_fe(
1525  fe_indices,
1526  /*codim=*/dim - 2);
1527 
1528  // if we found the most dominating element, then use
1529  // this to eliminate some of the degrees of freedom
1530  // by identification. otherwise, the code that
1531  // computes hanging node constraints will have to
1532  // deal with it by computing appropriate constraints
1533  // along this face/edge
1534  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
1535  {
1536  // loop over the indices of all the finite
1537  // elements that are not dominating, and
1538  // identify their dofs to the most dominating
1539  // one
1540  for (const auto &other_fe_index : fe_indices)
1541  if (other_fe_index != most_dominating_fe_index)
1542  {
1543  ensure_existence_of_dof_identities<2>(
1544  dof_handler.get_fe(most_dominating_fe_index),
1545  dof_handler.get_fe(other_fe_index),
1546  quad_dof_identities[most_dominating_fe_index]
1547  [other_fe_index]);
1548 
1549  DoFIdentities &identities =
1550  *quad_dof_identities[most_dominating_fe_index]
1551  [other_fe_index];
1552  for (const auto &identity : identities)
1553  {
1554  const types::global_dof_index master_dof_index =
1555  quad->dof_index(identity.first,
1556  most_dominating_fe_index);
1557  const types::global_dof_index slave_dof_index =
1558  quad->dof_index(identity.second,
1559  other_fe_index);
1560 
1561  // check if we are on an interface between
1562  // a locally owned and a ghost cell on which
1563  // we need to work on.
1564  //
1565  // all degrees of freedom belonging to
1566  // dominating fe indices or to a processor with
1567  // a higher rank have been set at this point
1568  // (either in Phase 2, or after the first ghost
1569  // exchange in Phase 5). thus, we only have to
1570  // set the indices of degrees of freedom that
1571  // have been previously flagged invalid.
1572  if ((slave_dof_index ==
1574  (master_dof_index !=
1576  quad->set_dof_index(identity.second,
1577  master_dof_index,
1578  other_fe_index);
1579  }
1580  }
1581  }
1582  }
1583 
1584  // finally restore the user flags
1585  const_cast<::Triangulation<dim, spacedim> &>(
1586  dof_handler.get_triangulation())
1587  .load_user_flags_quad(user_flags);
1588  }
1589 
1590 
1591 
1604  template <int dim, int spacedim>
1605  static void
1607  const DoFHandler<dim, spacedim> &)
1608  {}
1609 
1610 
1611  template <int dim, int spacedim>
1612  static void
1614  hp::DoFHandler<dim, spacedim> &dof_handler)
1615  {
1616  {
1617  Threads::TaskGroup<> tasks;
1618 
1619  tasks += Threads::new_task([&]() {
1620  merge_invalid_vertex_dofs_on_ghost_interfaces(dof_handler);
1621  });
1622 
1623  if (dim > 1)
1624  {
1625  tasks += Threads::new_task([&]() {
1626  merge_invalid_line_dofs_on_ghost_interfaces(dof_handler);
1627  });
1628  }
1629 
1630  if (dim > 2)
1631  {
1632  tasks += Threads::new_task([&]() {
1633  merge_invalid_quad_dofs_on_ghost_interfaces(dof_handler);
1634  });
1635  }
1636 
1637  tasks.join_all();
1638  }
1639 
1640  update_all_active_cell_dof_indices_caches(dof_handler);
1641  }
1642 
1643 
1644 
1651  template <class DoFHandlerType>
1654  DoFHandlerType & dof_handler)
1655  {
1656  Assert(dof_handler.get_triangulation().n_levels() > 0,
1657  ExcMessage("Empty triangulation"));
1658 
1659  // Step 1: distribute dofs on all cells, but definitely
1660  // exclude artificial cells
1661  types::global_dof_index next_free_dof = 0;
1662 
1663  std::vector<types::global_dof_index> dof_indices;
1664 
1665  for (auto cell : dof_handler.active_cell_iterators())
1666  if (!cell->is_artificial())
1667  if ((subdomain_id == numbers::invalid_subdomain_id) ||
1668  (cell->subdomain_id() == subdomain_id))
1669  {
1670  dof_indices.resize(cell->get_fe().dofs_per_cell);
1671 
1672  // circumvent cache
1673  internal::DoFAccessorImplementation::Implementation::
1674  get_dof_indices(*cell,
1675  dof_indices,
1676  cell->active_fe_index());
1677 
1678  for (auto &dof_index : dof_indices)
1679  if (dof_index == numbers::invalid_dof_index)
1680  dof_index = next_free_dof++;
1681 
1682  cell->set_dof_indices(dof_indices);
1683  }
1684 
1685  update_all_active_cell_dof_indices_caches(dof_handler);
1686 
1687  return next_free_dof;
1688  }
1689 
1690 
1691 
1705  template <class DoFHandlerType>
1706  static void
1708  std::vector<types::global_dof_index> &renumbering,
1710  const DoFHandlerType & dof_handler)
1711  {
1712  std::vector<types::global_dof_index> local_dof_indices;
1713 
1714  for (const auto &cell : dof_handler.active_cell_iterators())
1715  if (cell->is_ghost() && (cell->subdomain_id() < subdomain_id))
1716  {
1717  // we found a neighboring ghost cell whose subdomain
1718  // is "stronger" than our own subdomain
1719 
1720  // delete all dofs that live there and that we have
1721  // previously assigned a number to (i.e. the ones on
1722  // the interface)
1723  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
1724  cell->get_dof_indices(local_dof_indices);
1725  for (const auto &local_dof_index : local_dof_indices)
1726  if (local_dof_index != numbers::invalid_dof_index)
1727  renumbering[local_dof_index] = numbers::invalid_dof_index;
1728  }
1729  }
1730 
1731 
1732 
1733  /* -------------- distribute_mg_dofs functionality ------------- */
1734 
1735 
1736 
1737  template <class DoFHandlerType>
1740  DoFHandlerType & dof_handler,
1741  const unsigned int level)
1742  {
1743  Assert(DoFHandlerType::is_hp_dof_handler == false,
1744  ExcInternalError());
1745 
1746  const unsigned int dim = DoFHandlerType::dimension;
1747  const unsigned int spacedim = DoFHandlerType::space_dimension;
1748 
1749  const ::Triangulation<dim, spacedim> &tria =
1750  dof_handler.get_triangulation();
1751  Assert(tria.n_levels() > 0, ExcMessage("Empty triangulation"));
1752  if (level >= tria.n_levels())
1753  return 0; // this is allowed for multigrid
1754 
1755  types::global_dof_index next_free_dof = 0;
1756 
1757  std::vector<types::global_dof_index> dof_indices;
1758 
1759  for (auto cell : dof_handler.cell_iterators_on_level(level))
1760  if ((level_subdomain_id == numbers::invalid_subdomain_id) ||
1761  (cell->level_subdomain_id() == level_subdomain_id))
1762  {
1763  dof_indices.resize(cell->get_fe().dofs_per_cell);
1764 
1765  cell->get_mg_dof_indices(dof_indices);
1766 
1767  for (auto &dof_index : dof_indices)
1768  if (dof_index == numbers::invalid_dof_index)
1769  dof_index = next_free_dof++;
1770 
1771  cell->set_mg_dof_indices(dof_indices);
1772  }
1773 
1774  return next_free_dof;
1775  }
1776 
1777 
1778 
1779  /* --------------------- renumber_dofs functionality ---------------- */
1780 
1781 
1789  template <int dim, int spacedim>
1790  static void
1792  const std::vector<types::global_dof_index> &new_numbers,
1793  const IndexSet & indices_we_care_about,
1794  DoFHandler<dim, spacedim> & dof_handler,
1795  const bool check_validity)
1796  {
1797  // we can not use cell iterators in this function since then
1798  // we would renumber the dofs on the interface of two cells
1799  // more than once. Anyway, this way it's not only more
1800  // correct but also faster; note, however, that dof numbers
1801  // may be invalid_dof_index, namely when the appropriate
1802  // vertex/line/etc is unused
1803  for (std::vector<types::global_dof_index>::iterator i =
1804  dof_handler.vertex_dofs.begin();
1805  i != dof_handler.vertex_dofs.end();
1806  ++i)
1807  if (*i != numbers::invalid_dof_index)
1808  *i = (indices_we_care_about.size() == 0) ?
1809  (new_numbers[*i]) :
1810  (new_numbers[indices_we_care_about.index_within_set(*i)]);
1811  else if (check_validity)
1812  // if index is invalid_dof_index: check if this one
1813  // really is unused
1814  Assert(dof_handler.get_triangulation().vertex_used(
1815  (i - dof_handler.vertex_dofs.begin()) /
1816  dof_handler.get_fe().dofs_per_vertex) == false,
1817  ExcInternalError());
1818  }
1819 
1820 
1821 
1829  template <int dim, int spacedim>
1830  static void
1832  const std::vector<types::global_dof_index> &new_numbers,
1833  const IndexSet & indices_we_care_about,
1834  DoFHandler<dim, spacedim> & dof_handler)
1835  {
1836  for (unsigned int level = 0; level < dof_handler.levels.size();
1837  ++level)
1838  for (std::vector<types::global_dof_index>::iterator i =
1839  dof_handler.levels[level]->dof_object.dofs.begin();
1840  i != dof_handler.levels[level]->dof_object.dofs.end();
1841  ++i)
1842  if (*i != numbers::invalid_dof_index)
1843  *i =
1844  ((indices_we_care_about.size() == 0) ?
1845  new_numbers[*i] :
1846  new_numbers[indices_we_care_about.index_within_set(*i)]);
1847  }
1848 
1849 
1850 
1858  template <int spacedim>
1859  static void
1861  const std::vector<types::global_dof_index> & /*new_numbers*/,
1862  const IndexSet & /*indices_we_care_about*/,
1863  DoFHandler<1, spacedim> & /*dof_handler*/)
1864  {
1865  // nothing to do in 1d since there are no separate faces
1866  }
1867 
1868 
1869 
1870  template <int spacedim>
1871  static void
1873  const std::vector<types::global_dof_index> &new_numbers,
1874  const IndexSet & indices_we_care_about,
1875  DoFHandler<2, spacedim> & dof_handler)
1876  {
1877  // treat dofs on lines
1878  for (std::vector<types::global_dof_index>::iterator i =
1879  dof_handler.faces->lines.dofs.begin();
1880  i != dof_handler.faces->lines.dofs.end();
1881  ++i)
1882  if (*i != numbers::invalid_dof_index)
1883  *i = ((indices_we_care_about.size() == 0) ?
1884  new_numbers[*i] :
1885  new_numbers[indices_we_care_about.index_within_set(*i)]);
1886  }
1887 
1888 
1889 
1890  template <int spacedim>
1891  static void
1893  const std::vector<types::global_dof_index> &new_numbers,
1894  const IndexSet & indices_we_care_about,
1895  DoFHandler<3, spacedim> & dof_handler)
1896  {
1897  // treat dofs on lines
1898  for (std::vector<types::global_dof_index>::iterator i =
1899  dof_handler.faces->lines.dofs.begin();
1900  i != dof_handler.faces->lines.dofs.end();
1901  ++i)
1902  if (*i != numbers::invalid_dof_index)
1903  *i = ((indices_we_care_about.size() == 0) ?
1904  new_numbers[*i] :
1905  new_numbers[indices_we_care_about.index_within_set(*i)]);
1906 
1907  // treat dofs on quads
1908  for (std::vector<types::global_dof_index>::iterator i =
1909  dof_handler.faces->quads.dofs.begin();
1910  i != dof_handler.faces->quads.dofs.end();
1911  ++i)
1912  if (*i != numbers::invalid_dof_index)
1913  *i = ((indices_we_care_about.size() == 0) ?
1914  new_numbers[*i] :
1915  new_numbers[indices_we_care_about.index_within_set(*i)]);
1916  }
1917 
1918 
1919 
1920  template <int dim, int spacedim>
1921  static void
1923  const std::vector<types::global_dof_index> &new_numbers,
1924  const IndexSet & indices_we_care_about,
1925  hp::DoFHandler<dim, spacedim> & dof_handler,
1926  const bool check_validity)
1927  {
1928  for (unsigned int vertex_index = 0;
1929  vertex_index < dof_handler.get_triangulation().n_vertices();
1930  ++vertex_index)
1931  {
1932  const unsigned int n_active_fe_indices =
1933  ::internal::DoFAccessorImplementation::Implementation::
1934  n_active_vertex_fe_indices(dof_handler, vertex_index);
1935 
1936  // if this vertex is unused, then we really ought not to have
1937  // allocated any space for it, i.e., n_active_fe_indices should be
1938  // zero, and there is no space to actually store dof indices for
1939  // this vertex
1940  if (dof_handler.get_triangulation().vertex_used(vertex_index) ==
1941  false)
1942  Assert(n_active_fe_indices == 0, ExcInternalError());
1943 
1944  // otherwise the vertex is used; it may still not hold any dof
1945  // indices if it is located on an artificial cell and not adjacent
1946  // to a ghost cell, but in that case there is simply nothing for
1947  // us to do
1948  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1949  {
1950  const unsigned int fe_index =
1951  ::internal::DoFAccessorImplementation::
1952  Implementation::nth_active_vertex_fe_index(dof_handler,
1953  vertex_index,
1954  f);
1955 
1956  for (unsigned int d = 0;
1957  d < dof_handler.get_fe(fe_index).dofs_per_vertex;
1958  ++d)
1959  {
1960  const types::global_dof_index old_dof_index =
1961  ::internal::DoFAccessorImplementation::
1962  Implementation::get_vertex_dof_index(dof_handler,
1963  vertex_index,
1964  fe_index,
1965  d);
1966 
1967  // if check_validity was set, then we are to verify that
1968  // the previous indices were all valid. this really should
1969  // be the case: we allocated space for these vertex dofs,
1970  // i.e., at least one adjacent cell has a valid
1971  // active_fe_index, so there are DoFs that really live
1972  // on this vertex. if check_validity is set, then we
1973  // must make sure that they have been set to something
1974  // useful
1975  if (check_validity)
1976  Assert(old_dof_index != numbers::invalid_dof_index,
1977  ExcInternalError());
1978 
1979  if (old_dof_index != numbers::invalid_dof_index)
1980  {
1981  // In the following blocks, we first check whether
1982  // we were given an IndexSet of DoFs to touch. If not
1983  // (the first 'if' case here), then we are in the
1984  // sequential case and are allowed to touch all DoFs.
1985  //
1986  // If yes (the 'else' case), then we need to
1987  // distinguish whether the DoF whose number we want to
1988  // touch is in fact locally owned (i.e., is in the
1989  // index set) and then we can actually assign it a new
1990  // number; otherwise, we have encountered a
1991  // non-locally owned DoF for which we don't know the
1992  // new number yet and so set it to an invalid index.
1993  // This will later be fixed up after the first ghost
1994  // exchange phase when we unify hp DoFs on neighboring
1995  // cells.
1996  if (indices_we_care_about.size() == 0)
1998  Implementation::set_vertex_dof_index(
1999  dof_handler,
2000  vertex_index,
2001  fe_index,
2002  d,
2003  new_numbers[old_dof_index]);
2004  else
2005  {
2006  if (indices_we_care_about.is_element(
2007  old_dof_index))
2009  Implementation::set_vertex_dof_index(
2010  dof_handler,
2011  vertex_index,
2012  fe_index,
2013  d,
2014  new_numbers[indices_we_care_about
2015  .index_within_set(
2016  old_dof_index)]);
2017  else
2018  ::internal::DoFAccessorImplementation::
2019  Implementation::set_vertex_dof_index(
2020  dof_handler,
2021  vertex_index,
2022  fe_index,
2023  d,
2025  }
2026  }
2027  }
2028  }
2029  }
2030  }
2031 
2032 
2033 
2034  template <int dim, int spacedim>
2035  static void
2037  const std::vector<types::global_dof_index> &new_numbers,
2038  const IndexSet & indices_we_care_about,
2039  hp::DoFHandler<dim, spacedim> & dof_handler)
2040  {
2041  for (const auto &cell : dof_handler.active_cell_iterators())
2042  if (!cell->is_artificial())
2043  {
2044  const unsigned int fe_index = cell->active_fe_index();
2045 
2046  for (unsigned int d = 0;
2047  d < dof_handler.get_fe(fe_index)
2048  .template n_dofs_per_object<dim>();
2049  ++d)
2050  {
2051  const types::global_dof_index old_dof_index =
2052  cell->dof_index(d, fe_index);
2053  if (old_dof_index != numbers::invalid_dof_index)
2054  {
2055  // In the following blocks, we first check whether
2056  // we were given an IndexSet of DoFs to touch. If not
2057  // (the first 'if' case here), then we are in the
2058  // sequential case and are allowed to touch all DoFs.
2059  //
2060  // If yes (the 'else' case), then we need to distinguish
2061  // whether the DoF whose number we want to touch is in
2062  // fact locally owned (i.e., is in the index set) and
2063  // then we can actually assign it a new number;
2064  // otherwise, we have encountered a non-locally owned
2065  // DoF for which we don't know the new number yet and so
2066  // set it to an invalid index. This will later be fixed
2067  // up after the first ghost exchange phase when we unify
2068  // hp DoFs on neighboring cells.
2069  if (indices_we_care_about.size() == 0)
2070  cell->set_dof_index(d,
2071  new_numbers[old_dof_index],
2072  fe_index);
2073  else
2074  {
2075  if (indices_we_care_about.is_element(old_dof_index))
2076  cell->set_dof_index(
2077  d,
2078  new_numbers[indices_we_care_about
2079  .index_within_set(old_dof_index)],
2080  fe_index);
2081  else
2082  cell->set_dof_index(d,
2084  fe_index);
2085  }
2086  }
2087  }
2088  }
2089  }
2090 
2091 
2092 
2093  template <int spacedim>
2094  static void
2096  const std::vector<types::global_dof_index> & /*new_numbers*/,
2097  const IndexSet & /*indices_we_care_about*/,
2098  hp::DoFHandler<1, spacedim> & /*dof_handler*/)
2099  {
2100  // nothing to do in 1d since there are no separate faces -- we've
2101  // already taken care of this when dealing with the vertices
2102  }
2103 
2104 
2105 
2106  template <int spacedim>
2107  static void
2109  const std::vector<types::global_dof_index> &new_numbers,
2110  const IndexSet & indices_we_care_about,
2111  hp::DoFHandler<2, spacedim> & dof_handler)
2112  {
2113  const unsigned int dim = 2;
2114 
2115  // deal with DoFs on lines
2116  {
2117  // save user flags on lines so we can use them to mark lines
2118  // we've already treated
2119  std::vector<bool> saved_line_user_flags;
2120  const_cast<::Triangulation<dim, spacedim> &>(
2121  dof_handler.get_triangulation())
2122  .save_user_flags_line(saved_line_user_flags);
2123  const_cast<::Triangulation<dim, spacedim> &>(
2124  dof_handler.get_triangulation())
2125  .clear_user_flags_line();
2126 
2127  for (const auto &cell : dof_handler.active_cell_iterators())
2128  if (!cell->is_artificial())
2129  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
2130  ++l)
2131  if (cell->line(l)->user_flag_set() == false)
2132  {
2133  const typename hp::DoFHandler<dim,
2134  spacedim>::line_iterator
2135  line = cell->line(l);
2136  line->set_user_flag();
2137 
2138  const unsigned int n_active_fe_indices =
2139  line->n_active_fe_indices();
2140 
2141  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2142  {
2143  const unsigned int fe_index =
2144  line->nth_active_fe_index(f);
2145 
2146  for (unsigned int d = 0;
2147  d < dof_handler.get_fe(fe_index).dofs_per_line;
2148  ++d)
2149  {
2150  const types::global_dof_index old_dof_index =
2151  line->dof_index(d, fe_index);
2152  if (old_dof_index != numbers::invalid_dof_index)
2153  {
2154  // In the following blocks, we first check
2155  // whether we were given an IndexSet of DoFs
2156  // to touch. If not (the first 'if' case
2157  // here), then we are in the sequential case
2158  // and are allowed to touch all DoFs.
2159  //
2160  // If yes (the 'else' case), then we need to
2161  // distinguish whether the DoF whose number we
2162  // want to touch is in fact locally owned
2163  // (i.e., is in the index set) and then we can
2164  // actually assign it a new number; otherwise,
2165  // we have encountered a non-locally owned DoF
2166  // for which we don't know the new number yet
2167  // and so set it to an invalid index. This
2168  // will later be fixed up after the first
2169  // ghost exchange phase when we unify hp DoFs
2170  // on neighboring cells.
2171  if (indices_we_care_about.size() == 0)
2172  line->set_dof_index(
2173  d, new_numbers[old_dof_index], fe_index);
2174  else
2175  {
2176  if (indices_we_care_about.is_element(
2177  old_dof_index))
2178  line->set_dof_index(
2179  d,
2180  new_numbers[indices_we_care_about
2182  old_dof_index)],
2183  fe_index);
2184  else
2185  line->set_dof_index(
2186  d,
2188  fe_index);
2189  }
2190  }
2191  }
2192  }
2193  }
2194 
2195  // at the end, restore the user
2196  // flags for the lines
2197  const_cast<::Triangulation<dim, spacedim> &>(
2198  dof_handler.get_triangulation())
2199  .load_user_flags_line(saved_line_user_flags);
2200  }
2201  }
2202 
2203 
2204 
2205  template <int spacedim>
2206  static void
2208  const std::vector<types::global_dof_index> &new_numbers,
2209  const IndexSet & indices_we_care_about,
2210  hp::DoFHandler<3, spacedim> & dof_handler)
2211  {
2212  const unsigned int dim = 3;
2213 
2214  // deal with DoFs on lines
2215  {
2216  // save user flags on lines so we can use them to mark lines
2217  // we've already treated
2218  std::vector<bool> saved_line_user_flags;
2219  const_cast<::Triangulation<dim, spacedim> &>(
2220  dof_handler.get_triangulation())
2221  .save_user_flags_line(saved_line_user_flags);
2222  const_cast<::Triangulation<dim, spacedim> &>(
2223  dof_handler.get_triangulation())
2224  .clear_user_flags_line();
2225 
2226  for (const auto &cell : dof_handler.active_cell_iterators())
2227  if (!cell->is_artificial())
2228  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
2229  ++l)
2230  if (cell->line(l)->user_flag_set() == false)
2231  {
2232  const typename hp::DoFHandler<dim,
2233  spacedim>::line_iterator
2234  line = cell->line(l);
2235  line->set_user_flag();
2236 
2237  const unsigned int n_active_fe_indices =
2238  line->n_active_fe_indices();
2239 
2240  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2241  {
2242  const unsigned int fe_index =
2243  line->nth_active_fe_index(f);
2244 
2245  for (unsigned int d = 0;
2246  d < dof_handler.get_fe(fe_index).dofs_per_line;
2247  ++d)
2248  {
2249  const types::global_dof_index old_dof_index =
2250  line->dof_index(d, fe_index);
2251  if (old_dof_index != numbers::invalid_dof_index)
2252  {
2253  // In the following blocks, we first check
2254  // whether we were given an IndexSet of DoFs
2255  // to touch. If not (the first 'if' case
2256  // here), then we are in the sequential case
2257  // and are allowed to touch all DoFs.
2258  //
2259  // If yes (the 'else' case), then we need to
2260  // distinguish whether the DoF whose number we
2261  // want to touch is in fact locally owned
2262  // (i.e., is in the index set) and then we can
2263  // actually assign it a new number; otherwise,
2264  // we have encountered a non-locally owned DoF
2265  // for which we don't know the new number yet
2266  // and so set it to an invalid index. This
2267  // will later be fixed up after the first
2268  // ghost exchange phase when we unify hp DoFs
2269  // on neighboring cells.
2270  if (indices_we_care_about.size() == 0)
2271  line->set_dof_index(
2272  d, new_numbers[old_dof_index], fe_index);
2273  else if (indices_we_care_about.is_element(
2274  old_dof_index))
2275  line->set_dof_index(
2276  d,
2277  new_numbers[indices_we_care_about
2279  old_dof_index)],
2280  fe_index);
2281  else
2282  line->set_dof_index(
2283  d, numbers::invalid_dof_index, fe_index);
2284  }
2285  }
2286  }
2287  }
2288 
2289  // at the end, restore the user
2290  // flags for the lines
2291  const_cast<::Triangulation<dim, spacedim> &>(
2292  dof_handler.get_triangulation())
2293  .load_user_flags_line(saved_line_user_flags);
2294  }
2295 
2296  // then deal with dofs on quads
2297  {
2298  std::vector<bool> saved_quad_user_flags;
2299  const_cast<::Triangulation<dim, spacedim> &>(
2300  dof_handler.get_triangulation())
2301  .save_user_flags_quad(saved_quad_user_flags);
2302  const_cast<::Triangulation<dim, spacedim> &>(
2303  dof_handler.get_triangulation())
2304  .clear_user_flags_quad();
2305 
2306  for (const auto &cell : dof_handler.active_cell_iterators())
2307  if (!cell->is_artificial())
2308  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell;
2309  ++q)
2310  if (cell->quad(q)->user_flag_set() == false)
2311  {
2312  const typename hp::DoFHandler<dim,
2313  spacedim>::quad_iterator
2314  quad = cell->quad(q);
2315  quad->set_user_flag();
2316 
2317  const unsigned int n_active_fe_indices =
2318  quad->n_active_fe_indices();
2319 
2320  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2321  {
2322  const unsigned int fe_index =
2323  quad->nth_active_fe_index(f);
2324 
2325  for (unsigned int d = 0;
2326  d < dof_handler.get_fe(fe_index).dofs_per_quad;
2327  ++d)
2328  {
2329  const types::global_dof_index old_dof_index =
2330  quad->dof_index(d, fe_index);
2331  if (old_dof_index != numbers::invalid_dof_index)
2332  {
2333  // In the following blocks, we first check
2334  // whether we were given an IndexSet of DoFs
2335  // to touch. If not (the first 'if' case
2336  // here), then we are in the sequential case
2337  // and are allowed to touch all DoFs.
2338  //
2339  // If yes (the 'else' case), then we need to
2340  // distinguish whether the DoF whose number we
2341  // want to touch is in fact locally owned
2342  // (i.e., is in the index set) and then we can
2343  // actually assign it a new number; otherwise,
2344  // we have encountered a non-locally owned DoF
2345  // for which we don't know the new number yet
2346  // and so set it to an invalid index. This
2347  // will later be fixed up after the first
2348  // ghost exchange phase when we unify hp DoFs
2349  // on neighboring cells.
2350  if (indices_we_care_about.size() == 0)
2351  quad->set_dof_index(
2352  d, new_numbers[old_dof_index], fe_index);
2353  else
2354  {
2355  if (indices_we_care_about.is_element(
2356  old_dof_index))
2357  quad->set_dof_index(
2358  d,
2359  new_numbers[indices_we_care_about
2361  old_dof_index)],
2362  fe_index);
2363  else
2364  quad->set_dof_index(
2365  d,
2367  fe_index);
2368  }
2369  }
2370  }
2371  }
2372  }
2373 
2374  // at the end, restore the user flags for the quads
2375  const_cast<::Triangulation<dim, spacedim> &>(
2376  dof_handler.get_triangulation())
2377  .load_user_flags_quad(saved_quad_user_flags);
2378  }
2379  }
2380 
2381 
2382 
2394  template <class DoFHandlerType>
2395  static void
2396  renumber_dofs(const std::vector<types::global_dof_index> &new_numbers,
2397  const IndexSet &indices_we_care_about,
2398  DoFHandlerType &dof_handler,
2399  const bool check_validity)
2400  {
2401  if (DoFHandlerType::dimension == 1)
2402  Assert(indices_we_care_about == IndexSet(0), ExcNotImplemented());
2403 
2404  // renumber DoF indices on vertices, cells, and faces. this
2405  // can be done in parallel because the respective functions
2406  // work on separate data structures
2407  Threads::TaskGroup<> tasks;
2408  tasks += Threads::new_task([&]() {
2409  renumber_vertex_dofs(new_numbers,
2410  indices_we_care_about,
2411  dof_handler,
2412  check_validity);
2413  });
2414  tasks += Threads::new_task([&]() {
2415  renumber_face_dofs(new_numbers, indices_we_care_about, dof_handler);
2416  });
2417  tasks += Threads::new_task([&]() {
2418  renumber_cell_dofs(new_numbers, indices_we_care_about, dof_handler);
2419  });
2420  tasks.join_all();
2421 
2422  // update the cache used for cell dof indices
2423  update_all_active_cell_dof_indices_caches(dof_handler);
2424  }
2425 
2426 
2427 
2428  /* --------------------- renumber_mg_dofs functionality ----------------
2429  */
2430 
2438  template <int dim, int spacedim>
2439  static void
2441  const std::vector<::types::global_dof_index> &new_numbers,
2442  const IndexSet & indices_we_care_about,
2443  DoFHandler<dim, spacedim> &dof_handler,
2444  const unsigned int level,
2445  const bool check_validity)
2446  {
2447  (void)check_validity;
2448  Assert(level < dof_handler.get_triangulation().n_levels(),
2449  ExcInternalError());
2450 
2451  for (typename std::vector<
2452  typename DoFHandler<dim, spacedim>::MGVertexDoFs>::iterator i =
2453  dof_handler.mg_vertex_dofs.begin();
2454  i != dof_handler.mg_vertex_dofs.end();
2455  ++i)
2456  // if the present vertex lives on the current level
2457  if ((i->get_coarsest_level() <= level) &&
2458  (i->get_finest_level() >= level))
2459  for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_vertex;
2460  ++d)
2461  {
2463  i->get_index(level,
2464  d,
2465  dof_handler.get_fe().dofs_per_vertex);
2466 
2467  if (idx != numbers::invalid_dof_index)
2468  {
2469  Assert(check_validity == false ||
2470  (indices_we_care_about.size() > 0 ?
2471  indices_we_care_about.is_element(idx) :
2472  (idx < new_numbers.size())),
2473  ExcInternalError());
2474  i->set_index(level,
2475  d,
2476  dof_handler.get_fe().dofs_per_vertex,
2477  (indices_we_care_about.size() == 0) ?
2478  (new_numbers[idx]) :
2479  (new_numbers[indices_we_care_about
2480  .index_within_set(idx)]));
2481  }
2482  }
2483  }
2484 
2485 
2486 
2494  template <int dim, int spacedim>
2495  static void
2497  const std::vector<::types::global_dof_index> &new_numbers,
2498  const IndexSet & indices_we_care_about,
2499  DoFHandler<dim, spacedim> &dof_handler,
2500  const unsigned int level)
2501  {
2502  for (std::vector<types::global_dof_index>::iterator i =
2503  dof_handler.mg_levels[level]->dof_object.dofs.begin();
2504  i != dof_handler.mg_levels[level]->dof_object.dofs.end();
2505  ++i)
2506  {
2507  if (*i != numbers::invalid_dof_index)
2508  {
2509  Assert((indices_we_care_about.size() > 0 ?
2510  indices_we_care_about.is_element(*i) :
2511  (*i < new_numbers.size())),
2512  ExcInternalError());
2513  *i =
2514  (indices_we_care_about.size() == 0) ?
2515  (new_numbers[*i]) :
2516  (new_numbers[indices_we_care_about.index_within_set(*i)]);
2517  }
2518  }
2519  }
2520 
2521 
2522 
2530  template <int spacedim>
2531  static void
2533  const std::vector<types::global_dof_index> & /*new_numbers*/,
2534  const IndexSet & /*indices_we_care_about*/,
2535  DoFHandler<1, spacedim> & /*dof_handler*/,
2536  const unsigned int /*level*/,
2537  const bool /*check_validity*/)
2538  {
2539  // nothing to do in 1d because there are no separate faces
2540  }
2541 
2542 
2543 
2544  template <int spacedim>
2545  static void
2547  const std::vector<::types::global_dof_index> &new_numbers,
2548  const IndexSet & indices_we_care_about,
2549  DoFHandler<2, spacedim> &dof_handler,
2550  const unsigned int level,
2551  const bool check_validity)
2552  {
2553  if (dof_handler.get_fe().dofs_per_line > 0)
2554  {
2555  // save user flags as they will be modified
2556  std::vector<bool> user_flags;
2557  dof_handler.get_triangulation().save_user_flags(user_flags);
2558  const_cast<::Triangulation<2, spacedim> &>(
2559  dof_handler.get_triangulation())
2560  .clear_user_flags();
2561 
2562  // flag all lines adjacent to cells of the current
2563  // level, as those lines logically belong to the same
2564  // level as the cell, at least for for isotropic
2565  // refinement
2567  endc = dof_handler.end(level);
2568  for (cell = dof_handler.begin(level); cell != endc; ++cell)
2569  if (cell->level_subdomain_id() !=
2571  for (const unsigned int line :
2573  cell->face(line)->set_user_flag();
2574 
2575  for (typename DoFHandler<2, spacedim>::cell_iterator cell =
2576  dof_handler.begin();
2577  cell != dof_handler.end();
2578  ++cell)
2579  for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell;
2580  ++l)
2581  if (cell->line(l)->user_flag_set())
2582  {
2583  for (unsigned int d = 0;
2584  d < dof_handler.get_fe().dofs_per_line;
2585  ++d)
2586  {
2588  cell->line(l)->mg_dof_index(level, d);
2589  if (check_validity)
2591  ExcInternalError());
2592 
2593  if (idx != numbers::invalid_dof_index)
2594  cell->line(l)->set_mg_dof_index(
2595  level,
2596  d,
2597  ((indices_we_care_about.size() == 0) ?
2598  new_numbers[idx] :
2599  new_numbers[indices_we_care_about
2600  .index_within_set(idx)]));
2601  }
2602  cell->line(l)->clear_user_flag();
2603  }
2604  // finally, restore user flags
2605  const_cast<::Triangulation<2, spacedim> &>(
2606  dof_handler.get_triangulation())
2607  .load_user_flags(user_flags);
2608  }
2609  }
2610 
2611 
2612 
2613  template <int spacedim>
2614  static void
2616  const std::vector<::types::global_dof_index> &new_numbers,
2617  const IndexSet & indices_we_care_about,
2618  DoFHandler<3, spacedim> &dof_handler,
2619  const unsigned int level,
2620  const bool check_validity)
2621  {
2622  if (dof_handler.get_fe().dofs_per_line > 0 ||
2623  dof_handler.get_fe().dofs_per_quad > 0)
2624  {
2625  // save user flags as they will be modified
2626  std::vector<bool> user_flags;
2627  dof_handler.get_triangulation().save_user_flags(user_flags);
2628  const_cast<::Triangulation<3, spacedim> &>(
2629  dof_handler.get_triangulation())
2630  .clear_user_flags();
2631 
2632  // flag all lines adjacent to cells of the current
2633  // level, as those lines logically belong to the same
2634  // level as the cell, at least for isotropic refinement
2636  endc = dof_handler.end(level);
2637  for (cell = dof_handler.begin(level); cell != endc; ++cell)
2638  if (cell->level_subdomain_id() !=
2640  for (unsigned int line = 0;
2641  line < GeometryInfo<3>::lines_per_cell;
2642  ++line)
2643  cell->line(line)->set_user_flag();
2644 
2645  for (typename DoFHandler<3, spacedim>::cell_iterator cell =
2646  dof_handler.begin();
2647  cell != dof_handler.end();
2648  ++cell)
2649  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell;
2650  ++l)
2651  if (cell->line(l)->user_flag_set())
2652  {
2653  for (unsigned int d = 0;
2654  d < dof_handler.get_fe().dofs_per_line;
2655  ++d)
2656  {
2658  cell->line(l)->mg_dof_index(level, d);
2659  if (check_validity)
2661  ExcInternalError());
2662 
2663  if (idx != numbers::invalid_dof_index)
2664  cell->line(l)->set_mg_dof_index(
2665  level,
2666  d,
2667  ((indices_we_care_about.size() == 0) ?
2668  new_numbers[idx] :
2669  new_numbers[indices_we_care_about
2670  .index_within_set(idx)]));
2671  }
2672  cell->line(l)->clear_user_flag();
2673  }
2674 
2675  // flag all quads adjacent to cells of the current level, as
2676  // those quads logically belong to the same level as the cell,
2677  // at least for isotropic refinement
2678  for (cell = dof_handler.begin(level); cell != endc; ++cell)
2679  if (cell->level_subdomain_id() !=
2681  for (unsigned int quad = 0;
2682  quad < GeometryInfo<3>::quads_per_cell;
2683  ++quad)
2684  cell->quad(quad)->set_user_flag();
2685 
2686  for (typename DoFHandler<3, spacedim>::cell_iterator cell =
2687  dof_handler.begin();
2688  cell != dof_handler.end();
2689  ++cell)
2690  for (unsigned int l = 0; l < GeometryInfo<3>::quads_per_cell;
2691  ++l)
2692  if (cell->quad(l)->user_flag_set())
2693  {
2694  for (unsigned int d = 0;
2695  d < dof_handler.get_fe().dofs_per_quad;
2696  ++d)
2697  {
2699  cell->quad(l)->mg_dof_index(level, d);
2700  if (check_validity)
2702  ExcInternalError());
2703 
2704  if (idx != numbers::invalid_dof_index)
2705  cell->quad(l)->set_mg_dof_index(
2706  level,
2707  d,
2708  ((indices_we_care_about.size() == 0) ?
2709  new_numbers[idx] :
2710  new_numbers[indices_we_care_about
2711  .index_within_set(idx)]));
2712  }
2713  cell->quad(l)->clear_user_flag();
2714  }
2715 
2716  // finally, restore user flags
2717  const_cast<::Triangulation<3, spacedim> &>(
2718  dof_handler.get_triangulation())
2719  .load_user_flags(user_flags);
2720  }
2721  }
2722 
2723 
2724 
2725  template <int dim, int spacedim>
2726  static void
2728  const std::vector<::types::global_dof_index> &new_numbers,
2729  const IndexSet & indices_we_care_about,
2730  DoFHandler<dim, spacedim> &dof_handler,
2731  const unsigned int level,
2732  const bool check_validity)
2733  {
2734  Assert(level < dof_handler.get_triangulation().n_global_levels(),
2735  ExcInternalError());
2736 
2737  // renumber DoF indices on vertices, cells, and faces. this
2738  // can be done in parallel because the respective functions
2739  // work on separate data structures
2740  Threads::TaskGroup<> tasks;
2741  tasks += Threads::new_task([&]() {
2742  renumber_vertex_mg_dofs(new_numbers,
2743  indices_we_care_about,
2744  dof_handler,
2745  level,
2746  check_validity);
2747  });
2748  tasks += Threads::new_task([&]() {
2749  renumber_face_mg_dofs(new_numbers,
2750  indices_we_care_about,
2751  dof_handler,
2752  level,
2753  check_validity);
2754  });
2755  tasks += Threads::new_task([&]() {
2756  renumber_cell_mg_dofs(new_numbers,
2757  indices_we_care_about,
2758  dof_handler,
2759  level);
2760  });
2761  tasks.join_all();
2762  }
2763 
2764 
2765 
2766  template <int dim, int spacedim>
2767  static void
2769  const std::vector<::types::global_dof_index> & /*new_numbers*/,
2770  const IndexSet & /*indices_we_care_about*/,
2771  hp::DoFHandler<dim, spacedim> & /*dof_handler*/,
2772  const unsigned int /*level*/,
2773  const bool /*check_validity*/)
2774  {
2775  Assert(false, ExcNotImplemented());
2776  }
2777  };
2778 
2779 
2780 
2781  /* --------------------- class Sequential ---------------- */
2782 
2783 
2784 
2785  template <class DoFHandlerType>
2786  Sequential<DoFHandlerType>::Sequential(DoFHandlerType &dof_handler)
2787  : dof_handler(&dof_handler)
2788  {}
2789 
2790 
2791 
2792  template <class DoFHandlerType>
2793  NumberCache
2795  {
2796  const types::global_dof_index n_initial_dofs =
2798  *dof_handler);
2799 
2800  const types::global_dof_index n_dofs =
2802  n_initial_dofs,
2803  /*check_validity=*/true);
2804 
2805  // return a sequential, complete index set
2806  return NumberCache(n_dofs);
2807  }
2808 
2809 
2810 
2811  template <class DoFHandlerType>
2812  std::vector<NumberCache>
2814  {
2815  std::vector<bool> user_flags;
2816  dof_handler->get_triangulation().save_user_flags(user_flags);
2817 
2818  const_cast<::Triangulation<DoFHandlerType::dimension,
2819  DoFHandlerType::space_dimension> &>(
2820  dof_handler->get_triangulation())
2821  .clear_user_flags();
2822 
2823  std::vector<NumberCache> number_caches;
2824  number_caches.reserve(dof_handler->get_triangulation().n_levels());
2825  for (unsigned int level = 0;
2826  level < dof_handler->get_triangulation().n_levels();
2827  ++level)
2828  {
2829  // first distribute dofs on this level
2830  const types::global_dof_index n_level_dofs =
2833 
2834  // then add a complete, sequential index set
2835  number_caches.emplace_back(n_level_dofs);
2836  }
2837 
2838  const_cast<::Triangulation<DoFHandlerType::dimension,
2839  DoFHandlerType::space_dimension> &>(
2840  dof_handler->get_triangulation())
2841  .load_user_flags(user_flags);
2842 
2843  return number_caches;
2844  }
2845 
2846 
2847 
2848  template <class DoFHandlerType>
2849  NumberCache
2851  const std::vector<types::global_dof_index> &new_numbers) const
2852  {
2853  Implementation::renumber_dofs(new_numbers,
2854  IndexSet(0),
2855  *dof_handler,
2856  /*check_validity=*/true);
2857 
2858  // return a sequential, complete index set. take into account that the
2859  // number of DoF indices may in fact be smaller than there were before
2860  // if some previously separately numbered dofs have been identified.
2861  // this is, for example, what the hp::DoFHandler does: it first
2862  // enumerates all DoFs on cells independently, and then unifies
2863  // some located at vertices or faces; this leaves us with fewer
2864  // DoFs than there were before, so use the largest index as
2865  // the one to determine the size of the index space
2866  return NumberCache(
2867  *std::max_element(new_numbers.begin(), new_numbers.end()) + 1);
2868  }
2869 
2870 
2871 
2872  template <class DoFHandlerType>
2873  NumberCache
2875  const unsigned int level,
2876  const std::vector<types::global_dof_index> &new_numbers) const
2877  {
2879  new_numbers, IndexSet(0), *dof_handler, level, true);
2880 
2881  // return a sequential, complete index set
2882  return NumberCache(new_numbers.size());
2883  }
2884 
2885 
2886  /* --------------------- class ParallelShared ---------------- */
2887 
2888 
2889  template <class DoFHandlerType>
2891  DoFHandlerType &dof_handler)
2892  : dof_handler(&dof_handler)
2893  {}
2894 
2895 
2896 
2897  namespace
2898  {
2907  template <class DoFHandlerType>
2908  std::vector<types::subdomain_id>
2909  get_dof_subdomain_association(const DoFHandlerType & dof_handler,
2910  const types::global_dof_index n_dofs,
2911  const unsigned int n_procs)
2912  {
2913  (void)n_procs;
2914  std::vector<types::subdomain_id> subdomain_association(
2916  std::vector<types::global_dof_index> local_dof_indices;
2917  local_dof_indices.reserve(DoFTools::max_dofs_per_cell(dof_handler));
2918 
2919  // loop over all cells and record which subdomain a DoF belongs to.
2920  // give to the smaller subdomain_id in case it is on an interface
2921  typename DoFHandlerType::active_cell_iterator
2922  cell = dof_handler.begin_active(),
2923  endc = dof_handler.end();
2924  for (; cell != endc; ++cell)
2925  {
2926  // get the owner of the cell; note that we have made sure above
2927  // that all cells are either locally owned or ghosts (not
2928  // artificial), so this call will always yield the true owner
2929  const types::subdomain_id subdomain_id = cell->subdomain_id();
2930  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
2931  local_dof_indices.resize(dofs_per_cell);
2932  cell->get_dof_indices(local_dof_indices);
2933 
2934  // set subdomain ids. if dofs already have their values set then
2935  // they must be on partition interfaces. In that case assign them
2936  // to the processor with the smaller subdomain id.
2937  for (unsigned int i = 0; i < dofs_per_cell; ++i)
2938  if (subdomain_association[local_dof_indices[i]] ==
2940  subdomain_association[local_dof_indices[i]] = subdomain_id;
2941  else if (subdomain_association[local_dof_indices[i]] >
2942  subdomain_id)
2943  {
2944  subdomain_association[local_dof_indices[i]] = subdomain_id;
2945  }
2946  }
2947 
2948  Assert(std::find(subdomain_association.begin(),
2949  subdomain_association.end(),
2951  subdomain_association.end(),
2952  ExcInternalError());
2953 
2954  Assert(*std::max_element(subdomain_association.begin(),
2955  subdomain_association.end()) < n_procs,
2956  ExcInternalError());
2957 
2958  return subdomain_association;
2959  }
2960 
2961 
2968  template <class DoFHandlerType>
2969  std::vector<types::subdomain_id>
2970  get_dof_level_subdomain_association(
2971  const DoFHandlerType & dof_handler,
2972  const types::global_dof_index n_dofs_on_level,
2973  const unsigned int n_procs,
2974  const unsigned int level)
2975  {
2976  (void)n_procs;
2977  std::vector<types::subdomain_id> level_subdomain_association(
2978  n_dofs_on_level, numbers::invalid_subdomain_id);
2979  std::vector<types::global_dof_index> local_dof_indices;
2980  local_dof_indices.reserve(DoFTools::max_dofs_per_cell(dof_handler));
2981 
2982  // loop over all cells and record which subdomain a DoF belongs to.
2983  // interface goes to proccessor with smaller subdomain id
2984  typename DoFHandlerType::cell_iterator cell =
2985  dof_handler.begin(level),
2986  endc = dof_handler.end(level);
2987  for (; cell != endc; ++cell)
2988  {
2989  // get the owner of the cell; note that we have made sure above
2990  // that all cells are either locally owned or ghosts (not
2991  // artificial), so this call will always yield the true owner
2992  const types::subdomain_id level_subdomain_id =
2993  cell->level_subdomain_id();
2994  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
2995  local_dof_indices.resize(dofs_per_cell);
2996  cell->get_mg_dof_indices(local_dof_indices);
2997 
2998  // set level subdomain ids. if dofs already have their values set
2999  // then they must be on partition interfaces. In that case assign
3000  // them to the processor with the smaller subdomain id.
3001  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3002  if (level_subdomain_association[local_dof_indices[i]] ==
3004  level_subdomain_association[local_dof_indices[i]] =
3005  level_subdomain_id;
3006  else if (level_subdomain_association[local_dof_indices[i]] >
3007  level_subdomain_id)
3008  {
3009  level_subdomain_association[local_dof_indices[i]] =
3010  level_subdomain_id;
3011  }
3012  }
3013 
3014  Assert(std::find(level_subdomain_association.begin(),
3015  level_subdomain_association.end(),
3017  level_subdomain_association.end(),
3018  ExcInternalError());
3019 
3020  Assert(*std::max_element(level_subdomain_association.begin(),
3021  level_subdomain_association.end()) < n_procs,
3022  ExcInternalError());
3023 
3024  return level_subdomain_association;
3025  }
3026  } // namespace
3027 
3028 
3029 
3030  template <class DoFHandlerType>
3031  NumberCache
3033  {
3034  const unsigned int dim = DoFHandlerType::dimension;
3035  const unsigned int spacedim = DoFHandlerType::space_dimension;
3036 
3037  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3038  (dynamic_cast<
3039  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3040  &this->dof_handler->get_triangulation()));
3041  Assert(tr != nullptr, ExcInternalError());
3042 
3043  const unsigned int n_procs =
3044  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3045 
3046  // If the underlying shared::Tria allows artificial cells,
3047  // then save the current set of subdomain ids, and set
3048  // subdomain ids to the "true" owner of each cell. we later
3049  // restore these flags
3050  std::vector<types::subdomain_id> saved_subdomain_ids;
3051  if (tr->with_artificial_cells())
3052  {
3053  saved_subdomain_ids.resize(tr->n_active_cells());
3054 
3055  const std::vector<types::subdomain_id> &true_subdomain_ids =
3056  tr->get_true_subdomain_ids_of_cells();
3057 
3058  for (const auto &cell : tr->active_cell_iterators())
3059  {
3060  const unsigned int index = cell->active_cell_index();
3061  saved_subdomain_ids[index] = cell->subdomain_id();
3062  cell->set_subdomain_id(true_subdomain_ids[index]);
3063  }
3064  }
3065 
3066  // first let the sequential algorithm do its magic. it is going to
3067  // enumerate DoFs on all cells, regardless of owner
3068  const types::global_dof_index n_initial_dofs =
3070  *this->dof_handler);
3071 
3072  const types::global_dof_index n_dofs =
3074  n_initial_dofs,
3075  /*check_validity=*/true);
3076 
3077  // then re-enumerate them based on their subdomain association.
3078  // for this, we first have to identify for each current DoF
3079  // index which subdomain they belong to. ideally, we would
3080  // like to call DoFRenumbering::subdomain_wise(), but
3081  // because the NumberCache of the current DoFHandler is not
3082  // fully set up yet, we can't quite do that. also, that
3083  // function has to deal with other kinds of triangulations as
3084  // well, whereas we here know what kind of triangulation
3085  // we have and can simplify the code accordingly
3086  std::vector<types::global_dof_index> new_dof_indices(
3087  n_dofs, enumeration_dof_index);
3088  {
3089  // first get the association of each dof with a subdomain and
3090  // determine the total number of subdomain ids used
3091  const std::vector<types::subdomain_id> subdomain_association =
3092  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3093 
3094  // then renumber the subdomains by first looking at those belonging
3095  // to subdomain 0, then those of subdomain 1, etc. note that the
3096  // algorithm is stable, i.e. if two dofs i,j have i<j and belong to
3097  // the same subdomain, then they will be in this order also after
3098  // reordering
3099  types::global_dof_index next_free_index = 0;
3100  for (types::subdomain_id subdomain = 0; subdomain < n_procs;
3101  ++subdomain)
3102  for (types::global_dof_index i = 0; i < n_dofs; ++i)
3103  if (subdomain_association[i] == subdomain)
3104  {
3105  Assert(new_dof_indices[i] == enumeration_dof_index,
3106  ExcInternalError());
3107  new_dof_indices[i] = next_free_index;
3108  ++next_free_index;
3109  }
3110 
3111  // we should have numbered all dofs
3112  Assert(next_free_index == n_dofs, ExcInternalError());
3113  Assert(std::find(new_dof_indices.begin(),
3114  new_dof_indices.end(),
3115  enumeration_dof_index) == new_dof_indices.end(),
3116  ExcInternalError());
3117  }
3118  // finally do the renumbering. we can use the sequential
3119  // version of the function because we do things on all
3120  // cells and all cells have their subdomain ids and DoFs
3121  // correctly set
3122  Implementation::renumber_dofs(new_dof_indices,
3123  IndexSet(0),
3124  *this->dof_handler,
3125  /*check_validity=*/true);
3126 
3127  // update the number cache. for this, we first have to find the
3128  // subdomain association for each DoF again following renumbering, from
3129  // which we can then compute the IndexSets of locally owned DoFs for all
3130  // processors. all other fields then follow from this
3131  //
3132  // given the way we enumerate degrees of freedom, the locally owned
3133  // ranges must all be contiguous and consecutive. this makes filling
3134  // the IndexSets cheap. an assertion at the top verifies that this
3135  // assumption is true
3136  const std::vector<types::subdomain_id> subdomain_association =
3137  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3138 
3139  for (unsigned int i = 1; i < n_dofs; ++i)
3140  Assert(subdomain_association[i] >= subdomain_association[i - 1],
3141  ExcInternalError());
3142 
3143  std::vector<IndexSet> locally_owned_dofs_per_processor(
3144  n_procs, IndexSet(n_dofs));
3145  {
3146  // we know that the set of subdomain indices is contiguous from
3147  // the assertion above; find the start and end index for each
3148  // processor, taking into account that sometimes a processor
3149  // may not in fact have any DoFs at all. we do the latter
3150  // by just identifying contiguous ranges of subdomain_ids
3151  // and filling IndexSets for those subdomains; subdomains
3152  // that don't appear will lead to IndexSets that are simply
3153  // never touched and remain empty as initialized above.
3154  unsigned int start_index = 0;
3155  unsigned int end_index = 0;
3156  while (start_index < n_dofs)
3157  {
3158  while ((end_index) < n_dofs &&
3159  (subdomain_association[end_index] ==
3160  subdomain_association[start_index]))
3161  ++end_index;
3162 
3163  // we've now identified a range of same indices. set that
3164  // range in the corresponding IndexSet
3165  if (end_index > start_index)
3166  {
3167  const unsigned int subdomain_owner =
3168  subdomain_association[start_index];
3169  locally_owned_dofs_per_processor[subdomain_owner].add_range(
3170  start_index, end_index);
3171  }
3172 
3173  // then move on to thinking about the next range
3174  start_index = end_index;
3175  }
3176  }
3177 
3178  // finally, restore current subdomain ids
3179  if (tr->with_artificial_cells())
3180  for (const auto &cell : tr->active_cell_iterators())
3181  cell->set_subdomain_id(
3182  saved_subdomain_ids[cell->active_cell_index()]);
3183 
3184  // return a NumberCache object made up from the sets of locally
3185  // owned DoFs
3186  return NumberCache(
3187  locally_owned_dofs_per_processor,
3188  this->dof_handler->get_triangulation().locally_owned_subdomain());
3189  }
3190 
3191 
3192 
3193  template <class DoFHandlerType>
3194  std::vector<NumberCache>
3196  {
3197  const unsigned int dim = DoFHandlerType::dimension;
3198  const unsigned int spacedim = DoFHandlerType::space_dimension;
3199 
3200  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3201  (dynamic_cast<
3202  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3203  &this->dof_handler->get_triangulation()));
3204  Assert(tr != nullptr, ExcInternalError());
3205 
3206  const unsigned int n_procs =
3207  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3208  const unsigned int n_levels = tr->n_global_levels();
3209 
3210  std::vector<NumberCache> number_caches;
3211  number_caches.reserve(n_levels);
3212 
3213  // We create an index set for each level
3214  for (unsigned int lvl = 0; lvl < n_levels; ++lvl)
3215  {
3216  // If the underlying shared::Tria allows artificial cells,
3217  // then save the current set of level subdomain ids, and set
3218  // subdomain ids to the "true" owner of each cell. we later
3219  // restore these flags
3220  // Note: "allows_artificial_cells" is currently enforced for
3221  // MG computations.
3222  std::vector<types::subdomain_id> saved_level_subdomain_ids;
3223  saved_level_subdomain_ids.resize(tr->n_cells(lvl));
3224  {
3225  typename ::parallel::shared::Triangulation<dim, spacedim>::
3226  cell_iterator cell =
3227  this->dof_handler->get_triangulation().begin(
3228  lvl),
3229  endc =
3230  this->dof_handler->get_triangulation().end(lvl);
3231 
3232  const std::vector<types::subdomain_id> &true_level_subdomain_ids =
3233  tr->get_true_level_subdomain_ids_of_cells(lvl);
3234 
3235  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3236  {
3237  saved_level_subdomain_ids[index] = cell->level_subdomain_id();
3238  cell->set_level_subdomain_id(true_level_subdomain_ids[index]);
3239  }
3240  }
3241 
3242  // Next let the sequential algorithm do its magic. it is going to
3243  // enumerate DoFs on all cells on the given level, regardless of
3244  // owner
3245  const types::global_dof_index n_dofs_on_level =
3248 
3249  // then re-enumerate them based on their level subdomain
3250  // association. for this, we first have to identify for each current
3251  // DoF index which subdomain they belong to. ideally, we would like
3252  // to call DoFRenumbering::subdomain_wise(), but because the
3253  // NumberCache of the current DoFHandler is not fully set up yet, we
3254  // can't quite do that. also, that function has to deal with other
3255  // kinds of triangulations as well, whereas we here know what kind
3256  // of triangulation we have and can simplify the code accordingly
3257  std::vector<types::global_dof_index> new_dof_indices(
3258  n_dofs_on_level, numbers::invalid_dof_index);
3259  {
3260  // first get the association of each dof with a subdomain and
3261  // determine the total number of subdomain ids used
3262  const std::vector<types::subdomain_id>
3263  level_subdomain_association =
3264  get_dof_level_subdomain_association(*this->dof_handler,
3265  n_dofs_on_level,
3266  n_procs,
3267  lvl);
3268 
3269  // then renumber the subdomains by first looking at those
3270  // belonging to subdomain 0, then those of subdomain 1, etc. note
3271  // that the algorithm is stable, i.e. if two dofs i,j have i<j and
3272  // belong to the same subdomain, then they will be in this order
3273  // also after reordering
3274  types::global_dof_index next_free_index = 0;
3275  for (types::subdomain_id level_subdomain = 0;
3276  level_subdomain < n_procs;
3277  ++level_subdomain)
3278  for (types::global_dof_index i = 0; i < n_dofs_on_level; ++i)
3279  if (level_subdomain_association[i] == level_subdomain)
3280  {
3281  Assert(new_dof_indices[i] == numbers::invalid_dof_index,
3282  ExcInternalError());
3283  new_dof_indices[i] = next_free_index;
3284  ++next_free_index;
3285  }
3286 
3287  // we should have numbered all dofs
3288  Assert(next_free_index == n_dofs_on_level, ExcInternalError());
3289  Assert(std::find(new_dof_indices.begin(),
3290  new_dof_indices.end(),
3292  new_dof_indices.end(),
3293  ExcInternalError());
3294  }
3295 
3296  // finally do the renumbering. we can use the sequential
3297  // version of the function because we do things on all
3298  // cells and all cells have their subdomain ids and DoFs
3299  // correctly set
3301  new_dof_indices, IndexSet(0), *this->dof_handler, lvl, true);
3302 
3303  // update the number cache. for this, we first have to find the
3304  // level subdomain association for each DoF again following
3305  // renumbering, from which we can then compute the IndexSets of
3306  // locally owned DoFs for all processors. all other fields then
3307  // follow from this
3308  //
3309  // given the way we enumerate degrees of freedom, the locally owned
3310  // ranges must all be contiguous and consecutive. this makes filling
3311  // the IndexSets cheap. an assertion at the top verifies that this
3312  // assumption is true
3313  const std::vector<types::subdomain_id> level_subdomain_association =
3314  get_dof_level_subdomain_association(*this->dof_handler,
3315  n_dofs_on_level,
3316  n_procs,
3317  lvl);
3318 
3319  for (unsigned int i = 1; i < n_dofs_on_level; ++i)
3320  Assert(level_subdomain_association[i] >=
3321  level_subdomain_association[i - 1],
3322  ExcInternalError());
3323 
3324  std::vector<IndexSet> locally_owned_dofs_per_processor(
3325  n_procs, IndexSet(n_dofs_on_level));
3326  {
3327  // we know that the set of subdomain indices is contiguous from
3328  // the assertion above; find the start and end index for each
3329  // processor, taking into account that sometimes a processor
3330  // may not in fact have any DoFs at all. we do the latter
3331  // by just identifying contiguous ranges of level_subdomain_ids
3332  // and filling IndexSets for those subdomains; subdomains
3333  // that don't appear will lead to IndexSets that are simply
3334  // never touched and remain empty as initialized above.
3335  unsigned int start_index = 0;
3336  unsigned int end_index = 0;
3337  while (start_index < n_dofs_on_level)
3338  {
3339  while ((end_index) < n_dofs_on_level &&
3340  (level_subdomain_association[end_index] ==
3341  level_subdomain_association[start_index]))
3342  ++end_index;
3343 
3344  // we've now identified a range of same indices. set that
3345  // range in the corresponding IndexSet
3346  if (end_index > start_index)
3347  {
3348  const unsigned int level_subdomain_owner =
3349  level_subdomain_association[start_index];
3350  locally_owned_dofs_per_processor[level_subdomain_owner]
3351  .add_range(start_index, end_index);
3352  }
3353 
3354  // then move on to thinking about the next range
3355  start_index = end_index;
3356  }
3357  }
3358 
3359  // finally, restore current level subdomain ids
3360  {
3361  typename ::parallel::shared::Triangulation<dim, spacedim>::
3362  cell_iterator cell =
3363  this->dof_handler->get_triangulation().begin(
3364  lvl),
3365  endc =
3366  this->dof_handler->get_triangulation().end(lvl);
3367 
3368  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3369  cell->set_level_subdomain_id(saved_level_subdomain_ids[index]);
3370 
3371  // add NumberCache for current level
3372  number_caches.emplace_back(
3373  NumberCache(locally_owned_dofs_per_processor,
3374  this->dof_handler->get_triangulation()
3375  .locally_owned_subdomain()));
3376  }
3377  }
3378 
3379  return number_caches;
3380  }
3381 
3382 
3383 
3384  template <class DoFHandlerType>
3385  NumberCache
3387  const std::vector<types::global_dof_index> &new_numbers) const
3388  {
3389 #ifndef DEAL_II_WITH_MPI
3390  (void)new_numbers;
3391  Assert(false, ExcNotImplemented());
3392  return NumberCache();
3393 #else
3394  const unsigned int dim = DoFHandlerType::dimension;
3395  const unsigned int spacedim = DoFHandlerType::space_dimension;
3396 
3397  // Similar to distribute_dofs() we need to have a special treatment in
3398  // case artificial cells are present.
3399  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3400  (dynamic_cast<
3401  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3402  &this->dof_handler->get_triangulation()));
3403  Assert(tr != nullptr, ExcInternalError());
3404 
3405  typename ::parallel::shared::Triangulation<dim, spacedim>::
3406  active_cell_iterator
3407  cell = this->dof_handler->get_triangulation().begin_active(),
3408  endc = this->dof_handler->get_triangulation().end();
3409  std::vector<types::subdomain_id> current_subdomain_ids(
3410  tr->n_active_cells());
3411  const std::vector<types::subdomain_id> &true_subdomain_ids =
3412  tr->get_true_subdomain_ids_of_cells();
3413  if (tr->with_artificial_cells())
3414  for (unsigned int index = 0; cell != endc; cell++, index++)
3415  {
3416  current_subdomain_ids[index] = cell->subdomain_id();
3417  cell->set_subdomain_id(true_subdomain_ids[index]);
3418  }
3419 
3420  std::vector<types::global_dof_index> global_gathered_numbers(
3421  this->dof_handler->n_dofs(), 0);
3422  // as we call DoFRenumbering::subdomain_wise (*dof_handler) from
3423  // distribute_dofs(), we need to support sequential-like input.
3424  // Distributed-like input from, for example, component_wise renumbering
3425  // is also supported.
3426  if (new_numbers.size() == this->dof_handler->n_dofs())
3427  {
3428  global_gathered_numbers = new_numbers;
3429  }
3430  else
3431  {
3432  Assert(new_numbers.size() ==
3433  this->dof_handler->locally_owned_dofs().n_elements(),
3434  ExcInternalError());
3435  const unsigned int n_cpu =
3436  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3437  std::vector<types::global_dof_index> gathered_new_numbers(
3438  this->dof_handler->n_dofs(), 0);
3439  Assert(Utilities::MPI::this_mpi_process(tr->get_communicator()) ==
3440  this->dof_handler->get_triangulation()
3441  .locally_owned_subdomain(),
3442  ExcInternalError())
3443 
3444  // gather new numbers among processors into one vector
3445  {
3446  std::vector<types::global_dof_index> new_numbers_copy(
3447  new_numbers);
3448 
3449  // store the number of elements that are to be received from each
3450  // process
3451  std::vector<int> rcounts(n_cpu);
3452 
3454  // set rcounts based on new_numbers:
3455  int cur_count = new_numbers_copy.size();
3456  int ierr = MPI_Allgather(&cur_count,
3457  1,
3458  MPI_INT,
3459  rcounts.data(),
3460  1,
3461  MPI_INT,
3462  tr->get_communicator());
3463  AssertThrowMPI(ierr);
3464 
3465  // compute the displacements (relative to recvbuf)
3466  // at which to place the incoming data from process i
3467  std::vector<int> displacements(n_cpu);
3468  for (unsigned int i = 0; i < n_cpu; i++)
3469  {
3470  displacements[i] = shift;
3471  shift += rcounts[i];
3472  }
3473  Assert(new_numbers_copy.size() ==
3474  static_cast<unsigned int>(
3476  tr->get_communicator())]),
3477  ExcInternalError());
3478  ierr = MPI_Allgatherv(new_numbers_copy.data(),
3479  new_numbers_copy.size(),
3481  gathered_new_numbers.data(),
3482  rcounts.data(),
3483  displacements.data(),
3485  tr->get_communicator());
3486  AssertThrowMPI(ierr);
3487  }
3488 
3489  // put new numbers according to the current
3490  // locally_owned_dofs_per_processor IndexSets
3492  // flag_1 and flag_2 are
3493  // used to control that there is a
3494  // one-to-one relation between old and new DoFs.
3495  std::vector<unsigned int> flag_1(this->dof_handler->n_dofs(), 0);
3496  std::vector<unsigned int> flag_2(this->dof_handler->n_dofs(), 0);
3497  for (unsigned int i = 0; i < n_cpu; i++)
3498  {
3499  const IndexSet iset =
3500  this->dof_handler->locally_owned_dofs_per_processor()[i];
3501  for (types::global_dof_index ind = 0; ind < iset.n_elements();
3502  ind++)
3503  {
3504  const types::global_dof_index target =
3505  iset.nth_index_in_set(ind);
3507  gathered_new_numbers[shift + ind];
3508  Assert(target < this->dof_handler->n_dofs(),
3509  ExcInternalError());
3510  Assert(value < this->dof_handler->n_dofs(),
3511  ExcInternalError());
3512  global_gathered_numbers[target] = value;
3513  flag_1[target]++;
3514  flag_2[value]++;
3515  }
3516  shift += iset.n_elements();
3517  }
3518 
3519  Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1,
3520  ExcInternalError());
3521  Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1,
3522  ExcInternalError());
3523  Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1,
3524  ExcInternalError());
3525  Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1,
3526  ExcInternalError());
3527  }
3528 
3529  // let the sequential algorithm do its magic; ignore the
3530  // return type, but reconstruct the number cache based on
3531  // which DoFs each process owns
3532  Implementation::renumber_dofs(global_gathered_numbers,
3533  IndexSet(0),
3534  *this->dof_handler,
3535  /*check_validity=*/true);
3536 
3537  const NumberCache number_cache(
3539  this->dof_handler->get_triangulation().locally_owned_subdomain());
3540 
3541  // restore artificial cells
3542  cell = tr->begin_active();
3543  if (tr->with_artificial_cells())
3544  for (unsigned int index = 0; cell != endc; cell++, index++)
3545  cell->set_subdomain_id(current_subdomain_ids[index]);
3546 
3547  return number_cache;
3548 #endif
3549  }
3550 
3551 
3552 
3553  template <class DoFHandlerType>
3554  NumberCache
3556  const unsigned int /*level*/,
3557  const std::vector<types::global_dof_index> & /*new_numbers*/) const
3558  {
3559  // multigrid is not currently implemented for shared triangulations
3560  Assert(false, ExcNotImplemented());
3561 
3562  return {};
3563  }
3564 
3565 
3566 
3567  /* --------------------- class ParallelDistributed ---------------- */
3568 
3569 #ifdef DEAL_II_WITH_MPI
3570 
3571  namespace
3572  {
3573  template <int dim, int spacedim>
3574  void
3575  get_mg_dofindices_recursively(
3576  const ::parallel::DistributedTriangulationBase<dim, spacedim>
3577  &tria,
3579  & dealii_cell,
3580  const typename CellId::binary_type & quadrant,
3581  std::vector<::types::global_dof_index> &dof_numbers_and_indices)
3582  {
3583  if (dealii_cell->id() == CellId(quadrant))
3584  {
3585  // why would somebody request a cell that is not ours?
3586  Assert(dealii_cell->level_subdomain_id() ==
3587  tria.locally_owned_subdomain(),
3588  ExcInternalError());
3589 
3590  std::vector<::types::global_dof_index> local_dof_indices(
3591  dealii_cell->get_fe().dofs_per_cell);
3592  dealii_cell->get_mg_dof_indices(local_dof_indices);
3593 
3594  dof_numbers_and_indices.push_back(
3595  dealii_cell->get_fe().dofs_per_cell);
3596  dof_numbers_and_indices.insert(dof_numbers_and_indices.end(),
3597  local_dof_indices.begin(),
3598  local_dof_indices.end());
3599  return; // we are done
3600  }
3601 
3602  if (dealii_cell->is_active())
3603  return;
3604 
3605  if (!dealii_cell->id().is_ancestor_of(CellId(quadrant)))
3606  return;
3607 
3608  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
3609  ++c)
3610  get_mg_dofindices_recursively<dim, spacedim>(
3611  tria, dealii_cell->child(c), quadrant, dof_numbers_and_indices);
3612  }
3613 
3614 
3615 
3616  template <int dim, int spacedim>
3617  void
3618  find_marked_mg_ghost_cells_recursively(
3619  const typename ::parallel::
3620  DistributedTriangulationBase<dim, spacedim> &tria,
3621  const unsigned int tree_index,
3623  &dealii_cell,
3624  std::map<
3626  std::vector<std::pair<unsigned int, typename CellId::binary_type>>>
3627  &neighbor_cell_list)
3628  {
3629  // recurse...
3630  if (dealii_cell->has_children())
3631  {
3632  for (unsigned int c = 0;
3633  c < GeometryInfo<dim>::max_children_per_cell;
3634  ++c)
3635  find_marked_mg_ghost_cells_recursively<dim, spacedim>(
3636  tria, tree_index, dealii_cell->child(c), neighbor_cell_list);
3637  }
3638 
3639  if (dealii_cell->user_flag_set() &&
3640  dealii_cell->level_subdomain_id() !=
3641  tria.locally_owned_subdomain())
3642  {
3643  neighbor_cell_list[dealii_cell->level_subdomain_id()]
3644  .emplace_back(tree_index,
3645  dealii_cell->id().template to_binary<spacedim>());
3646  }
3647  }
3648 
3649 
3650 
3651  template <int dim, int spacedim>
3652  void
3653  set_mg_dofindices_recursively(
3654  const ::parallel::DistributedTriangulationBase<dim, spacedim>
3655  &tria,
3657  & dealii_cell,
3658  const typename CellId::binary_type &quadrant,
3659  ::types::global_dof_index * dofs)
3660  {
3661  if (dealii_cell->id() == CellId(quadrant))
3662  {
3663  Assert(dealii_cell->level_subdomain_id() !=
3665  ExcInternalError());
3666 
3667  // update dof indices of cell
3668  std::vector<::types::global_dof_index> dof_indices(
3669  dealii_cell->get_fe().dofs_per_cell);
3670  dealii_cell->get_mg_dof_indices(dof_indices);
3671 
3672  bool complete = true;
3673  for (unsigned int i = 0; i < dof_indices.size(); ++i)
3674  if (dofs[i] != numbers::invalid_dof_index)
3675  {
3676  Assert((dof_indices[i] == (numbers::invalid_dof_index)) ||
3677  (dof_indices[i] == dofs[i]),
3678  ExcInternalError());
3679  dof_indices[i] = dofs[i];
3680  }
3681  else
3682  complete = false;
3683 
3684  if (!complete)
3685  const_cast<
3687  dealii_cell)
3688  ->set_user_flag();
3689  else
3690  const_cast<
3692  dealii_cell)
3693  ->clear_user_flag();
3694 
3695  const_cast<
3697  dealii_cell)
3698  ->set_mg_dof_indices(dof_indices);
3699  return;
3700  }
3701 
3702  if (dealii_cell->is_active())
3703  return;
3704 
3705  if (!dealii_cell->id().is_ancestor_of(CellId(quadrant)))
3706  return;
3707 
3708  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
3709  ++c)
3710  set_mg_dofindices_recursively<dim, spacedim>(tria,
3711  dealii_cell->child(c),
3712  quadrant,
3713  dofs);
3714  }
3715 
3716 
3717 
3718  template <int dim, int spacedim, class DoFHandlerType>
3719  void
3720  communicate_mg_ghost_cells(
3721  const typename ::parallel::
3722  DistributedTriangulationBase<dim, spacedim> &tria,
3723  DoFHandlerType & dof_handler)
3724  {
3725  using QuadrantBufferType =
3726  std::vector<std::pair<unsigned int, typename CellId::binary_type>>;
3727  // build list of cells to request for each neighbor
3728  std::set<::types::subdomain_id> level_ghost_owners =
3729  tria.level_ghost_owners();
3730  std::map<::types::subdomain_id, QuadrantBufferType>
3731  neighbor_cell_list;
3732  for (const auto level_ghost_owner : level_ghost_owners)
3733  neighbor_cell_list[level_ghost_owner] = {};
3734 
3735  for (typename DoFHandlerType::level_cell_iterator cell =
3736  dof_handler.begin(0);
3737  cell != dof_handler.end(0);
3738  ++cell)
3739  {
3741  try
3742  {
3743  coarse_cell_id = cell->id().get_coarse_cell_id();
3744  }
3745  catch (...)
3746  {
3747  // In the case of parallel::fullydistributed::Triangulation,
3748  // a dummy cell throws an exception which is caught here.
3749  // We ignore this cell here.
3750  continue;
3751  };
3752 
3753  find_marked_mg_ghost_cells_recursively<dim, spacedim>(
3754  tria, coarse_cell_id, cell, neighbor_cell_list);
3755  }
3756  Assert(level_ghost_owners.size() == neighbor_cell_list.size(),
3757  ExcInternalError());
3758 
3759 
3760  // Before sending & receiving, make sure we protect this section with
3761  // a mutex:
3762  static Utilities::MPI::CollectiveMutex mutex;
3764  mutex, tria.get_communicator());
3765 
3766  const int mpi_tag = Utilities::MPI::internal::Tags::
3768  const int mpi_tag_reply = Utilities::MPI::internal::Tags::
3770 
3771  //* send our requests:
3772  std::vector<MPI_Request> requests(level_ghost_owners.size());
3773  {
3774  unsigned int idx = 0;
3775  for (const auto &it : neighbor_cell_list)
3776  {
3777  // send the data about the relevant cells
3778  const int ierr =
3779  MPI_Isend(it.second.data(),
3780  it.second.size() * sizeof(it.second[0]),
3781  MPI_BYTE,
3782  it.first,
3783  mpi_tag,
3784  tria.get_communicator(),
3785  &requests[idx]);
3786  AssertThrowMPI(ierr);
3787  ++idx;
3788  }
3789  }
3790 
3791  //* receive requests and reply with the ghost indices
3792  std::vector<QuadrantBufferType> quadrant_data_to_send(
3793  level_ghost_owners.size());
3794  std::vector<std::vector<types::global_dof_index>>
3795  send_dof_numbers_and_indices(level_ghost_owners.size());
3796  std::vector<MPI_Request> reply_requests(level_ghost_owners.size());
3797 
3798  for (unsigned int idx = 0; idx < level_ghost_owners.size(); ++idx)
3799  {
3800  MPI_Status status;
3801  int ierr = MPI_Probe(MPI_ANY_SOURCE,
3802  mpi_tag,
3803  tria.get_communicator(),
3804  &status);
3805  AssertThrowMPI(ierr);
3806 
3807  int len;
3808  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3809  AssertThrowMPI(ierr);
3810  Assert(len % sizeof(quadrant_data_to_send[idx][0]) == 0,
3811  ExcInternalError());
3812 
3813  const unsigned int n_cells =
3814  len / sizeof(quadrant_data_to_send[idx][0]);
3815  quadrant_data_to_send[idx].resize(n_cells);
3816 
3817  ierr = MPI_Recv(quadrant_data_to_send[idx].data(),
3818  len,
3819  MPI_BYTE,
3820  status.MPI_SOURCE,
3821  status.MPI_TAG,
3822  tria.get_communicator(),
3823  &status);
3824  AssertThrowMPI(ierr);
3825 
3826  // store the dof indices for each cell
3827  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells);
3828  ++c)
3829  {
3830  const auto temp =
3831  CellId(quadrant_data_to_send[idx][c].first, 0, nullptr)
3832  .to_cell(tria);
3833 
3834  typename DoFHandlerType::level_cell_iterator cell(
3835  &dof_handler.get_triangulation(),
3836  0,
3837  temp->index(),
3838  &dof_handler);
3839 
3840  get_mg_dofindices_recursively<dim, spacedim>(
3841  tria,
3842  cell,
3843  quadrant_data_to_send[idx][c].second,
3844  send_dof_numbers_and_indices[idx]);
3845  }
3846 
3847  // send reply
3848  ierr = MPI_Isend(send_dof_numbers_and_indices[idx].data(),
3849  send_dof_numbers_and_indices[idx].size(),
3851  status.MPI_SOURCE,
3852  mpi_tag_reply,
3853  tria.get_communicator(),
3854  &reply_requests[idx]);
3855  AssertThrowMPI(ierr);
3856  }
3857 
3858  //* finally receive the replies
3859  for (unsigned int idx = 0; idx < level_ghost_owners.size(); ++idx)
3860  {
3861  MPI_Status status;
3862  int ierr = MPI_Probe(MPI_ANY_SOURCE,
3863  mpi_tag_reply,
3864  tria.get_communicator(),
3865  &status);
3866  AssertThrowMPI(ierr);
3867  int len;
3868  ierr = MPI_Get_count(&status, DEAL_II_DOF_INDEX_MPI_TYPE, &len);
3869  const QuadrantBufferType &quadrants =
3870  neighbor_cell_list[status.MPI_SOURCE];
3871  AssertThrowMPI(ierr);
3872  Assert((len > 0 && !quadrants.empty()) ||
3873  (len == 0 && quadrants.empty()),
3874  ExcInternalError());
3875  std::vector<types::global_dof_index>
3876  receive_dof_numbers_and_indices(len);
3877 
3878  ierr = MPI_Recv(receive_dof_numbers_and_indices.data(),
3879  len,
3881  status.MPI_SOURCE,
3882  status.MPI_TAG,
3883  tria.get_communicator(),
3884  &status);
3885  AssertThrowMPI(ierr);
3886 
3887  // set the dof indices for each cell
3889  receive_dof_numbers_and_indices.data();
3890  for (const auto &it : quadrants)
3891  {
3892  const auto temp = CellId(it.first, 0, nullptr).to_cell(tria);
3893 
3894  typename DoFHandlerType::level_cell_iterator cell(
3895  &tria, 0, temp->index(), &dof_handler);
3896 
3897  Assert(cell->get_fe().dofs_per_cell == dofs[0],
3898  ExcInternalError());
3899 
3900  set_mg_dofindices_recursively<dim, spacedim>(tria,
3901  cell,
3902  it.second,
3903  dofs + 1);
3904  dofs += 1 + dofs[0];
3905  }
3906  Assert(dofs == receive_dof_numbers_and_indices.data() +
3907  receive_dof_numbers_and_indices.size(),
3908  ExcInternalError());
3909  }
3910 
3911  // complete all sends, so that we can safely destroy the
3912  // buffers.
3913  if (requests.size() > 0)
3914  {
3915  const int ierr = MPI_Waitall(requests.size(),
3916  requests.data(),
3917  MPI_STATUSES_IGNORE);
3918  AssertThrowMPI(ierr);
3919  }
3920  if (reply_requests.size() > 0)
3921  {
3922  const int ierr = MPI_Waitall(reply_requests.size(),
3923  reply_requests.data(),
3924  MPI_STATUSES_IGNORE);
3925  AssertThrowMPI(ierr);
3926  }
3927  }
3928 
3929 
3930 
3931  template <int spacedim>
3932  void
3933  communicate_mg_ghost_cells(const typename ::parallel::
3934  distributed::Triangulation<1, spacedim> &,
3936  {
3937  Assert(false, ExcNotImplemented());
3938  }
3939 
3940 
3941 
3942  template <int spacedim>
3943  void
3944  communicate_mg_ghost_cells(const typename ::parallel::
3945  distributed::Triangulation<1, spacedim> &,
3947  {
3948  Assert(false, ExcNotImplemented());
3949  }
3950 
3951 
3952 
3971  template <class DoFHandlerType>
3972  void
3973  communicate_dof_indices_on_marked_cells(
3974  const DoFHandlerType &dof_handler,
3975  const std::map<unsigned int, std::set<::types::subdomain_id>> &)
3976  {
3977 # ifndef DEAL_II_WITH_MPI
3979  Assert(false, ExcNotImplemented());
3980 # else
3981  const unsigned int dim = DoFHandlerType::dimension;
3982  const unsigned int spacedim = DoFHandlerType::space_dimension;
3983 
3984  // define functions that pack data on cells that are ghost cells
3985  // somewhere else, and unpack data on cells where we get information
3986  // from elsewhere
3987  auto pack =
3988  [](const typename DoFHandlerType::active_cell_iterator &cell)
3989  -> std_cxx17::optional<std::vector<types::global_dof_index>> {
3990  Assert(cell->is_locally_owned(), ExcInternalError());
3991 
3992  // first see whether we need to do anything at all on this cell.
3993  // this is determined by whether the user_flag is set on the
3994  // cell that indicates that the *complete* set of DoF indices
3995  // has not been sent
3996  if (cell->user_flag_set())
3997  {
3998  // get dof indices for the current cell
3999  std::vector<types::global_dof_index> local_dof_indices(
4000  cell->get_fe().dofs_per_cell);
4001  cell->get_dof_indices(local_dof_indices);
4002 
4003  // now see if there are dof indices that were previously
4004  // unknown. this can only happen in phase 1, and in
4005  // that case we know that the user flag must have been set
4006  //
4007  // in any case, if the cell *is* complete, we do not
4008  // need to send the data any more in the next phase. indicate
4009  // this by removing the user flag
4010  if (std::find(local_dof_indices.begin(),
4011  local_dof_indices.end(),
4013  local_dof_indices.end())
4014  {
4015  Assert(cell->user_flag_set(), ExcInternalError());
4016  }
4017  else
4018  cell->clear_user_flag();
4019 
4020  return local_dof_indices;
4021  }
4022  else
4023  {
4024  // the fact that the user flag wasn't set means that there is
4025  // nothing we need to send that hasn't been sent so far.
4026  // so return an empty array, but also verify that indeed
4027  // the cell is complete
4028 # ifdef DEBUG
4029  std::vector<types::global_dof_index> local_dof_indices(
4030  cell->get_fe().dofs_per_cell);
4031  cell->get_dof_indices(local_dof_indices);
4032 
4033  const bool is_complete =
4034  (std::find(local_dof_indices.begin(),
4035  local_dof_indices.end(),
4037  local_dof_indices.end());
4038  Assert(is_complete, ExcInternalError());
4039 # endif
4040  return std_cxx17::optional<
4041  std::vector<types::global_dof_index>>();
4042  }
4043  };
4044 
4045  auto unpack =
4046  [](const typename DoFHandlerType::active_cell_iterator &cell,
4047  const std::vector<types::global_dof_index> &received_dof_indices)
4048  -> void {
4049  // this function should only be called on ghost cells, and
4050  // on top of that, only on cells that have not been
4051  // completed -- which we indicate via the user flag.
4052  // check both
4053  Assert(cell->is_ghost(), ExcInternalError());
4054  Assert(cell->user_flag_set(), ExcInternalError());
4055 
4056  // if we just got an incomplete array of DoF indices, then we must
4057  // be in the first ghost exchange and the user flag must have been
4058  // set. we tested that already above.
4059  //
4060  // if we did get a complete array, then we may be in the first
4061  // or second ghost exchange, but in any case we need not exchange
4062  // another time. so delete the user flag
4063  const bool is_complete = (std::find(received_dof_indices.begin(),
4064  received_dof_indices.end(),
4066  received_dof_indices.end());
4067  if (is_complete)
4068  cell->clear_user_flag();
4069 
4070  // in any case, set the DoF indices on this cell. some
4071  // of the ones we received may still be invalid because
4072  // the sending processor did not know them yet, so we
4073  // need to merge the ones we get with those that are
4074  // already set here and may have already been known. for
4075  // those that we already know *and* get, they must obviously
4076  // agree
4077  //
4078  // before getting the local dof indices, we need to update the
4079  // cell dof indices cache because we may have set dof indices
4080  // on a neighboring ghost cell before this one, which may have
4081  // affected the dof indices we know about the current cell
4082  std::vector<types::global_dof_index> local_dof_indices(
4083  cell->get_fe().dofs_per_cell);
4084  cell->update_cell_dof_indices_cache();
4085  cell->get_dof_indices(local_dof_indices);
4086 
4087  for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
4088  if (local_dof_indices[i] == numbers::invalid_dof_index)
4089  local_dof_indices[i] = received_dof_indices[i];
4090  else
4091  // we already know the dof index. check that there
4092  // is no conflict
4093  Assert((received_dof_indices[i] ==
4095  (received_dof_indices[i] == local_dof_indices[i]),
4096  ExcInternalError());
4097 
4098  const_cast<typename DoFHandlerType::active_cell_iterator &>(cell)
4099  ->set_dof_indices(local_dof_indices);
4100  };
4101 
4103  std::vector<types::global_dof_index>,
4104  DoFHandlerType>(dof_handler, pack, unpack);
4105 
4106  // finally update the cell DoF indices caches to make sure
4107  // our internal data structures are consistent
4108  update_all_active_cell_dof_indices_caches(dof_handler);
4109 
4110 
4111  // have a barrier so that sends between two calls to this
4112  // function are not mixed up.
4113  //
4114  // this is necessary because above we just see if there are
4115  // messages and then receive them, without discriminating
4116  // where they come from and whether they were sent in phase
4117  // 1 or 2 (the function is called twice in a row). the need
4118  // for a global communication step like this barrier could
4119  // be avoided by receiving messages specifically from those
4120  // processors from which we expect messages, and by using
4121  // different tags for phase 1 and 2, but the cost of a
4122  // barrier is negligible compared to everything else we do
4123  // here
4124  if (const auto *triangulation =
4125  dynamic_cast<const ::parallel::
4126  DistributedTriangulationBase<dim, spacedim> *>(
4127  &dof_handler.get_triangulation()))
4128  {
4129  const int ierr = MPI_Barrier(triangulation->get_communicator());
4130  AssertThrowMPI(ierr);
4131  }
4132  else
4133  {
4134  Assert(false,
4135  ExcMessage(
4136  "The function communicate_dof_indices_on_marked_cells() "
4137  "only works with parallel distributed triangulations."));
4138  }
4139 # endif
4140  }
4141 
4142 
4143 
4144  } // namespace
4145 
4146 #endif // DEAL_II_WITH_MPI
4147 
4148 
4149 
4150  template <class DoFHandlerType>
4152  DoFHandlerType &dof_handler)
4153  : dof_handler(&dof_handler)
4154  {}
4155 
4156 
4157 
4158  template <class DoFHandlerType>
4159  NumberCache
4161  {
4162 #ifndef DEAL_II_WITH_MPI
4163  Assert(false, ExcNotImplemented());
4164  return NumberCache();
4165 #else
4166  const unsigned int dim = DoFHandlerType::dimension;
4167  const unsigned int spacedim = DoFHandlerType::space_dimension;
4168 
4170  *triangulation =
4171  (dynamic_cast<
4173  const_cast<::Triangulation<dim, spacedim> *>(
4174  &dof_handler->get_triangulation())));
4175  Assert(triangulation != nullptr, ExcInternalError());
4176 
4178  triangulation->locally_owned_subdomain();
4179 
4180 
4181  /*
4182  The following algorithm has a number of stages that are all
4183  documented in the paper that describes the parallel::distributed
4184  functionality:
4185 
4186  1/ locally enumerate dofs on locally owned cells
4187  2/ eliminate dof duplicates on all cells.
4188  un-numerate those that are on interfaces with ghost
4189  cells and that we don't own based on the tie-breaking
4190  criterion. unify dofs afterwards.
4191  3/ unify dofs and re-enumerate the remaining valid ones.
4192  the end result is that we only enumerate locally owned
4193  DoFs
4194  4/ shift indices so that each processor has a unique
4195  range of indices
4196  5/ for all locally owned cells that are ghost
4197  cells somewhere else, send our own DoF indices
4198  to the appropriate set of other processors.
4199  overwrite invalid DoF indices on ghost interfaces
4200  with the corresponding valid ones that we now know.
4201  6/ send DoF indices again to get the correct indices
4202  on ghost cells that we may not have known earlier
4203  */
4204 
4205  // --------- Phase 1: enumerate dofs on locally owned cells
4206  const types::global_dof_index n_initial_local_dofs =
4208 
4209  // --------- Phase 2: eliminate dof duplicates on all cells:
4210  // - un-numerate dofs on interfaces to ghost cells
4211  // that we don't own
4212  // - in case of hp::DoFHandler, unify dofs
4213  std::vector<::types::global_dof_index> renumbering(
4214  n_initial_local_dofs, enumeration_dof_index);
4215 
4216  // first, we invalidate degrees of freedom that belong to processors
4217  // of a lower rank, from which we will receive the final (and lower)
4218  // degrees of freedom later.
4221  renumbering, subdomain_id, *dof_handler);
4222 
4223  // then, we identify DoF duplicates if a hp::DoFHandler is used
4224  std::vector<std::map<types::global_dof_index, types::global_dof_index>>
4225  all_constrained_indices(dim);
4226  Implementation::compute_dof_identities(all_constrained_indices,
4227  *dof_handler);
4228 
4229  // --------- Phase 3: re-enumerate the valid degrees of freedom
4230  // consecutively. thus, we finally receive the
4231  // correct number of locally owned DoFs after
4232  // this step.
4233  //
4234  // the order in which we handle Phases 2 and 3 is important,
4235  // since we want to clarify ownership of degrees of freedom before
4236  // we actually unify and enumerate their indices. otherwise, we could
4237  // end up having a degee of freedom to which only invalid indices will
4238  // be assigned.
4239  const types::global_dof_index n_locally_owned_dofs =
4241  renumbering, all_constrained_indices, *dof_handler);
4242 
4243  // --------- Phase 4: shift indices so that each processor has a unique
4244  // range of indices
4245  ::types::global_dof_index my_shift = 0;
4246  const int ierr =
4247  MPI_Exscan(DEAL_II_MPI_CONST_CAST(&n_locally_owned_dofs),
4248  &my_shift,
4249  1,
4251  MPI_SUM,
4252  triangulation->get_communicator());
4253  AssertThrowMPI(ierr);
4254 
4255  // make dof indices globally consecutive
4256  for (auto &new_index : renumbering)
4257  if (new_index != numbers::invalid_dof_index)
4258  new_index += my_shift;
4259 
4260  // now re-enumerate all dofs to this shifted and condensed
4261  // numbering form. we renumber some dofs as invalid, so
4262  // choose the nocheck-version.
4263  Implementation::renumber_dofs(renumbering,
4264  IndexSet(0),
4265  *dof_handler,
4266  /*check_validity=*/false);
4267 
4268  // now a little bit of housekeeping
4269  const ::types::global_dof_index n_global_dofs =
4270  Utilities::MPI::sum(n_locally_owned_dofs,
4271  triangulation->get_communicator());
4272 
4273  NumberCache number_cache;
4274  number_cache.n_global_dofs = n_global_dofs;
4275  number_cache.n_locally_owned_dofs = n_locally_owned_dofs;
4276  number_cache.locally_owned_dofs = IndexSet(n_global_dofs);
4277  number_cache.locally_owned_dofs.add_range(my_shift,
4278  my_shift +
4279  n_locally_owned_dofs);
4280  number_cache.locally_owned_dofs.compress();
4281 
4282  // this ends the phase where we enumerate degrees of freedom on
4283  // each processor. what is missing is communicating DoF indices
4284  // on ghost cells
4285 
4286  // --------- Phase 5: for all locally owned cells that are ghost
4287  // cells somewhere else, send our own DoF indices
4288  // to the appropriate set of other processors
4289  {
4290  std::vector<bool> user_flags;
4291  triangulation->save_user_flags(user_flags);
4292  triangulation->clear_user_flags();
4293 
4294  // figure out which cells are ghost cells on which we have
4295  // to exchange DoF indices
4296  const std::map<unsigned int, std::set<::types::subdomain_id>>
4299 
4300  // mark all cells that either have to send data (locally
4301  // owned cells that are adjacent to ghost neighbors in some
4302  // way) or receive data (all ghost cells) via the user flags
4303  for (const auto &cell : dof_handler->active_cell_iterators())
4304  if (cell->is_locally_owned())
4305  {
4306  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
4307  if (vertices_with_ghost_neighbors.find(cell->vertex_index(
4308  v)) != vertices_with_ghost_neighbors.end())
4309  {
4310  cell->set_user_flag();
4311  break;
4312  }
4313  }
4314  else if (cell->is_ghost())
4315  cell->set_user_flag();
4316 
4317 
4318 
4319  // Send and receive cells. After this, only the local cells
4320  // are marked, that received new data. This has to be
4321  // communicated in a second communication step.
4322  //
4323  // as explained in the 'distributed' paper, this has to be
4324  // done twice
4325  communicate_dof_indices_on_marked_cells(
4326  *dof_handler, vertices_with_ghost_neighbors);
4327 
4328  // in case of hp::DoFHandlers, we may have received valid
4329  // indices of degrees of freedom that are dominated by a fe
4330  // object adjacent to a ghost interface.
4331  // thus, we overwrite the remaining invalid indices with
4332  // the valid ones in this step.
4334  *dof_handler);
4335 
4336  // --------- Phase 6: all locally owned cells have their correct
4337  // DoF indices set. however, some ghost cells
4338  // may still have invalid ones. thus, exchange
4339  // one more time.
4340  communicate_dof_indices_on_marked_cells(
4341  *dof_handler, vertices_with_ghost_neighbors);
4342 
4343  // at this point, we must have taken care of the data transfer
4344  // on all cells we had previously marked. verify this
4345 # ifdef DEBUG
4346  for (const auto &cell : dof_handler->active_cell_iterators())
4347  Assert(cell->user_flag_set() == false, ExcInternalError());
4348 # endif
4349 
4350  triangulation->load_user_flags(user_flags);
4351  }
4352 
4353 # ifdef DEBUG
4354  // check that we are really done
4355  {
4356  std::vector<::types::global_dof_index> local_dof_indices;
4357 
4358  for (const auto &cell : dof_handler->active_cell_iterators())
4359  if (!cell->is_artificial())
4360  {
4361  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
4362  cell->get_dof_indices(local_dof_indices);
4363  if (local_dof_indices.end() !=
4364  std::find(local_dof_indices.begin(),
4365  local_dof_indices.end(),
4367  {
4368  if (cell->is_ghost())
4369  {
4370  Assert(false,
4371  ExcMessage(
4372  "A ghost cell ended up with incomplete "
4373  "DoF index information. This should not "
4374  "have happened!"));
4375  }
4376  else
4377  {
4378  Assert(
4379  false,
4380  ExcMessage(
4381  "A locally owned cell ended up with incomplete "
4382  "DoF index information. This should not "
4383  "have happened!"));
4384  }
4385  }
4386  }
4387  }
4388 # endif // DEBUG
4389  return number_cache;
4390 #endif // DEAL_II_WITH_MPI
4391  }
4392 
4393 
4394 
4395  template <class DoFHandlerType>
4396  std::vector<NumberCache>
4398  {
4399 #ifndef DEAL_II_WITH_MPI
4400  Assert(false, ExcNotImplemented());
4401  return std::vector<NumberCache>();
4402 #else
4403  const unsigned int dim = DoFHandlerType::dimension;
4404  const unsigned int spacedim = DoFHandlerType::space_dimension;
4405 
4407  *triangulation =
4408  (dynamic_cast<
4410  const_cast<::Triangulation<dim, spacedim> *>(
4411  &dof_handler->get_triangulation())));
4412  Assert(triangulation != nullptr, ExcInternalError());
4413 
4415  ExcMessage(
4416  "Multigrid DoFs can only be distributed on a parallel "
4417  "Triangulation if the flag construct_multigrid_hierarchy "
4418  "is set in the constructor."));
4419 
4420  // loop over all levels that exist globally (across all
4421  // processors), even if the current processor does not in fact
4422  // have any cells on that level or if the local part of the
4423  // Triangulation has fewer levels. we need to do this because
4424  // we need to communicate across all processors on all levels
4425  const unsigned int n_levels = triangulation->n_global_levels();
4426  std::vector<NumberCache> number_caches;
4427  number_caches.reserve(n_levels);
4428  for (unsigned int level = 0; level < n_levels; ++level)
4429  {
4430  NumberCache level_number_cache;
4431 
4432  //* 1. distribute on own subdomain
4433  const unsigned int n_initial_local_dofs =
4435  triangulation->locally_owned_subdomain(), *dof_handler, level);
4436 
4437  //* 2. iterate over ghostcells and kill dofs that are not
4438  // owned by us
4439  std::vector<::types::global_dof_index> renumbering(
4440  n_initial_local_dofs);
4441  for (::types::global_dof_index i = 0; i < renumbering.size();
4442  ++i)
4443  renumbering[i] = i;
4444 
4445  if (level < triangulation->n_levels())
4446  {
4447  std::vector<::types::global_dof_index> local_dof_indices;
4448 
4449  typename DoFHandlerType::level_cell_iterator
4450  cell = dof_handler->begin(level),
4451  endc = dof_handler->end(level);
4452 
4453  for (; cell != endc; ++cell)
4454  if (cell->level_subdomain_id() !=
4456  (cell->level_subdomain_id() <
4457  triangulation->locally_owned_subdomain()))
4458  {
4459  // we found a neighboring ghost cell whose
4460  // subdomain is "stronger" than our own
4461  // subdomain
4462 
4463  // delete all dofs that live there and that we
4464  // have previously assigned a number to
4465  // (i.e. the ones on the interface)
4466  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
4467  cell->get_mg_dof_indices(local_dof_indices);
4468  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
4469  ++i)
4470  if (local_dof_indices[i] != numbers::invalid_dof_index)
4471  renumbering[local_dof_indices[i]] =
4473  }
4474  }
4475 
4476  level_number_cache.n_locally_owned_dofs = 0;
4477  for (types::global_dof_index &index : renumbering)
4478  if (index != numbers::invalid_dof_index)
4479  index = level_number_cache.n_locally_owned_dofs++;
4480 
4481  //* 3. communicate local dofcount and shift ids to make
4482  // them unique
4483  ::types::global_dof_index my_shift = 0;
4484  int ierr = MPI_Exscan(DEAL_II_MPI_CONST_CAST(
4485  &level_number_cache.n_locally_owned_dofs),
4486  &my_shift,
4487  1,
4489  MPI_SUM,
4490  triangulation->get_communicator());
4491  AssertThrowMPI(ierr);
4492 
4493  // The last processor knows about the total number of dofs, so we
4494  // can use a cheaper broadcast rather than an MPI_Allreduce via
4495  // MPI::sum().
4496  level_number_cache.n_global_dofs =
4497  my_shift + level_number_cache.n_locally_owned_dofs;
4498  ierr = MPI_Bcast(&level_number_cache.n_global_dofs,
4499  1,
4502  triangulation->get_communicator()) -
4503  1,
4504  triangulation->get_communicator());
4505 
4506  // shift indices
4507  for (types::global_dof_index &index : renumbering)
4508  if (index != numbers::invalid_dof_index)
4509  index += my_shift;
4510 
4511  // now re-enumerate all dofs to this shifted and condensed
4512  // numbering form. we renumber some dofs as invalid, so
4513  // choose the nocheck-version of the function
4514  //
4515  // of course there is nothing for us to renumber if the
4516  // level we are currently dealing with doesn't even exist
4517  // within the current triangulation, so skip renumbering
4518  // in that case
4519  if (level < triangulation->n_levels())
4521  renumbering, IndexSet(0), *dof_handler, level, false);
4522 
4523  // now a little bit of housekeeping
4524  level_number_cache.locally_owned_dofs =
4525  IndexSet(level_number_cache.n_global_dofs);
4526  level_number_cache.locally_owned_dofs.add_range(
4527  my_shift, my_shift + level_number_cache.n_locally_owned_dofs);
4528  level_number_cache.locally_owned_dofs.compress();
4529 
4530  number_caches.emplace_back(level_number_cache);
4531  }
4532 
4533 
4534  //* communicate ghost DoFs
4535  // We mark all ghost cells by setting the user_flag and then request
4536  // these cells from the corresponding owners. As this information
4537  // can be incomplete,
4538  {
4539  std::vector<bool> user_flags;
4540  triangulation->save_user_flags(user_flags);
4541  triangulation->clear_user_flags();
4542 
4543  // mark all ghost cells for transfer
4544  {
4545  typename DoFHandlerType::level_cell_iterator cell,
4546  endc = dof_handler->end();
4547  for (cell = dof_handler->begin(); cell != endc; ++cell)
4548  if (cell->level_subdomain_id() !=
4550  !cell->is_locally_owned_on_level())
4551  cell->set_user_flag();
4552  }
4553 
4554  // Phase 1. Request all marked cells from corresponding owners. If we
4555  // managed to get every DoF, remove the user_flag, otherwise we
4556  // will request them again in the step below.
4557  communicate_mg_ghost_cells(*triangulation, *dof_handler);
4558 
4559  // have a barrier so that sends from above and below this
4560  // place are not mixed up.
4561  //
4562  // this is necessary because above we just see if there are
4563  // messages and then receive them, without discriminating
4564  // where they come from and whether they were sent in phase
4565  // 1 or 2 in communicate_mg_ghost_cells() on another
4566  // processor. the need for a global communication step like
4567  // this barrier could be avoided by receiving messages
4568  // specifically from those processors from which we expect
4569  // messages, and by using different tags for phase 1 and 2,
4570  // but the cost of a barrier is negligible compared to
4571  // everything else we do here
4572  const int ierr = MPI_Barrier(triangulation->get_communicator());
4573  AssertThrowMPI(ierr);
4574 
4575  // Phase 2, only request the cells that were not completed
4576  // in Phase 1.
4577  communicate_mg_ghost_cells(*triangulation, *dof_handler);
4578 
4579 # ifdef DEBUG
4580  // make sure we have removed all flags:
4581  {
4582  typename DoFHandlerType::level_cell_iterator cell,
4583  endc = dof_handler->end();
4584  for (cell = dof_handler->begin(); cell != endc; ++cell)
4585  if (cell->level_subdomain_id() !=
4587  !cell->is_locally_owned_on_level())
4588  Assert(cell->user_flag_set() == false, ExcInternalError());
4589  }
4590 # endif
4591 
4592  triangulation->load_user_flags(user_flags);
4593  }
4594 
4595 
4596 
4597 # ifdef DEBUG
4598  // check that we are really done
4599  {
4600  std::vector<::types::global_dof_index> local_dof_indices;
4601  typename DoFHandlerType::level_cell_iterator cell,
4602  endc = dof_handler->end();
4603 
4604  for (cell = dof_handler->begin(); cell != endc; ++cell)
4605  if (cell->level_subdomain_id() !=
4607  {
4608  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
4609  cell->get_mg_dof_indices(local_dof_indices);
4610  if (local_dof_indices.end() !=
4611  std::find(local_dof_indices.begin(),
4612  local_dof_indices.end(),
4614  {
4615  Assert(false, ExcMessage("not all DoFs got distributed!"));
4616  }
4617  }
4618  }
4619 # endif // DEBUG
4620 
4621  return number_caches;
4622 
4623 #endif // DEAL_II_WITH_MPI
4624  }
4625 
4626 
4627  template <class DoFHandlerType>
4628  NumberCache
4630  const std::vector<::types::global_dof_index> &new_numbers) const
4631  {
4632  (void)new_numbers;
4633 
4634  Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
4635  ExcInternalError());
4636 
4637 #ifndef DEAL_II_WITH_MPI
4638  Assert(false, ExcNotImplemented());
4639  return NumberCache();
4640 #else
4641  const unsigned int dim = DoFHandlerType::dimension;
4642  const unsigned int spacedim = DoFHandlerType::space_dimension;
4643 
4645  *triangulation =
4646  (dynamic_cast<
4648  const_cast<::Triangulation<dim, spacedim> *>(
4649  &dof_handler->get_triangulation())));
4650  Assert(triangulation != nullptr, ExcInternalError());
4651 
4652 
4653  // We start by checking whether only the numbering within the MPI
4654  // ranks changed. In that case, we can apply the renumbering with some
4655  // local renumbering only (this is similar to the renumber_mg_dofs()
4656  // function below)
4657  const bool locally_owned_set_changes =
4658  std::any_of(new_numbers.cbegin(),
4659  new_numbers.cend(),
4660  [this](const types::global_dof_index i) {
4661  return dof_handler->locally_owned_dofs().is_element(
4662  i) == false;
4663  });
4664 
4665  if (Utilities::MPI::sum(static_cast<unsigned int>(
4666  locally_owned_set_changes),
4667  triangulation->get_communicator()) == 0)
4668  {
4669  // Since only the order within the local subdomains has changed,
4670  // all we need to do is to propagate the knowledge about the
4671  // numbers from the locally owned dofs (given by the new_numbers
4672  // array) to all ghosted dofs on neighboring processors. We can do
4673  // this by ghost layer exchange routines as in parallel vectors:
4674  // We create an IndexSet for the relevant dofs and then export
4675  // into an array of those values via Utilities::MPI::Partitioner.
4676  IndexSet relevant_dofs;
4678  relevant_dofs);
4679  std::vector<types::global_dof_index> ghosted_new_numbers(
4680  relevant_dofs.n_elements());
4681  {
4682  Utilities::MPI::Partitioner partitioner(
4683  dof_handler->locally_owned_dofs(),
4684  relevant_dofs,
4685  triangulation->get_communicator());
4686 
4687  // choose some number that makes it unlikely to get conflicts
4688  // with other ongoing non-blocking communication (there
4689  // shouldn't be any at this place in most programs).
4690  const unsigned int communication_channel = 19;
4691  std::vector<types::global_dof_index> temp_array(
4692  partitioner.n_import_indices());
4693  std::vector<MPI_Request> requests;
4694  partitioner.export_to_ghosted_array_start(
4695  communication_channel,
4696  make_array_view(new_numbers),
4697  make_array_view(temp_array),
4699  ghosted_new_numbers.data() + new_numbers.size(),
4700  partitioner.n_ghost_indices()),
4701  requests);
4702  partitioner.export_to_ghosted_array_finish(
4704  ghosted_new_numbers.data() + new_numbers.size(),
4705  partitioner.n_ghost_indices()),
4706  requests);
4707 
4708  // we need to fill the indices of the locally owned part into
4709  // the new numbers array, which is not provided by the parallel
4710  // partitioner. their right position is somewhere in the middle
4711  // of the array, so we first copy the ghosted part from smaller
4712  // ranks to the front, then insert the data in the middle.
4713  unsigned int n_ghosts_on_smaller_ranks = 0;
4714  for (std::pair<unsigned int, unsigned int> t :
4715  partitioner.ghost_targets())
4716  {
4717  if (t.first > partitioner.this_mpi_process())
4718  break;
4719  n_ghosts_on_smaller_ranks += t.second;
4720  }
4721  if (n_ghosts_on_smaller_ranks > 0)
4722  {
4723  Assert(ghosted_new_numbers.data() != nullptr,
4724  ExcInternalError());
4725  std::memmove(ghosted_new_numbers.data(),
4726  ghosted_new_numbers.data() + new_numbers.size(),
4727  sizeof(types::global_dof_index) *
4728  n_ghosts_on_smaller_ranks);
4729  }
4730  if (new_numbers.size() > 0)
4731  {
4732  Assert(new_numbers.data() != nullptr, ExcInternalError());
4733  std::memcpy(ghosted_new_numbers.data() +
4734  n_ghosts_on_smaller_ranks,
4735  new_numbers.data(),
4736  sizeof(types::global_dof_index) *
4737  new_numbers.size());
4738  }
4739  }
4740 
4741  // In case we do not carry any relevant dof (but only some remote
4742  // processor), we do not need to call the renumbering. We call the
4743  // version without validity check because vertex dofs will be
4744  // set already in the artificial region.
4745  if (relevant_dofs.n_elements() > 0)
4746  Implementation::renumber_dofs(ghosted_new_numbers,
4747  relevant_dofs,
4748  *dof_handler,
4749  /*check_validity=*/false);
4750 
4751  NumberCache number_cache;
4752  number_cache.locally_owned_dofs = dof_handler->locally_owned_dofs();
4753  number_cache.n_global_dofs = dof_handler->n_dofs();
4754  number_cache.n_locally_owned_dofs =
4755  number_cache.locally_owned_dofs.n_elements();
4756  return number_cache;
4757  }
4758  else
4759  {
4760  // Now back to the more complicated case
4761  //
4762  // First figure out the new set of locally owned DoF indices.
4763  // If we own no DoFs, we still need to go through this function,
4764  // but we can skip this calculation.
4765  //
4766  // The IndexSet::add_indices() function is substantially more
4767  // efficient if the set of indices is already sorted because
4768  // it can then insert ranges instead of individual elements.
4769  // consequently, pre-sort the array of new indices
4770  IndexSet my_locally_owned_new_dof_indices(dof_handler->n_dofs());
4771  if (dof_handler->n_locally_owned_dofs() > 0)
4772  {
4773  std::vector<::types::global_dof_index>
4774  new_numbers_sorted = new_numbers;
4775  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4776 
4777  my_locally_owned_new_dof_indices.add_indices(
4778  new_numbers_sorted.begin(), new_numbers_sorted.end());
4779  my_locally_owned_new_dof_indices.compress();
4780 
4781  Assert(my_locally_owned_new_dof_indices.n_elements() ==
4782  new_numbers.size(),
4783  ExcInternalError());
4784  }
4785 
4786  // delete all knowledge of DoF indices that are not locally
4787  // owned. we do so by getting DoF indices on cells, checking
4788  // whether they are locally owned, if not, setting them to
4789  // an invalid value, and then setting them again on the current
4790  // cell
4791  //
4792  // DoFs we (i) know about, and (ii) don't own locally must be
4793  // located either on ghost cells, or on the interface between a
4794  // locally owned cell and a ghost cell. In any case, it is
4795  // sufficient to kill them only from the ghost side cell, so loop
4796  // only over ghost cells
4797  {
4798  std::vector<::types::global_dof_index> local_dof_indices;
4799 
4800  for (auto cell : dof_handler->active_cell_iterators())
4801  if (cell->is_ghost())
4802  {
4803  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
4804  cell->get_dof_indices(local_dof_indices);
4805 
4806  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
4807  ++i)
4808  // delete a DoF index if it has not already been deleted
4809  // (e.g., by visiting a neighboring cell, if it is on the
4810  // boundary), and if we don't own it
4811  if ((local_dof_indices[i] !=
4813  (!dof_handler->locally_owned_dofs().is_element(
4814  local_dof_indices[i])))
4815  local_dof_indices[i] = numbers::invalid_dof_index;
4816 
4817  cell->set_dof_indices(local_dof_indices);
4818  }
4819  }
4820 
4821 
4822  // renumber. Skip when there is nothing to do because we own no DoF.
4823  if (dof_handler->locally_owned_dofs().n_elements() > 0)
4824  Implementation::renumber_dofs(new_numbers,
4825  dof_handler->locally_owned_dofs(),
4826  *dof_handler,
4827  /*check_validity=*/false);
4828 
4829  // Communicate newly assigned DoF indices to other processors
4830  // and get the same information for our own ghost cells.
4831  //
4832  // This is the same as phase 5+6 in the distribute_dofs() algorithm,
4833  // taking into account that we have to unify a few DoFs in between
4834  // then communication phases if we do hp numbering
4835  {
4836  std::vector<bool> user_flags;
4837  triangulation->save_user_flags(user_flags);
4838  triangulation->clear_user_flags();
4839 
4840  // mark all own cells for transfer
4841  for (const auto &cell : dof_handler->active_cell_iterators())
4842  if (!cell->is_artificial())
4843  cell->set_user_flag();
4844 
4845  // figure out which cells are ghost cells on which we have
4846  // to exchange DoF indices
4847  const std::map<unsigned int,
4848  std::set<::types::subdomain_id>>
4851  *triangulation);
4852 
4853 
4854  // Send and receive cells. After this, only the local cells
4855  // are marked, that received new data. This has to be
4856  // communicated in a second communication step.
4857  //
4858  // as explained in the 'distributed' paper, this has to be
4859  // done twice
4860  communicate_dof_indices_on_marked_cells(
4861  *dof_handler, vertices_with_ghost_neighbors);
4862 
4863  // in case of hp::DoFHandlers, we may have received valid
4864  // indices of degrees of freedom that are dominated by a fe
4865  // object adjacent to a ghost interface.
4866  // thus, we overwrite the remaining invalid indices with
4867  // the valid ones in this step.
4869  *dof_handler);
4870 
4871  communicate_dof_indices_on_marked_cells(
4872  *dof_handler, vertices_with_ghost_neighbors);
4873 
4874  triangulation->load_user_flags(user_flags);
4875  }
4876 
4877  NumberCache number_cache;
4878  number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
4879  number_cache.n_global_dofs = dof_handler->n_dofs();
4880  number_cache.n_locally_owned_dofs =
4881  number_cache.locally_owned_dofs.n_elements();
4882  return number_cache;
4883  }
4884 #endif
4885  }
4886 
4887 
4888 
4889  template <class DoFHandlerType>
4890  NumberCache
4892  const unsigned int level,
4893  const std::vector<types::global_dof_index> &new_numbers) const
4894  {
4895  // we only implement the case where the multigrid numbers are
4896  // renumbered within the processor's partition, rather than the most
4897  // general case
4898  const IndexSet index_set = dof_handler->locally_owned_mg_dofs(level);
4899 
4900 #ifdef DEAL_II_WITH_MPI
4901 
4902  constexpr int dim = DoFHandlerType::dimension;
4903  constexpr int spacedim = DoFHandlerType::space_dimension;
4904  const ::parallel::TriangulationBase<dim, spacedim> *tr =
4905  (dynamic_cast<const ::parallel::TriangulationBase<dim, spacedim>
4906  *>(&this->dof_handler->get_triangulation()));
4907  Assert(tr != nullptr, ExcInternalError());
4908 
4909  const unsigned int my_rank =
4910  Utilities::MPI::this_mpi_process(tr->get_communicator());
4911 
4912 # ifdef DEBUG
4913  for (types::global_dof_index i : new_numbers)
4914  {
4915  Assert(index_set.is_element(i),
4917  "Renumberings that change the locally owned mg dofs "
4918  "partitioning are currently not implemented for "
4919  "the multigrid levels"));
4920  }
4921 # endif
4922 
4923  // we need to access all locally relevant degrees of freedom. we
4924  // use Utilities::MPI::Partitioner for handling the data exchange
4925  // of the new numbers, which is simply the extraction of ghost data
4926  IndexSet relevant_dofs;
4928  level,
4929  relevant_dofs);
4930  std::vector<types::global_dof_index> ghosted_new_numbers(
4931  relevant_dofs.n_elements());
4932  {
4933  Utilities::MPI::Partitioner partitioner(index_set,
4934  relevant_dofs,
4935  tr->get_communicator());
4936  std::vector<types::global_dof_index> temp_array(
4937  partitioner.n_import_indices());
4938  const unsigned int communication_channel = 17;
4939  std::vector<MPI_Request> requests;
4940  partitioner.export_to_ghosted_array_start(
4941  communication_channel,
4942  make_array_view(new_numbers),
4943  make_array_view(temp_array),
4944  ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
4945  new_numbers.size(),
4946  partitioner.n_ghost_indices()),
4947  requests);
4948  partitioner.export_to_ghosted_array_finish(
4949  ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
4950  new_numbers.size(),
4951  partitioner.n_ghost_indices()),
4952  requests);
4953 
4954  // we need to fill the indices of the locally owned part into the
4955  // new numbers array. their right position is somewhere in the
4956  // middle of the array, so we first copy the ghosted part from
4957  // smaller ranks to the front, then insert the data in the middle.
4958  unsigned int n_ghosts_on_smaller_ranks = 0;
4959  for (std::pair<unsigned int, unsigned int> t :
4960  partitioner.ghost_targets())
4961  {
4962  if (t.first > my_rank)
4963  break;
4964  n_ghosts_on_smaller_ranks += t.second;
4965  }
4966  if (n_ghosts_on_smaller_ranks > 0)
4967  {
4968  Assert(ghosted_new_numbers.data() != nullptr, ExcInternalError());
4969  std::memmove(ghosted_new_numbers.data(),
4970  ghosted_new_numbers.data() + new_numbers.size(),
4971  sizeof(types::global_dof_index) *
4972  n_ghosts_on_smaller_ranks);
4973  }
4974  if (new_numbers.size() > 0)
4975  {
4976  Assert(new_numbers.data() != nullptr, ExcInternalError());
4977  std::memcpy(ghosted_new_numbers.data() +
4978  n_ghosts_on_smaller_ranks,
4979  new_numbers.data(),
4980  sizeof(types::global_dof_index) * new_numbers.size());
4981  }
4982  }
4983 
4984  // in case we do not own any of the given level (but only some remote
4985  // processor), we do not need to call the renumbering
4986  if (level < this->dof_handler->get_triangulation().n_levels() &&
4987  relevant_dofs.n_elements() > 0)
4989  ghosted_new_numbers, relevant_dofs, *dof_handler, level, true);
4990 #else
4991  (void)new_numbers;
4992  Assert(false, ExcNotImplemented());
4993 #endif
4994 
4995  NumberCache number_cache;
4996  number_cache.locally_owned_dofs = index_set;
4997  number_cache.n_global_dofs = dof_handler->n_dofs();
4998  number_cache.n_locally_owned_dofs =
4999  number_cache.locally_owned_dofs.n_elements();
5000  return number_cache;
5001  }
5002  } // namespace Policy
5003  } // namespace DoFHandlerImplementation
5004 } // namespace internal
5005 
5006 
5007 
5008 /*-------------- Explicit Instantiations -------------------------------*/
5009 #include "dof_handler_policy.inst"
5010 
5011 
void get_active_fe_indices(std::vector< unsigned int > &active_fe_indices) const
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1329
static void renumber_vertex_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
const Triangulation< dim, spacedim > & get_triangulation() const
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
static const unsigned int invalid_unsigned_int
Definition: types.h:191
static void renumber_face_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 2, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandlerType &dof_handler)
Definition: dof_tools.cc:1375
const types::subdomain_id invalid_subdomain_id
Definition: types.h:285
void clear_user_flags()
Definition: tria.cc:11390
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:918
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 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 DoFHandlerType &)
Triangulation< dim, spacedim >::cell_iterator to_cell(const Triangulation< dim, spacedim > &tria) const
Definition: cell_id.cc:159
static types::global_dof_index distribute_dofs(const types::subdomain_id subdomain_id, DoFHandlerType &dof_handler)
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 3, spacedim > &dof_handler)
Task< RT > new_task(const std::function< RT()> &function)
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, hp::DoFHandler< 3, spacedim > &dof_handler)
void load_user_flags(std::istream &in)
Definition: tria.cc:11450
dof_handler_policy.cc: communicate_mg_ghost_cells()
Definition: mpi_tags.h:66
cell_iterator end() const
Definition: dof_handler.cc:947
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:933
static types::global_dof_index unify_dof_indices(const DoFHandler< dim, spacedim > &, const unsigned int n_dofs_before_identification, const bool)
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, hp::DoFHandler< 2, spacedim > &dof_handler)
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition: types.h:86
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const hp::DoFHandler< 3, spacedim > &dof_handler)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1520
Point< 2 > second
Definition: grid_out.cc:4353
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 > &, const IndexSet &, hp::DoFHandler< 1, spacedim > &)
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:79
static void renumber_vertex_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, hp::DoFHandler< dim, spacedim > &dof_handler, const bool check_validity)
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition: dof_handler.h:1341
static void merge_invalid_line_dofs_on_ghost_interfaces(hp::DoFHandler< dim, spacedim > &dof_handler)
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
void extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1215
static void merge_invalid_dof_indices_on_ghost_interfaces(const DoFHandler< dim, spacedim > &)
static types::global_dof_index unify_dof_indices(hp::DoFHandler< dim, spacedim > &dof_handler, const unsigned int n_dofs_before_identification, const bool check_validity)
static void renumber_cell_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, hp::DoFHandler< dim, spacedim > &dof_handler)
virtual std::vector< NumberCache > distribute_mg_dofs() const override
static void merge_invalid_quad_dofs_on_ghost_interfaces(hp::DoFHandler< 3, spacedim > &dof_handler)
static void renumber_face_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &)
size_type size() const
Definition: index_set.h:1635
virtual std::vector< NumberCache > distribute_mg_dofs() const override
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::quad_iterator quad_iterator
Definition: dof_handler.h:270
virtual std::vector< NumberCache > distribute_mg_dofs() const override
virtual bool is_multilevel_hierarchy_constructed() const =0
unsigned int subdomain_id
Definition: types.h:43
static void merge_invalid_quad_dofs_on_ghost_interfaces(hp::DoFHandler< dim, spacedim > &)
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1408
IteratorRange< active_cell_iterator > active_cell_iterators() const
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1922
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
static void compute_dof_identities(const std::vector< std::map< types::global_dof_index, types::global_dof_index >> &, const DoFHandler< dim, spacedim > &)
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:607
unsigned int max_dofs_per_cell(const DoFHandler< dim, spacedim > &dh)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:361
unsigned int level
Definition: grid_out.cc:4355
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
static void renumber_cell_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler)
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:393
void set_dof_index(const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1173
types::subdomain_id locally_owned_subdomain() const override
Definition: tria_base.cc:288
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:911
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const hp::DoFHandler< dim, spacedim > &dof_handler)
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > faces
Definition: dof_handler.h:1349
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1209
const hp::FECollection< dim, spacedim > & get_fe_collection() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< typename ::internal::p4est::types< dim >::quadrant > quadrants
Definition: tria.cc:2241
static void merge_invalid_line_dofs_on_ghost_interfaces(hp::DoFHandler< 1, spacedim > &)
Definition: hp.h:117
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > levels
Definition: dof_handler.h:1337
std::vector< unsigned int > tree_index
Definition: tria.cc:2237
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12834
static std::map< types::global_dof_index, types::global_dof_index > compute_vertex_dof_identities(const hp::DoFHandler< dim, spacedim > &dof_handler)
static void merge_invalid_vertex_dofs_on_ghost_interfaces(hp::DoFHandler< dim, spacedim > &dof_handler)
virtual const MPI_Comm & get_communicator() const
Definition: tria_base.cc:138
Definition: cell_id.h:69
static void renumber_mg_dofs(const std::vector<::types::global_dof_index > &, const IndexSet &, hp::DoFHandler< dim, spacedim > &, const unsigned int, const bool)
Point< 2 > first
Definition: grid_out.cc:4352
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const hp::DoFHandler< 1, spacedim > &)
static void compute_dof_identities(std::vector< std::map< types::global_dof_index, types::global_dof_index >> &all_constrained_indices, const hp::DoFHandler< dim, spacedim > &dof_handler)
#define DEAL_II_MPI_CONST_CAST(expr)
Definition: mpi.h:76
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1674
unsigned int global_dof_index
Definition: types.h:76
const types::subdomain_id artificial_subdomain_id
Definition: types.h:302
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
void compress() const
Definition: index_set.h:1643
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:922
static void renumber_face_mg_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &, const unsigned int, const bool)
dof_handler_policy.cc: communicate_mg_ghost_cells()
Definition: mpi_tags.h:69
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1696
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1051
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:360
virtual unsigned int n_global_levels() const override
Definition: tria_base.cc:124
const Triangulation< dim, spacedim > & get_triangulation() const
void export_to_ghosted_array_start(const unsigned int communication_channel, const ArrayView< const Number, MemorySpaceType > &locally_owned_array, const ArrayView< Number, MemorySpaceType > &temporary_storage, const ArrayView< Number, MemorySpaceType > &ghost_array, std::vector< MPI_Request > &requests) const
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1353
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
void save_user_flags(std::ostream &out) const
Definition: tria.cc:11401
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index) const
static ::ExceptionBase & ExcNotImplemented()
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1184
static unsigned int n_threads()
Definition: table.h:695
bool is_element(const size_type index) const
Definition: index_set.h:1766
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const hp::DoFHandler< dim, spacedim > &)
static types::global_dof_index distribute_dofs_on_level(const types::subdomain_id level_subdomain_id, DoFHandlerType &dof_handler, const unsigned int level)
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:357
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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 DoFHandlerType &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_face_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 3, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
static void renumber_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
const types::global_dof_index invalid_dof_index
Definition: types.h:206
std::map< unsigned int, std::set<::types::subdomain_id > > * vertices_with_ghost_neighbors
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)
static void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandlerType &dof_handler, const bool check_validity)
global_cell_index coarse_cell_id
Definition: types.h:114
static const bool value
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
size_type n_elements() const
Definition: index_set.h:1833
static void merge_invalid_dof_indices_on_ghost_interfaces(hp::DoFHandler< dim, spacedim > &dof_handler)
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1881
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:817
std::vector< types::global_dof_index > vertex_dofs
Definition: dof_handler.h:1323
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
typename ActiveSelector::line_iterator line_iterator
Definition: dof_handler.h:246
static ::ExceptionBase & ExcInternalError()
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:5530