Reference documentation for deal.II version 9.2.0
\(\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 
17 #include <deal.II/base/geometry_info.h>
18 #include <deal.II/base/memory_consumption.h>
19 #include <deal.II/base/partitioner.h>
20 #include <deal.II/base/thread_management.h>
21 #include <deal.II/base/utilities.h>
22 #include <deal.II/base/work_stream.h>
23 
24 #include <deal.II/distributed/shared_tria.h>
25 #include <deal.II/distributed/tria.h>
26 
27 #include <deal.II/dofs/dof_accessor.h>
28 #include <deal.II/dofs/dof_handler.h>
29 #include <deal.II/dofs/dof_handler_policy.h>
30 
31 #include <deal.II/fe/fe.h>
32 
33 #include <deal.II/grid/grid_tools.h>
34 #include <deal.II/grid/tria.h>
35 #include <deal.II/grid/tria_iterator.h>
36 
37 #include <algorithm>
38 #include <memory>
39 #include <numeric>
40 #include <set>
41 
42 DEAL_II_NAMESPACE_OPEN
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 =
71  numbers::invalid_dof_index - 1;
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,
106  2 * MultithreadInfo::n_threads(),
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,
144  2 * MultithreadInfo::n_threads(),
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_cxx14::make_unique<DoFIdentities>(
179  fe1.hp_vertex_dof_identities(fe2));
180  break;
181  }
182 
183  case 1:
184  {
185  identities = std_cxx14::make_unique<DoFIdentities>(
186  fe1.hp_line_dof_identities(fe2));
187  break;
188  }
189 
190  case 2:
191  {
192  identities = std_cxx14::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 
218  struct Implementation
219  {
220  /* -------------- distribute_dofs functionality ------------- */
221 
226  template <int dim, int spacedim>
227  static std::map<types::global_dof_index, types::global_dof_index>
228  compute_vertex_dof_identities(
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...
240  ::Table<2, std::unique_ptr<DoFIdentities>>
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 ==
272  numbers::invalid_unsigned_int)
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 !=
344  numbers::invalid_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>
370  compute_line_dof_identities(const hp::DoFHandler<1, spacedim> &)
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>
378  compute_line_dof_identities(
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  {
418  const typename hp::DoFHandler<dim, spacedim>::line_iterator
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 =
503  numbers::invalid_unsigned_int;
504 
505  if (dominating_fe_index !=
506  numbers::invalid_unsigned_int)
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  {
523  const types::global_dof_index
524  master_dof_index = line->dof_index(
525  j, dominating_fe_index);
526  const types::global_dof_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 !=
545  numbers::invalid_dof_index)
546  {
547  if (master_dof_index !=
548  numbers::invalid_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] =
607  numbers::invalid_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 !=
644  numbers::invalid_unsigned_int)
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  {
663  const types::global_dof_index
664  master_dof_index = line->dof_index(
665  identity.first,
666  most_dominating_fe_index);
667  const types::global_dof_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 !=
686  numbers::invalid_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 !=
696  numbers::invalid_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>
730  compute_quad_dof_identities(const hp::DoFHandler<dim, spacedim> &)
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>
743  compute_quad_dof_identities(
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  {
780  const typename hp::DoFHandler<dim, spacedim>::quad_iterator
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 !=
833  numbers::invalid_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 !=
843  numbers::invalid_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
876  compute_dof_identities(
877  const std::vector<
878  std::map<types::global_dof_index, types::global_dof_index>> &,
879  const DoFHandler<dim, spacedim> &)
880  {}
881 
882 
883  template <int dim, int spacedim>
884  static void
885  compute_dof_identities(std::vector<std::map<types::global_dof_index,
886  types::global_dof_index>>
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>
943  static types::global_dof_index
944  enumerate_dof_indices_for_renumbering(
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>
1009  static types::global_dof_index
1010  unify_dof_indices(const DoFHandler<dim, 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>
1019  static types::global_dof_index
1020  unify_dof_indices(hp::DoFHandler<dim, spacedim> &dof_handler,
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
1051  merge_invalid_vertex_dofs_on_ghost_interfaces(
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...
1060  ::Table<2, std::unique_ptr<DoFIdentities>>
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 !=
1108  numbers::invalid_unsigned_int)
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  {
1138  const types::global_dof_index
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);
1146  const types::global_dof_index
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 ==
1168  numbers::invalid_dof_index) &&
1169  (master_dof_index !=
1170  numbers::invalid_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>
1193  static void merge_invalid_line_dofs_on_ghost_interfaces(
1194  hp::DoFHandler<1, spacedim> &)
1195  {}
1196 
1197 
1198  template <int dim, int spacedim>
1199  static void
1200  merge_invalid_line_dofs_on_ghost_interfaces(
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  {
1245  const typename hp::DoFHandler<dim, spacedim>::line_iterator
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 =
1314  numbers::invalid_unsigned_int;
1315 
1316  if (dominating_fe_index !=
1317  numbers::invalid_unsigned_int)
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  {
1334  const types::global_dof_index
1335  master_dof_index = line->dof_index(
1336  j, dominating_fe_index);
1337  const types::global_dof_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 ==
1355  numbers::invalid_dof_index) &&
1356  (master_dof_index !=
1357  numbers::invalid_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 !=
1395  numbers::invalid_unsigned_int)
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  {
1414  const types::global_dof_index
1415  master_dof_index = line->dof_index(
1416  identity.first,
1417  most_dominating_fe_index);
1418  const types::global_dof_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 ==
1436  numbers::invalid_dof_index) &&
1437  (master_dof_index !=
1438  numbers::invalid_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
1462  merge_invalid_quad_dofs_on_ghost_interfaces(
1463  hp::DoFHandler<dim, spacedim> &)
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>
1473  static void merge_invalid_quad_dofs_on_ghost_interfaces(
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  {
1514  const typename hp::DoFHandler<dim, spacedim>::quad_iterator
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 ==
1573  numbers::invalid_dof_index) &&
1574  (master_dof_index !=
1575  numbers::invalid_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
1606  merge_invalid_dof_indices_on_ghost_interfaces(
1607  const DoFHandler<dim, spacedim> &)
1608  {}
1609 
1610 
1611  template <int dim, int spacedim>
1612  static void
1613  merge_invalid_dof_indices_on_ghost_interfaces(
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>
1652  static types::global_dof_index
1653  distribute_dofs(const types::subdomain_id subdomain_id,
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  typename DoFHandlerType::active_cell_iterator
1663  cell = dof_handler.begin_active(),
1664  endc = dof_handler.end();
1665 
1666  for (; cell != endc; ++cell)
1667  if (!cell->is_artificial())
1668  if ((subdomain_id == numbers::invalid_subdomain_id) ||
1669  (cell->subdomain_id() == subdomain_id))
1670  {
1671  std::vector<types::global_dof_index> dof_indices(
1672  cell->get_fe().dofs_per_cell);
1673 
1674  internal::DoFAccessorImplementation::get_dof_indices(
1675  *cell, dof_indices, cell->active_fe_index());
1676 
1677  for (auto &dof_index : dof_indices)
1678  if (dof_index == numbers::invalid_dof_index)
1679  dof_index = next_free_dof++;
1680 
1681  cell->set_dof_indices(dof_indices);
1682  }
1683 
1684  update_all_active_cell_dof_indices_caches(dof_handler);
1685 
1686  return next_free_dof;
1687  }
1688 
1689 
1690 
1704  template <class DoFHandlerType>
1705  static void
1706  invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering(
1707  std::vector<types::global_dof_index> &renumbering,
1708  const types::subdomain_id subdomain_id,
1709  const DoFHandlerType & dof_handler)
1710  {
1711  std::vector<types::global_dof_index> local_dof_indices;
1712 
1713  for (const auto &cell : dof_handler.active_cell_iterators())
1714  if (cell->is_ghost() && (cell->subdomain_id() < subdomain_id))
1715  {
1716  // we found a neighboring ghost cell whose subdomain
1717  // is "stronger" than our own subdomain
1718 
1719  // delete all dofs that live there and that we have
1720  // previously assigned a number to (i.e. the ones on
1721  // the interface)
1722  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
1723  cell->get_dof_indices(local_dof_indices);
1724  for (const auto &local_dof_index : local_dof_indices)
1725  if (local_dof_index != numbers::invalid_dof_index)
1726  renumbering[local_dof_index] = numbers::invalid_dof_index;
1727  }
1728  }
1729 
1730 
1731 
1732  /* -------------- distribute_mg_dofs functionality ------------- */
1733 
1734 
1742  template <int dim, int spacedim>
1743  static types::global_dof_index
1744  distribute_mg_dofs_on_cell(
1745  const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
1746  types::global_dof_index next_free_dof,
1747  const std::integral_constant<int, 1> &)
1748  {
1749  // distribute dofs of vertices
1750  if (cell->get_fe().dofs_per_vertex > 0)
1751  for (const unsigned int v : GeometryInfo<1>::vertex_indices())
1752  {
1753  typename DoFHandler<dim, spacedim>::level_cell_iterator
1754  neighbor = cell->neighbor(v);
1755 
1756  if (neighbor.state() == IteratorState::valid)
1757  {
1758  // has neighbor already been processed?
1759  if (neighbor->user_flag_set() &&
1760  (neighbor->level() == cell->level()))
1761  // copy dofs if the neighbor is on the same level (only
1762  // then are mg dofs the same)
1763  {
1764  if (v == 0)
1765  for (unsigned int d = 0;
1766  d < cell->get_fe().dofs_per_vertex;
1767  ++d)
1768  cell->set_mg_vertex_dof_index(
1769  cell->level(),
1770  0,
1771  d,
1772  neighbor->mg_vertex_dof_index(cell->level(),
1773  1,
1774  d));
1775  else
1776  for (unsigned int d = 0;
1777  d < cell->get_fe().dofs_per_vertex;
1778  ++d)
1779  cell->set_mg_vertex_dof_index(
1780  cell->level(),
1781  1,
1782  d,
1783  neighbor->mg_vertex_dof_index(cell->level(),
1784  0,
1785  d));
1786 
1787  // next neighbor
1788  continue;
1789  }
1790  }
1791 
1792  // otherwise: create dofs newly
1793  for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
1794  ++d)
1795  cell->set_mg_vertex_dof_index(cell->level(),
1796  v,
1797  d,
1798  next_free_dof++);
1799  }
1800 
1801  // dofs of line
1802  if (cell->get_fe().dofs_per_line > 0)
1803  for (unsigned int d = 0; d < cell->get_fe().dofs_per_line; ++d)
1804  cell->set_mg_dof_index(cell->level(), d, next_free_dof++);
1805 
1806  // note that this cell has been processed
1807  cell->set_user_flag();
1808 
1809  return next_free_dof;
1810  }
1811 
1812 
1813 
1814  template <int dim, int spacedim>
1815  static types::global_dof_index
1816  distribute_mg_dofs_on_cell(
1817  const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
1818  types::global_dof_index next_free_dof,
1819  const std::integral_constant<int, 2> &)
1820  {
1821  if (cell->get_fe().dofs_per_vertex > 0)
1822  // number dofs on vertices
1823  for (const unsigned int vertex : GeometryInfo<2>::vertex_indices())
1824  // check whether dofs for this
1825  // vertex have been distributed
1826  // (only check the first dof)
1827  if (cell->mg_vertex_dof_index(cell->level(), vertex, 0) ==
1828  numbers::invalid_dof_index)
1829  for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
1830  ++d)
1831  cell->set_mg_vertex_dof_index(cell->level(),
1832  vertex,
1833  d,
1834  next_free_dof++);
1835 
1836  // for the four sides
1837  if (cell->get_fe().dofs_per_line > 0)
1838  for (const unsigned int side : GeometryInfo<2>::face_indices())
1839  {
1840  typename DoFHandler<dim, spacedim>::line_iterator line =
1841  cell->line(side);
1842 
1843  // distribute dofs if necessary: check whether line dof is
1844  // already numbered (check only first dof)
1845  if (line->mg_dof_index(cell->level(), 0) ==
1846  numbers::invalid_dof_index)
1847  // if not: distribute dofs
1848  for (unsigned int d = 0; d < cell->get_fe().dofs_per_line;
1849  ++d)
1850  line->set_mg_dof_index(cell->level(), d, next_free_dof++);
1851  }
1852 
1853 
1854  // dofs of quad
1855  if (cell->get_fe().dofs_per_quad > 0)
1856  for (unsigned int d = 0; d < cell->get_fe().dofs_per_quad; ++d)
1857  cell->set_mg_dof_index(cell->level(), d, next_free_dof++);
1858 
1859 
1860  // note that this cell has been processed
1861  cell->set_user_flag();
1862 
1863  return next_free_dof;
1864  }
1865 
1866 
1867 
1868  template <int dim, int spacedim>
1869  static types::global_dof_index
1870  distribute_mg_dofs_on_cell(
1871  const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
1872  types::global_dof_index next_free_dof,
1873  const std::integral_constant<int, 3> &)
1874  {
1875  if (cell->get_fe().dofs_per_vertex > 0)
1876  // number dofs on vertices
1877  for (const unsigned int vertex : GeometryInfo<3>::vertex_indices())
1878  // check whether dofs for this vertex have been distributed
1879  // (only check the first dof)
1880  if (cell->mg_vertex_dof_index(cell->level(), vertex, 0) ==
1881  numbers::invalid_dof_index)
1882  for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
1883  ++d)
1884  cell->set_mg_vertex_dof_index(cell->level(),
1885  vertex,
1886  d,
1887  next_free_dof++);
1888 
1889  // for the lines
1890  if (cell->get_fe().dofs_per_line > 0)
1891  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
1892  {
1893  typename DoFHandler<dim, spacedim>::line_iterator line =
1894  cell->line(l);
1895 
1896  // distribute dofs if necessary:
1897  // check whether line dof is already
1898  // numbered (check only first dof)
1899  if (line->mg_dof_index(cell->level(), 0) ==
1900  numbers::invalid_dof_index)
1901  // if not: distribute dofs
1902  for (unsigned int d = 0; d < cell->get_fe().dofs_per_line;
1903  ++d)
1904  line->set_mg_dof_index(cell->level(), d, next_free_dof++);
1905  }
1906 
1907  // for the quads
1908  if (cell->get_fe().dofs_per_quad > 0)
1909  for (unsigned int q = 0; q < GeometryInfo<3>::quads_per_cell; ++q)
1910  {
1911  typename DoFHandler<dim, spacedim>::quad_iterator quad =
1912  cell->quad(q);
1913 
1914  // distribute dofs if necessary:
1915  // check whether line dof is already
1916  // numbered (check only first dof)
1917  if (quad->mg_dof_index(cell->level(), 0) ==
1918  numbers::invalid_dof_index)
1919  // if not: distribute dofs
1920  for (unsigned int d = 0; d < cell->get_fe().dofs_per_quad;
1921  ++d)
1922  quad->set_mg_dof_index(cell->level(), d, next_free_dof++);
1923  }
1924 
1925 
1926  // dofs of cell
1927  if (cell->get_fe().dofs_per_hex > 0)
1928  for (unsigned int d = 0; d < cell->get_fe().dofs_per_hex; ++d)
1929  cell->set_mg_dof_index(cell->level(), d, next_free_dof++);
1930 
1931 
1932  // note that this cell has been processed
1933  cell->set_user_flag();
1934 
1935  return next_free_dof;
1936  }
1937 
1938 
1939 
1940  // same for the hp::DoFHandler
1941  template <int spacedim>
1942  static types::global_dof_index
1943  distribute_mg_dofs_on_cell(
1944  const hp::DoFHandler<1, spacedim> &dof_handler,
1945  const typename hp::DoFHandler<1, spacedim>::active_cell_iterator
1946  & cell,
1947  types::global_dof_index next_free_dof)
1948  {
1949  (void)dof_handler;
1950  (void)cell;
1951  (void)next_free_dof;
1952  return 0;
1953  }
1954 
1955 
1956 
1957  template <int spacedim>
1958  static types::global_dof_index
1959  distribute_mg_dofs_on_cell(
1960  const hp::DoFHandler<2, spacedim> &dof_handler,
1961  const typename hp::DoFHandler<2, spacedim>::active_cell_iterator
1962  & cell,
1963  types::global_dof_index next_free_dof)
1964  {
1965  (void)dof_handler;
1966  (void)cell;
1967  (void)next_free_dof;
1968  return 0;
1969  }
1970 
1971 
1972 
1973  template <int spacedim>
1974  static types::global_dof_index
1975  distribute_mg_dofs_on_cell(
1976  const hp::DoFHandler<3, spacedim> &dof_handler,
1977  const typename hp::DoFHandler<3, spacedim>::active_cell_iterator
1978  & cell,
1979  types::global_dof_index next_free_dof)
1980  {
1981  (void)dof_handler;
1982  (void)cell;
1983  (void)next_free_dof;
1984  return 0;
1985  }
1986 
1987 
1988 
1989  template <class DoFHandlerType>
1990  static types::global_dof_index
1991  distribute_dofs_on_level(const types::subdomain_id level_subdomain_id,
1992  DoFHandlerType & dof_handler,
1993  const unsigned int level)
1994  {
1995  const unsigned int dim = DoFHandlerType::dimension;
1996  const unsigned int spacedim = DoFHandlerType::space_dimension;
1997 
1998  const ::Triangulation<dim, spacedim> &tria =
1999  dof_handler.get_triangulation();
2000  Assert(tria.n_levels() > 0, ExcMessage("Empty triangulation"));
2001  if (level >= tria.n_levels())
2002  return 0; // this is allowed for multigrid
2003 
2004  // Clear user flags because we will need them. But first we save
2005  // them and make sure that we restore them later such that at
2006  // the end of this function the Triangulation will be in the
2007  // same state as it was at the beginning of this function.
2008  std::vector<bool> user_flags;
2009  tria.save_user_flags(user_flags);
2010  const_cast<::Triangulation<dim, spacedim> &>(tria)
2011  .clear_user_flags();
2012 
2013  types::global_dof_index next_free_dof = 0;
2014  typename DoFHandler<dim, spacedim>::level_cell_iterator
2015  cell = dof_handler.begin(level),
2016  endc = dof_handler.end(level);
2017 
2018  for (; cell != endc; ++cell)
2019  if ((level_subdomain_id == numbers::invalid_subdomain_id) ||
2020  (cell->level_subdomain_id() == level_subdomain_id))
2021  next_free_dof =
2022  Implementation::distribute_mg_dofs_on_cell<dim, spacedim>(
2023  cell, next_free_dof, std::integral_constant<int, dim>());
2024 
2025  // finally restore the user flags
2026  const_cast<::Triangulation<dim, spacedim> &>(tria)
2027  .load_user_flags(user_flags);
2028 
2029  return next_free_dof;
2030  }
2031 
2032 
2033 
2034  /* --------------------- renumber_dofs functionality ---------------- */
2035 
2036 
2044  template <int dim, int spacedim>
2045  static void
2046  renumber_vertex_dofs(
2047  const std::vector<types::global_dof_index> &new_numbers,
2048  const IndexSet & indices_we_care_about,
2049  DoFHandler<dim, spacedim> & dof_handler,
2050  const bool check_validity)
2051  {
2052  // we can not use cell iterators in this function since then
2053  // we would renumber the dofs on the interface of two cells
2054  // more than once. Anyway, this way it's not only more
2055  // correct but also faster; note, however, that dof numbers
2056  // may be invalid_dof_index, namely when the appropriate
2057  // vertex/line/etc is unused
2058  for (std::vector<types::global_dof_index>::iterator i =
2059  dof_handler.vertex_dofs.begin();
2060  i != dof_handler.vertex_dofs.end();
2061  ++i)
2062  if (*i != numbers::invalid_dof_index)
2063  *i = (indices_we_care_about.size() == 0) ?
2064  (new_numbers[*i]) :
2065  (new_numbers[indices_we_care_about.index_within_set(*i)]);
2066  else if (check_validity)
2067  // if index is invalid_dof_index: check if this one
2068  // really is unused
2069  Assert(dof_handler.get_triangulation().vertex_used(
2070  (i - dof_handler.vertex_dofs.begin()) /
2071  dof_handler.get_fe().dofs_per_vertex) == false,
2072  ExcInternalError());
2073  }
2074 
2075 
2076 
2084  template <int dim, int spacedim>
2085  static void
2086  renumber_cell_dofs(
2087  const std::vector<types::global_dof_index> &new_numbers,
2088  const IndexSet & indices_we_care_about,
2089  DoFHandler<dim, spacedim> & dof_handler)
2090  {
2091  for (unsigned int level = 0; level < dof_handler.levels.size();
2092  ++level)
2093  for (std::vector<types::global_dof_index>::iterator i =
2094  dof_handler.levels[level]->dof_object.dofs.begin();
2095  i != dof_handler.levels[level]->dof_object.dofs.end();
2096  ++i)
2097  if (*i != numbers::invalid_dof_index)
2098  *i =
2099  ((indices_we_care_about.size() == 0) ?
2100  new_numbers[*i] :
2101  new_numbers[indices_we_care_about.index_within_set(*i)]);
2102  }
2103 
2104 
2105 
2113  template <int spacedim>
2114  static void
2115  renumber_face_dofs(
2116  const std::vector<types::global_dof_index> & /*new_numbers*/,
2117  const IndexSet & /*indices_we_care_about*/,
2118  DoFHandler<1, spacedim> & /*dof_handler*/)
2119  {
2120  // nothing to do in 1d since there are no separate faces
2121  }
2122 
2123 
2124 
2125  template <int spacedim>
2126  static void
2127  renumber_face_dofs(
2128  const std::vector<types::global_dof_index> &new_numbers,
2129  const IndexSet & indices_we_care_about,
2130  DoFHandler<2, spacedim> & dof_handler)
2131  {
2132  // treat dofs on lines
2133  for (std::vector<types::global_dof_index>::iterator i =
2134  dof_handler.faces->lines.dofs.begin();
2135  i != dof_handler.faces->lines.dofs.end();
2136  ++i)
2137  if (*i != numbers::invalid_dof_index)
2138  *i = ((indices_we_care_about.size() == 0) ?
2139  new_numbers[*i] :
2140  new_numbers[indices_we_care_about.index_within_set(*i)]);
2141  }
2142 
2143 
2144 
2145  template <int spacedim>
2146  static void
2147  renumber_face_dofs(
2148  const std::vector<types::global_dof_index> &new_numbers,
2149  const IndexSet & indices_we_care_about,
2150  DoFHandler<3, spacedim> & dof_handler)
2151  {
2152  // treat dofs on lines
2153  for (std::vector<types::global_dof_index>::iterator i =
2154  dof_handler.faces->lines.dofs.begin();
2155  i != dof_handler.faces->lines.dofs.end();
2156  ++i)
2157  if (*i != numbers::invalid_dof_index)
2158  *i = ((indices_we_care_about.size() == 0) ?
2159  new_numbers[*i] :
2160  new_numbers[indices_we_care_about.index_within_set(*i)]);
2161 
2162  // treat dofs on quads
2163  for (std::vector<types::global_dof_index>::iterator i =
2164  dof_handler.faces->quads.dofs.begin();
2165  i != dof_handler.faces->quads.dofs.end();
2166  ++i)
2167  if (*i != numbers::invalid_dof_index)
2168  *i = ((indices_we_care_about.size() == 0) ?
2169  new_numbers[*i] :
2170  new_numbers[indices_we_care_about.index_within_set(*i)]);
2171  }
2172 
2173 
2174 
2175  template <int dim, int spacedim>
2176  static void
2177  renumber_vertex_dofs(
2178  const std::vector<types::global_dof_index> &new_numbers,
2179  const IndexSet & indices_we_care_about,
2180  hp::DoFHandler<dim, spacedim> & dof_handler,
2181  const bool check_validity)
2182  {
2183  for (unsigned int vertex_index = 0;
2184  vertex_index < dof_handler.get_triangulation().n_vertices();
2185  ++vertex_index)
2186  {
2187  const unsigned int n_active_fe_indices =
2188  ::internal::DoFAccessorImplementation::Implementation::
2189  n_active_vertex_fe_indices(dof_handler, vertex_index);
2190 
2191  // if this vertex is unused, then we really ought not to have
2192  // allocated any space for it, i.e., n_active_fe_indices should be
2193  // zero, and there is no space to actually store dof indices for
2194  // this vertex
2195  if (dof_handler.get_triangulation().vertex_used(vertex_index) ==
2196  false)
2197  Assert(n_active_fe_indices == 0, ExcInternalError());
2198 
2199  // otherwise the vertex is used; it may still not hold any dof
2200  // indices if it is located on an artificial cell and not adjacent
2201  // to a ghost cell, but in that case there is simply nothing for
2202  // us to do
2203  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2204  {
2205  const unsigned int fe_index =
2206  ::internal::DoFAccessorImplementation::
2207  Implementation::nth_active_vertex_fe_index(dof_handler,
2208  vertex_index,
2209  f);
2210 
2211  for (unsigned int d = 0;
2212  d < dof_handler.get_fe(fe_index).dofs_per_vertex;
2213  ++d)
2214  {
2215  const types::global_dof_index old_dof_index =
2216  ::internal::DoFAccessorImplementation::
2217  Implementation::get_vertex_dof_index(dof_handler,
2218  vertex_index,
2219  fe_index,
2220  d);
2221 
2222  // if check_validity was set, then we are to verify that
2223  // the previous indices were all valid. this really should
2224  // be the case: we allocated space for these vertex dofs,
2225  // i.e., at least one adjacent cell has a valid
2226  // active_fe_index, so there are DoFs that really live
2227  // on this vertex. if check_validity is set, then we
2228  // must make sure that they have been set to something
2229  // useful
2230  if (check_validity)
2231  Assert(old_dof_index != numbers::invalid_dof_index,
2232  ExcInternalError());
2233 
2234  if (old_dof_index != numbers::invalid_dof_index)
2235  {
2236  // In the following blocks, we first check whether
2237  // we were given an IndexSet of DoFs to touch. If not
2238  // (the first 'if' case here), then we are in the
2239  // sequential case and are allowed to touch all DoFs.
2240  //
2241  // If yes (the 'else' case), then we need to
2242  // distinguish whether the DoF whose number we want to
2243  // touch is in fact locally owned (i.e., is in the
2244  // index set) and then we can actually assign it a new
2245  // number; otherwise, we have encountered a
2246  // non-locally owned DoF for which we don't know the
2247  // new number yet and so set it to an invalid index.
2248  // This will later be fixed up after the first ghost
2249  // exchange phase when we unify hp DoFs on neighboring
2250  // cells.
2251  if (indices_we_care_about.size() == 0)
2252  ::internal::DoFAccessorImplementation::
2253  Implementation::set_vertex_dof_index(
2254  dof_handler,
2255  vertex_index,
2256  fe_index,
2257  d,
2258  new_numbers[old_dof_index]);
2259  else
2260  {
2261  if (indices_we_care_about.is_element(
2262  old_dof_index))
2263  ::internal::DoFAccessorImplementation::
2264  Implementation::set_vertex_dof_index(
2265  dof_handler,
2266  vertex_index,
2267  fe_index,
2268  d,
2269  new_numbers[indices_we_care_about
2270  .index_within_set(
2271  old_dof_index)]);
2272  else
2273  ::internal::DoFAccessorImplementation::
2274  Implementation::set_vertex_dof_index(
2275  dof_handler,
2276  vertex_index,
2277  fe_index,
2278  d,
2279  numbers::invalid_dof_index);
2280  }
2281  }
2282  }
2283  }
2284  }
2285  }
2286 
2287 
2288 
2289  template <int dim, int spacedim>
2290  static void
2291  renumber_cell_dofs(
2292  const std::vector<types::global_dof_index> &new_numbers,
2293  const IndexSet & indices_we_care_about,
2294  hp::DoFHandler<dim, spacedim> & dof_handler)
2295  {
2296  for (const auto &cell : dof_handler.active_cell_iterators())
2297  if (!cell->is_artificial())
2298  {
2299  const unsigned int fe_index = cell->active_fe_index();
2300 
2301  for (unsigned int d = 0;
2302  d < dof_handler.get_fe(fe_index)
2303  .template n_dofs_per_object<dim>();
2304  ++d)
2305  {
2306  const types::global_dof_index old_dof_index =
2307  cell->dof_index(d, fe_index);
2308  if (old_dof_index != numbers::invalid_dof_index)
2309  {
2310  // In the following blocks, we first check whether
2311  // we were given an IndexSet of DoFs to touch. If not
2312  // (the first 'if' case here), then we are in the
2313  // sequential case and are allowed to touch all DoFs.
2314  //
2315  // If yes (the 'else' case), then we need to distinguish
2316  // whether the DoF whose number we want to touch is in
2317  // fact locally owned (i.e., is in the index set) and
2318  // then we can actually assign it a new number;
2319  // otherwise, we have encountered a non-locally owned
2320  // DoF for which we don't know the new number yet and so
2321  // set it to an invalid index. This will later be fixed
2322  // up after the first ghost exchange phase when we unify
2323  // hp DoFs on neighboring cells.
2324  if (indices_we_care_about.size() == 0)
2325  cell->set_dof_index(d,
2326  new_numbers[old_dof_index],
2327  fe_index);
2328  else
2329  {
2330  if (indices_we_care_about.is_element(old_dof_index))
2331  cell->set_dof_index(
2332  d,
2333  new_numbers[indices_we_care_about
2334  .index_within_set(old_dof_index)],
2335  fe_index);
2336  else
2337  cell->set_dof_index(d,
2338  numbers::invalid_dof_index,
2339  fe_index);
2340  }
2341  }
2342  }
2343  }
2344  }
2345 
2346 
2347 
2348  template <int spacedim>
2349  static void
2350  renumber_face_dofs(
2351  const std::vector<types::global_dof_index> & /*new_numbers*/,
2352  const IndexSet & /*indices_we_care_about*/,
2353  hp::DoFHandler<1, spacedim> & /*dof_handler*/)
2354  {
2355  // nothing to do in 1d since there are no separate faces -- we've
2356  // already taken care of this when dealing with the vertices
2357  }
2358 
2359 
2360 
2361  template <int spacedim>
2362  static void
2363  renumber_face_dofs(
2364  const std::vector<types::global_dof_index> &new_numbers,
2365  const IndexSet & indices_we_care_about,
2366  hp::DoFHandler<2, spacedim> & dof_handler)
2367  {
2368  const unsigned int dim = 2;
2369 
2370  // deal with DoFs on lines
2371  {
2372  // save user flags on lines so we can use them to mark lines
2373  // we've already treated
2374  std::vector<bool> saved_line_user_flags;
2375  const_cast<::Triangulation<dim, spacedim> &>(
2376  dof_handler.get_triangulation())
2377  .save_user_flags_line(saved_line_user_flags);
2378  const_cast<::Triangulation<dim, spacedim> &>(
2379  dof_handler.get_triangulation())
2380  .clear_user_flags_line();
2381 
2382  for (const auto &cell : dof_handler.active_cell_iterators())
2383  if (!cell->is_artificial())
2384  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
2385  ++l)
2386  if (cell->line(l)->user_flag_set() == false)
2387  {
2388  const typename hp::DoFHandler<dim,
2389  spacedim>::line_iterator
2390  line = cell->line(l);
2391  line->set_user_flag();
2392 
2393  const unsigned int n_active_fe_indices =
2394  line->n_active_fe_indices();
2395 
2396  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2397  {
2398  const unsigned int fe_index =
2399  line->nth_active_fe_index(f);
2400 
2401  for (unsigned int d = 0;
2402  d < dof_handler.get_fe(fe_index).dofs_per_line;
2403  ++d)
2404  {
2405  const types::global_dof_index old_dof_index =
2406  line->dof_index(d, fe_index);
2407  if (old_dof_index != numbers::invalid_dof_index)
2408  {
2409  // In the following blocks, we first check
2410  // whether we were given an IndexSet of DoFs
2411  // to touch. If not (the first 'if' case
2412  // here), then we are in the sequential case
2413  // and are allowed to touch all DoFs.
2414  //
2415  // If yes (the 'else' case), then we need to
2416  // distinguish whether the DoF whose number we
2417  // want to touch is in fact locally owned
2418  // (i.e., is in the index set) and then we can
2419  // actually assign it a new number; otherwise,
2420  // we have encountered a non-locally owned DoF
2421  // for which we don't know the new number yet
2422  // and so set it to an invalid index. This
2423  // will later be fixed up after the first
2424  // ghost exchange phase when we unify hp DoFs
2425  // on neighboring cells.
2426  if (indices_we_care_about.size() == 0)
2427  line->set_dof_index(
2428  d, new_numbers[old_dof_index], fe_index);
2429  else
2430  {
2431  if (indices_we_care_about.is_element(
2432  old_dof_index))
2433  line->set_dof_index(
2434  d,
2435  new_numbers[indices_we_care_about
2436  .index_within_set(
2437  old_dof_index)],
2438  fe_index);
2439  else
2440  line->set_dof_index(
2441  d,
2442  numbers::invalid_dof_index,
2443  fe_index);
2444  }
2445  }
2446  }
2447  }
2448  }
2449 
2450  // at the end, restore the user
2451  // flags for the lines
2452  const_cast<::Triangulation<dim, spacedim> &>(
2453  dof_handler.get_triangulation())
2454  .load_user_flags_line(saved_line_user_flags);
2455  }
2456  }
2457 
2458 
2459 
2460  template <int spacedim>
2461  static void
2462  renumber_face_dofs(
2463  const std::vector<types::global_dof_index> &new_numbers,
2464  const IndexSet & indices_we_care_about,
2465  hp::DoFHandler<3, spacedim> & dof_handler)
2466  {
2467  const unsigned int dim = 3;
2468 
2469  // deal with DoFs on lines
2470  {
2471  // save user flags on lines so we can use them to mark lines
2472  // we've already treated
2473  std::vector<bool> saved_line_user_flags;
2474  const_cast<::Triangulation<dim, spacedim> &>(
2475  dof_handler.get_triangulation())
2476  .save_user_flags_line(saved_line_user_flags);
2477  const_cast<::Triangulation<dim, spacedim> &>(
2478  dof_handler.get_triangulation())
2479  .clear_user_flags_line();
2480 
2481  for (const auto &cell : dof_handler.active_cell_iterators())
2482  if (!cell->is_artificial())
2483  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
2484  ++l)
2485  if (cell->line(l)->user_flag_set() == false)
2486  {
2487  const typename hp::DoFHandler<dim,
2488  spacedim>::line_iterator
2489  line = cell->line(l);
2490  line->set_user_flag();
2491 
2492  const unsigned int n_active_fe_indices =
2493  line->n_active_fe_indices();
2494 
2495  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2496  {
2497  const unsigned int fe_index =
2498  line->nth_active_fe_index(f);
2499 
2500  for (unsigned int d = 0;
2501  d < dof_handler.get_fe(fe_index).dofs_per_line;
2502  ++d)
2503  {
2504  const types::global_dof_index old_dof_index =
2505  line->dof_index(d, fe_index);
2506  if (old_dof_index != numbers::invalid_dof_index)
2507  {
2508  // In the following blocks, we first check
2509  // whether we were given an IndexSet of DoFs
2510  // to touch. If not (the first 'if' case
2511  // here), then we are in the sequential case
2512  // and are allowed to touch all DoFs.
2513  //
2514  // If yes (the 'else' case), then we need to
2515  // distinguish whether the DoF whose number we
2516  // want to touch is in fact locally owned
2517  // (i.e., is in the index set) and then we can
2518  // actually assign it a new number; otherwise,
2519  // we have encountered a non-locally owned DoF
2520  // for which we don't know the new number yet
2521  // and so set it to an invalid index. This
2522  // will later be fixed up after the first
2523  // ghost exchange phase when we unify hp DoFs
2524  // on neighboring cells.
2525  if (indices_we_care_about.size() == 0)
2526  line->set_dof_index(
2527  d, new_numbers[old_dof_index], fe_index);
2528  else if (indices_we_care_about.is_element(
2529  old_dof_index))
2530  line->set_dof_index(
2531  d,
2532  new_numbers[indices_we_care_about
2533  .index_within_set(
2534  old_dof_index)],
2535  fe_index);
2536  else
2537  line->set_dof_index(
2538  d, numbers::invalid_dof_index, fe_index);
2539  }
2540  }
2541  }
2542  }
2543 
2544  // at the end, restore the user
2545  // flags for the lines
2546  const_cast<::Triangulation<dim, spacedim> &>(
2547  dof_handler.get_triangulation())
2548  .load_user_flags_line(saved_line_user_flags);
2549  }
2550 
2551  // then deal with dofs on quads
2552  {
2553  std::vector<bool> saved_quad_user_flags;
2554  const_cast<::Triangulation<dim, spacedim> &>(
2555  dof_handler.get_triangulation())
2556  .save_user_flags_quad(saved_quad_user_flags);
2557  const_cast<::Triangulation<dim, spacedim> &>(
2558  dof_handler.get_triangulation())
2559  .clear_user_flags_quad();
2560 
2561  for (const auto &cell : dof_handler.active_cell_iterators())
2562  if (!cell->is_artificial())
2563  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell;
2564  ++q)
2565  if (cell->quad(q)->user_flag_set() == false)
2566  {
2567  const typename hp::DoFHandler<dim,
2568  spacedim>::quad_iterator
2569  quad = cell->quad(q);
2570  quad->set_user_flag();
2571 
2572  const unsigned int n_active_fe_indices =
2573  quad->n_active_fe_indices();
2574 
2575  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2576  {
2577  const unsigned int fe_index =
2578  quad->nth_active_fe_index(f);
2579 
2580  for (unsigned int d = 0;
2581  d < dof_handler.get_fe(fe_index).dofs_per_quad;
2582  ++d)
2583  {
2584  const types::global_dof_index old_dof_index =
2585  quad->dof_index(d, fe_index);
2586  if (old_dof_index != numbers::invalid_dof_index)
2587  {
2588  // In the following blocks, we first check
2589  // whether we were given an IndexSet of DoFs
2590  // to touch. If not (the first 'if' case
2591  // here), then we are in the sequential case
2592  // and are allowed to touch all DoFs.
2593  //
2594  // If yes (the 'else' case), then we need to
2595  // distinguish whether the DoF whose number we
2596  // want to touch is in fact locally owned
2597  // (i.e., is in the index set) and then we can
2598  // actually assign it a new number; otherwise,
2599  // we have encountered a non-locally owned DoF
2600  // for which we don't know the new number yet
2601  // and so set it to an invalid index. This
2602  // will later be fixed up after the first
2603  // ghost exchange phase when we unify hp DoFs
2604  // on neighboring cells.
2605  if (indices_we_care_about.size() == 0)
2606  quad->set_dof_index(
2607  d, new_numbers[old_dof_index], fe_index);
2608  else
2609  {
2610  if (indices_we_care_about.is_element(
2611  old_dof_index))
2612  quad->set_dof_index(
2613  d,
2614  new_numbers[indices_we_care_about
2615  .index_within_set(
2616  old_dof_index)],
2617  fe_index);
2618  else
2619  quad->set_dof_index(
2620  d,
2621  numbers::invalid_dof_index,
2622  fe_index);
2623  }
2624  }
2625  }
2626  }
2627  }
2628 
2629  // at the end, restore the user flags for the quads
2630  const_cast<::Triangulation<dim, spacedim> &>(
2631  dof_handler.get_triangulation())
2632  .load_user_flags_quad(saved_quad_user_flags);
2633  }
2634  }
2635 
2636 
2637 
2649  template <class DoFHandlerType>
2650  static void
2651  renumber_dofs(const std::vector<types::global_dof_index> &new_numbers,
2652  const IndexSet &indices_we_care_about,
2653  DoFHandlerType &dof_handler,
2654  const bool check_validity)
2655  {
2656  if (DoFHandlerType::dimension == 1)
2657  Assert(indices_we_care_about == IndexSet(0), ExcNotImplemented());
2658 
2659  // renumber DoF indices on vertices, cells, and faces. this
2660  // can be done in parallel because the respective functions
2661  // work on separate data structures
2662  Threads::TaskGroup<> tasks;
2663  tasks += Threads::new_task([&]() {
2664  renumber_vertex_dofs(new_numbers,
2665  indices_we_care_about,
2666  dof_handler,
2667  check_validity);
2668  });
2669  tasks += Threads::new_task([&]() {
2670  renumber_face_dofs(new_numbers, indices_we_care_about, dof_handler);
2671  });
2672  tasks += Threads::new_task([&]() {
2673  renumber_cell_dofs(new_numbers, indices_we_care_about, dof_handler);
2674  });
2675  tasks.join_all();
2676 
2677  // update the cache used for cell dof indices
2678  update_all_active_cell_dof_indices_caches(dof_handler);
2679  }
2680 
2681 
2682 
2683  /* --------------------- renumber_mg_dofs functionality ----------------
2684  */
2685 
2693  template <int dim, int spacedim>
2694  static void
2695  renumber_vertex_mg_dofs(
2696  const std::vector<::types::global_dof_index> &new_numbers,
2697  const IndexSet & indices_we_care_about,
2698  DoFHandler<dim, spacedim> &dof_handler,
2699  const unsigned int level,
2700  const bool check_validity)
2701  {
2702  (void)check_validity;
2703  Assert(level < dof_handler.get_triangulation().n_levels(),
2704  ExcInternalError());
2705 
2706  for (typename std::vector<
2707  typename DoFHandler<dim, spacedim>::MGVertexDoFs>::iterator i =
2708  dof_handler.mg_vertex_dofs.begin();
2709  i != dof_handler.mg_vertex_dofs.end();
2710  ++i)
2711  // if the present vertex lives on the current level
2712  if ((i->get_coarsest_level() <= level) &&
2713  (i->get_finest_level() >= level))
2714  for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_vertex;
2715  ++d)
2716  {
2717  const ::types::global_dof_index idx =
2718  i->get_index(level,
2719  d,
2720  dof_handler.get_fe().dofs_per_vertex);
2721 
2722  if (idx != numbers::invalid_dof_index)
2723  {
2724  Assert(check_validity == false ||
2725  (indices_we_care_about.size() > 0 ?
2726  indices_we_care_about.is_element(idx) :
2727  (idx < new_numbers.size())),
2728  ExcInternalError());
2729  i->set_index(level,
2730  d,
2731  dof_handler.get_fe().dofs_per_vertex,
2732  (indices_we_care_about.size() == 0) ?
2733  (new_numbers[idx]) :
2734  (new_numbers[indices_we_care_about
2735  .index_within_set(idx)]));
2736  }
2737  }
2738  }
2739 
2740 
2741 
2749  template <int dim, int spacedim>
2750  static void
2751  renumber_cell_mg_dofs(
2752  const std::vector<::types::global_dof_index> &new_numbers,
2753  const IndexSet & indices_we_care_about,
2754  DoFHandler<dim, spacedim> &dof_handler,
2755  const unsigned int level)
2756  {
2757  for (std::vector<types::global_dof_index>::iterator i =
2758  dof_handler.mg_levels[level]->dof_object.dofs.begin();
2759  i != dof_handler.mg_levels[level]->dof_object.dofs.end();
2760  ++i)
2761  {
2762  if (*i != numbers::invalid_dof_index)
2763  {
2764  Assert((indices_we_care_about.size() > 0 ?
2765  indices_we_care_about.is_element(*i) :
2766  (*i < new_numbers.size())),
2767  ExcInternalError());
2768  *i =
2769  (indices_we_care_about.size() == 0) ?
2770  (new_numbers[*i]) :
2771  (new_numbers[indices_we_care_about.index_within_set(*i)]);
2772  }
2773  }
2774  }
2775 
2776 
2777 
2785  template <int spacedim>
2786  static void
2787  renumber_face_mg_dofs(
2788  const std::vector<types::global_dof_index> & /*new_numbers*/,
2789  const IndexSet & /*indices_we_care_about*/,
2790  DoFHandler<1, spacedim> & /*dof_handler*/,
2791  const unsigned int /*level*/,
2792  const bool /*check_validity*/)
2793  {
2794  // nothing to do in 1d because there are no separate faces
2795  }
2796 
2797 
2798 
2799  template <int spacedim>
2800  static void
2801  renumber_face_mg_dofs(
2802  const std::vector<::types::global_dof_index> &new_numbers,
2803  const IndexSet & indices_we_care_about,
2804  DoFHandler<2, spacedim> &dof_handler,
2805  const unsigned int level,
2806  const bool check_validity)
2807  {
2808  if (dof_handler.get_fe().dofs_per_line > 0)
2809  {
2810  // save user flags as they will be modified
2811  std::vector<bool> user_flags;
2812  dof_handler.get_triangulation().save_user_flags(user_flags);
2813  const_cast<::Triangulation<2, spacedim> &>(
2814  dof_handler.get_triangulation())
2815  .clear_user_flags();
2816 
2817  // flag all lines adjacent to cells of the current
2818  // level, as those lines logically belong to the same
2819  // level as the cell, at least for for isotropic
2820  // refinement
2821  typename DoFHandler<2, spacedim>::level_cell_iterator cell,
2822  endc = dof_handler.end(level);
2823  for (cell = dof_handler.begin(level); cell != endc; ++cell)
2824  if (cell->level_subdomain_id() !=
2825  numbers::artificial_subdomain_id)
2826  for (const unsigned int line :
2827  GeometryInfo<2>::face_indices())
2828  cell->face(line)->set_user_flag();
2829 
2830  for (typename DoFHandler<2, spacedim>::cell_iterator cell =
2831  dof_handler.begin();
2832  cell != dof_handler.end();
2833  ++cell)
2834  for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell;
2835  ++l)
2836  if (cell->line(l)->user_flag_set())
2837  {
2838  for (unsigned int d = 0;
2839  d < dof_handler.get_fe().dofs_per_line;
2840  ++d)
2841  {
2842  const ::types::global_dof_index idx =
2843  cell->line(l)->mg_dof_index(level, d);
2844  if (check_validity)
2845  Assert(idx != numbers::invalid_dof_index,
2846  ExcInternalError());
2847 
2848  if (idx != numbers::invalid_dof_index)
2849  cell->line(l)->set_mg_dof_index(
2850  level,
2851  d,
2852  ((indices_we_care_about.size() == 0) ?
2853  new_numbers[idx] :
2854  new_numbers[indices_we_care_about
2855  .index_within_set(idx)]));
2856  }
2857  cell->line(l)->clear_user_flag();
2858  }
2859  // finally, restore user flags
2860  const_cast<::Triangulation<2, spacedim> &>(
2861  dof_handler.get_triangulation())
2862  .load_user_flags(user_flags);
2863  }
2864  }
2865 
2866 
2867 
2868  template <int spacedim>
2869  static void
2870  renumber_face_mg_dofs(
2871  const std::vector<::types::global_dof_index> &new_numbers,
2872  const IndexSet & indices_we_care_about,
2873  DoFHandler<3, spacedim> &dof_handler,
2874  const unsigned int level,
2875  const bool check_validity)
2876  {
2877  if (dof_handler.get_fe().dofs_per_line > 0 ||
2878  dof_handler.get_fe().dofs_per_quad > 0)
2879  {
2880  // save user flags as they will be modified
2881  std::vector<bool> user_flags;
2882  dof_handler.get_triangulation().save_user_flags(user_flags);
2883  const_cast<::Triangulation<3, spacedim> &>(
2884  dof_handler.get_triangulation())
2885  .clear_user_flags();
2886 
2887  // flag all lines adjacent to cells of the current
2888  // level, as those lines logically belong to the same
2889  // level as the cell, at least for isotropic refinement
2890  typename DoFHandler<3, spacedim>::level_cell_iterator cell,
2891  endc = dof_handler.end(level);
2892  for (cell = dof_handler.begin(level); cell != endc; ++cell)
2893  if (cell->level_subdomain_id() !=
2894  numbers::artificial_subdomain_id)
2895  for (unsigned int line = 0;
2896  line < GeometryInfo<3>::lines_per_cell;
2897  ++line)
2898  cell->line(line)->set_user_flag();
2899 
2900  for (typename DoFHandler<3, spacedim>::cell_iterator cell =
2901  dof_handler.begin();
2902  cell != dof_handler.end();
2903  ++cell)
2904  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell;
2905  ++l)
2906  if (cell->line(l)->user_flag_set())
2907  {
2908  for (unsigned int d = 0;
2909  d < dof_handler.get_fe().dofs_per_line;
2910  ++d)
2911  {
2912  const ::types::global_dof_index idx =
2913  cell->line(l)->mg_dof_index(level, d);
2914  if (check_validity)
2915  Assert(idx != numbers::invalid_dof_index,
2916  ExcInternalError());
2917 
2918  if (idx != numbers::invalid_dof_index)
2919  cell->line(l)->set_mg_dof_index(
2920  level,
2921  d,
2922  ((indices_we_care_about.size() == 0) ?
2923  new_numbers[idx] :
2924  new_numbers[indices_we_care_about
2925  .index_within_set(idx)]));
2926  }
2927  cell->line(l)->clear_user_flag();
2928  }
2929 
2930  // flag all quads adjacent to cells of the current level, as
2931  // those quads logically belong to the same level as the cell,
2932  // at least for isotropic refinement
2933  for (cell = dof_handler.begin(level); cell != endc; ++cell)
2934  if (cell->level_subdomain_id() !=
2935  numbers::artificial_subdomain_id)
2936  for (unsigned int quad = 0;
2937  quad < GeometryInfo<3>::quads_per_cell;
2938  ++quad)
2939  cell->quad(quad)->set_user_flag();
2940 
2941  for (typename DoFHandler<3, spacedim>::cell_iterator cell =
2942  dof_handler.begin();
2943  cell != dof_handler.end();
2944  ++cell)
2945  for (unsigned int l = 0; l < GeometryInfo<3>::quads_per_cell;
2946  ++l)
2947  if (cell->quad(l)->user_flag_set())
2948  {
2949  for (unsigned int d = 0;
2950  d < dof_handler.get_fe().dofs_per_quad;
2951  ++d)
2952  {
2953  const ::types::global_dof_index idx =
2954  cell->quad(l)->mg_dof_index(level, d);
2955  if (check_validity)
2956  Assert(idx != numbers::invalid_dof_index,
2957  ExcInternalError());
2958 
2959  if (idx != numbers::invalid_dof_index)
2960  cell->quad(l)->set_mg_dof_index(
2961  level,
2962  d,
2963  ((indices_we_care_about.size() == 0) ?
2964  new_numbers[idx] :
2965  new_numbers[indices_we_care_about
2966  .index_within_set(idx)]));
2967  }
2968  cell->quad(l)->clear_user_flag();
2969  }
2970 
2971  // finally, restore user flags
2972  const_cast<::Triangulation<3, spacedim> &>(
2973  dof_handler.get_triangulation())
2974  .load_user_flags(user_flags);
2975  }
2976  }
2977 
2978 
2979 
2980  template <int dim, int spacedim>
2981  static void
2982  renumber_mg_dofs(
2983  const std::vector<::types::global_dof_index> &new_numbers,
2984  const IndexSet & indices_we_care_about,
2985  DoFHandler<dim, spacedim> &dof_handler,
2986  const unsigned int level,
2987  const bool check_validity)
2988  {
2989  Assert(level < dof_handler.get_triangulation().n_global_levels(),
2990  ExcInternalError());
2991 
2992  // renumber DoF indices on vertices, cells, and faces. this
2993  // can be done in parallel because the respective functions
2994  // work on separate data structures
2995  Threads::TaskGroup<> tasks;
2996  tasks += Threads::new_task([&]() {
2997  renumber_vertex_mg_dofs(new_numbers,
2998  indices_we_care_about,
2999  dof_handler,
3000  level,
3001  check_validity);
3002  });
3003  tasks += Threads::new_task([&]() {
3004  renumber_face_mg_dofs(new_numbers,
3005  indices_we_care_about,
3006  dof_handler,
3007  level,
3008  check_validity);
3009  });
3010  tasks += Threads::new_task([&]() {
3011  renumber_cell_mg_dofs(new_numbers,
3012  indices_we_care_about,
3013  dof_handler,
3014  level);
3015  });
3016  tasks.join_all();
3017  }
3018 
3019 
3020 
3021  template <int dim, int spacedim>
3022  static void
3023  renumber_mg_dofs(
3024  const std::vector<::types::global_dof_index> & /*new_numbers*/,
3025  const IndexSet & /*indices_we_care_about*/,
3026  hp::DoFHandler<dim, spacedim> & /*dof_handler*/,
3027  const unsigned int /*level*/,
3028  const bool /*check_validity*/)
3029  {
3030  Assert(false, ExcNotImplemented());
3031  }
3032  };
3033 
3034 
3035 
3036  /* --------------------- class Sequential ---------------- */
3037 
3038 
3039 
3040  template <class DoFHandlerType>
3041  Sequential<DoFHandlerType>::Sequential(DoFHandlerType &dof_handler)
3042  : dof_handler(&dof_handler)
3043  {}
3044 
3045 
3046 
3047  template <class DoFHandlerType>
3048  NumberCache
3049  Sequential<DoFHandlerType>::distribute_dofs() const
3050  {
3051  const types::global_dof_index n_initial_dofs =
3052  Implementation::distribute_dofs(numbers::invalid_subdomain_id,
3053  *dof_handler);
3054 
3055  const types::global_dof_index n_dofs =
3056  Implementation::unify_dof_indices(*dof_handler,
3057  n_initial_dofs,
3058  /*check_validity=*/true);
3059 
3060  // return a sequential, complete index set
3061  return NumberCache(n_dofs);
3062  }
3063 
3064 
3065 
3066  template <class DoFHandlerType>
3067  std::vector<NumberCache>
3068  Sequential<DoFHandlerType>::distribute_mg_dofs() const
3069  {
3070  std::vector<bool> user_flags;
3071  dof_handler->get_triangulation().save_user_flags(user_flags);
3072 
3073  const_cast<::Triangulation<DoFHandlerType::dimension,
3074  DoFHandlerType::space_dimension> &>(
3075  dof_handler->get_triangulation())
3076  .clear_user_flags();
3077 
3078  std::vector<NumberCache> number_caches;
3079  number_caches.reserve(dof_handler->get_triangulation().n_levels());
3080  for (unsigned int level = 0;
3081  level < dof_handler->get_triangulation().n_levels();
3082  ++level)
3083  {
3084  // first distribute dofs on this level
3085  const types::global_dof_index n_level_dofs =
3086  Implementation::distribute_dofs_on_level(
3087  numbers::invalid_subdomain_id, *dof_handler, level);
3088 
3089  // then add a complete, sequential index set
3090  number_caches.emplace_back(n_level_dofs);
3091  }
3092 
3093  const_cast<::Triangulation<DoFHandlerType::dimension,
3094  DoFHandlerType::space_dimension> &>(
3095  dof_handler->get_triangulation())
3096  .load_user_flags(user_flags);
3097 
3098  return number_caches;
3099  }
3100 
3101 
3102 
3103  template <class DoFHandlerType>
3104  NumberCache
3105  Sequential<DoFHandlerType>::renumber_dofs(
3106  const std::vector<types::global_dof_index> &new_numbers) const
3107  {
3108  Implementation::renumber_dofs(new_numbers,
3109  IndexSet(0),
3110  *dof_handler,
3111  /*check_validity=*/true);
3112 
3113  // return a sequential, complete index set. take into account that the
3114  // number of DoF indices may in fact be smaller than there were before
3115  // if some previously separately numbered dofs have been identified.
3116  // this is, for example, what the hp::DoFHandler does: it first
3117  // enumerates all DoFs on cells independently, and then unifies
3118  // some located at vertices or faces; this leaves us with fewer
3119  // DoFs than there were before, so use the largest index as
3120  // the one to determine the size of the index space
3121  return NumberCache(
3122  *std::max_element(new_numbers.begin(), new_numbers.end()) + 1);
3123  }
3124 
3125 
3126 
3127  template <class DoFHandlerType>
3128  NumberCache
3129  Sequential<DoFHandlerType>::renumber_mg_dofs(
3130  const unsigned int level,
3131  const std::vector<types::global_dof_index> &new_numbers) const
3132  {
3133  Implementation::renumber_mg_dofs(
3134  new_numbers, IndexSet(0), *dof_handler, level, true);
3135 
3136  // return a sequential, complete index set
3137  return NumberCache(new_numbers.size());
3138  }
3139 
3140 
3141  /* --------------------- class ParallelShared ---------------- */
3142 
3143 
3144  template <class DoFHandlerType>
3145  ParallelShared<DoFHandlerType>::ParallelShared(
3146  DoFHandlerType &dof_handler)
3147  : dof_handler(&dof_handler)
3148  {}
3149 
3150 
3151 
3152  namespace
3153  {
3162  template <class DoFHandlerType>
3163  std::vector<types::subdomain_id>
3164  get_dof_subdomain_association(const DoFHandlerType & dof_handler,
3165  const types::global_dof_index n_dofs,
3166  const unsigned int n_procs)
3167  {
3168  (void)n_procs;
3169  std::vector<types::subdomain_id> subdomain_association(
3170  n_dofs, numbers::invalid_subdomain_id);
3171  std::vector<types::global_dof_index> local_dof_indices;
3172  local_dof_indices.reserve(DoFTools::max_dofs_per_cell(dof_handler));
3173 
3174  // loop over all cells and record which subdomain a DoF belongs to.
3175  // give to the smaller subdomain_id in case it is on an interface
3176  typename DoFHandlerType::active_cell_iterator
3177  cell = dof_handler.begin_active(),
3178  endc = dof_handler.end();
3179  for (; cell != endc; ++cell)
3180  {
3181  // get the owner of the cell; note that we have made sure above
3182  // that all cells are either locally owned or ghosts (not
3183  // artificial), so this call will always yield the true owner
3184  const types::subdomain_id subdomain_id = cell->subdomain_id();
3185  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
3186  local_dof_indices.resize(dofs_per_cell);
3187  cell->get_dof_indices(local_dof_indices);
3188 
3189  // set subdomain ids. if dofs already have their values set then
3190  // they must be on partition interfaces. In that case assign them
3191  // to the processor with the smaller subdomain id.
3192  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3193  if (subdomain_association[local_dof_indices[i]] ==
3194  numbers::invalid_subdomain_id)
3195  subdomain_association[local_dof_indices[i]] = subdomain_id;
3196  else if (subdomain_association[local_dof_indices[i]] >
3197  subdomain_id)
3198  {
3199  subdomain_association[local_dof_indices[i]] = subdomain_id;
3200  }
3201  }
3202 
3203  Assert(std::find(subdomain_association.begin(),
3204  subdomain_association.end(),
3205  numbers::invalid_subdomain_id) ==
3206  subdomain_association.end(),
3207  ExcInternalError());
3208 
3209  Assert(*std::max_element(subdomain_association.begin(),
3210  subdomain_association.end()) < n_procs,
3211  ExcInternalError());
3212 
3213  return subdomain_association;
3214  }
3215 
3216 
3223  template <class DoFHandlerType>
3224  std::vector<types::subdomain_id>
3225  get_dof_level_subdomain_association(
3226  const DoFHandlerType & dof_handler,
3227  const types::global_dof_index n_dofs_on_level,
3228  const unsigned int n_procs,
3229  const unsigned int level)
3230  {
3231  (void)n_procs;
3232  std::vector<types::subdomain_id> level_subdomain_association(
3233  n_dofs_on_level, numbers::invalid_subdomain_id);
3234  std::vector<types::global_dof_index> local_dof_indices;
3235  local_dof_indices.reserve(DoFTools::max_dofs_per_cell(dof_handler));
3236 
3237  // loop over all cells and record which subdomain a DoF belongs to.
3238  // interface goes to proccessor with smaller subdomain id
3239  typename DoFHandlerType::cell_iterator cell =
3240  dof_handler.begin(level),
3241  endc = dof_handler.end(level);
3242  for (; cell != endc; ++cell)
3243  {
3244  // get the owner of the cell; note that we have made sure above
3245  // that all cells are either locally owned or ghosts (not
3246  // artificial), so this call will always yield the true owner
3247  const types::subdomain_id level_subdomain_id =
3248  cell->level_subdomain_id();
3249  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
3250  local_dof_indices.resize(dofs_per_cell);
3251  cell->get_mg_dof_indices(local_dof_indices);
3252 
3253  // set level subdomain ids. if dofs already have their values set
3254  // then they must be on partition interfaces. In that case assign
3255  // them to the processor with the smaller subdomain id.
3256  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3257  if (level_subdomain_association[local_dof_indices[i]] ==
3258  numbers::invalid_subdomain_id)
3259  level_subdomain_association[local_dof_indices[i]] =
3260  level_subdomain_id;
3261  else if (level_subdomain_association[local_dof_indices[i]] >
3262  level_subdomain_id)
3263  {
3264  level_subdomain_association[local_dof_indices[i]] =
3265  level_subdomain_id;
3266  }
3267  }
3268 
3269  Assert(std::find(level_subdomain_association.begin(),
3270  level_subdomain_association.end(),
3271  numbers::invalid_subdomain_id) ==
3272  level_subdomain_association.end(),
3273  ExcInternalError());
3274 
3275  Assert(*std::max_element(level_subdomain_association.begin(),
3276  level_subdomain_association.end()) < n_procs,
3277  ExcInternalError());
3278 
3279  return level_subdomain_association;
3280  }
3281  } // namespace
3282 
3283 
3284 
3285  template <class DoFHandlerType>
3286  NumberCache
3287  ParallelShared<DoFHandlerType>::distribute_dofs() const
3288  {
3289  const unsigned int dim = DoFHandlerType::dimension;
3290  const unsigned int spacedim = DoFHandlerType::space_dimension;
3291 
3292  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3293  (dynamic_cast<
3294  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3295  &this->dof_handler->get_triangulation()));
3296  Assert(tr != nullptr, ExcInternalError());
3297 
3298  const unsigned int n_procs =
3299  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3300 
3301  // If the underlying shared::Tria allows artificial cells,
3302  // then save the current set of subdomain ids, and set
3303  // subdomain ids to the "true" owner of each cell. we later
3304  // restore these flags
3305  std::vector<types::subdomain_id> saved_subdomain_ids;
3306  if (tr->with_artificial_cells())
3307  {
3308  saved_subdomain_ids.resize(tr->n_active_cells());
3309 
3310  const std::vector<types::subdomain_id> &true_subdomain_ids =
3311  tr->get_true_subdomain_ids_of_cells();
3312 
3313  for (const auto &cell : tr->active_cell_iterators())
3314  {
3315  const unsigned int index = cell->active_cell_index();
3316  saved_subdomain_ids[index] = cell->subdomain_id();
3317  cell->set_subdomain_id(true_subdomain_ids[index]);
3318  }
3319  }
3320 
3321  // first let the sequential algorithm do its magic. it is going to
3322  // enumerate DoFs on all cells, regardless of owner
3323  const types::global_dof_index n_initial_dofs =
3324  Implementation::distribute_dofs(numbers::invalid_subdomain_id,
3325  *this->dof_handler);
3326 
3327  const types::global_dof_index n_dofs =
3328  Implementation::unify_dof_indices(*this->dof_handler,
3329  n_initial_dofs,
3330  /*check_validity=*/true);
3331 
3332  // then re-enumerate them based on their subdomain association.
3333  // for this, we first have to identify for each current DoF
3334  // index which subdomain they belong to. ideally, we would
3335  // like to call DoFRenumbering::subdomain_wise(), but
3336  // because the NumberCache of the current DoFHandler is not
3337  // fully set up yet, we can't quite do that. also, that
3338  // function has to deal with other kinds of triangulations as
3339  // well, whereas we here know what kind of triangulation
3340  // we have and can simplify the code accordingly
3341  std::vector<types::global_dof_index> new_dof_indices(
3342  n_dofs, enumeration_dof_index);
3343  {
3344  // first get the association of each dof with a subdomain and
3345  // determine the total number of subdomain ids used
3346  const std::vector<types::subdomain_id> subdomain_association =
3347  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3348 
3349  // then renumber the subdomains by first looking at those belonging
3350  // to subdomain 0, then those of subdomain 1, etc. note that the
3351  // algorithm is stable, i.e. if two dofs i,j have i<j and belong to
3352  // the same subdomain, then they will be in this order also after
3353  // reordering
3354  types::global_dof_index next_free_index = 0;
3355  for (types::subdomain_id subdomain = 0; subdomain < n_procs;
3356  ++subdomain)
3357  for (types::global_dof_index i = 0; i < n_dofs; ++i)
3358  if (subdomain_association[i] == subdomain)
3359  {
3360  Assert(new_dof_indices[i] == enumeration_dof_index,
3361  ExcInternalError());
3362  new_dof_indices[i] = next_free_index;
3363  ++next_free_index;
3364  }
3365 
3366  // we should have numbered all dofs
3367  Assert(next_free_index == n_dofs, ExcInternalError());
3368  Assert(std::find(new_dof_indices.begin(),
3369  new_dof_indices.end(),
3370  enumeration_dof_index) == new_dof_indices.end(),
3371  ExcInternalError());
3372  }
3373  // finally do the renumbering. we can use the sequential
3374  // version of the function because we do things on all
3375  // cells and all cells have their subdomain ids and DoFs
3376  // correctly set
3377  Implementation::renumber_dofs(new_dof_indices,
3378  IndexSet(0),
3379  *this->dof_handler,
3380  /*check_validity=*/true);
3381 
3382  // update the number cache. for this, we first have to find the
3383  // subdomain association for each DoF again following renumbering, from
3384  // which we can then compute the IndexSets of locally owned DoFs for all
3385  // processors. all other fields then follow from this
3386  //
3387  // given the way we enumerate degrees of freedom, the locally owned
3388  // ranges must all be contiguous and consecutive. this makes filling
3389  // the IndexSets cheap. an assertion at the top verifies that this
3390  // assumption is true
3391  const std::vector<types::subdomain_id> subdomain_association =
3392  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3393 
3394  for (unsigned int i = 1; i < n_dofs; ++i)
3395  Assert(subdomain_association[i] >= subdomain_association[i - 1],
3396  ExcInternalError());
3397 
3398  std::vector<IndexSet> locally_owned_dofs_per_processor(
3399  n_procs, IndexSet(n_dofs));
3400  {
3401  // we know that the set of subdomain indices is contiguous from
3402  // the assertion above; find the start and end index for each
3403  // processor, taking into account that sometimes a processor
3404  // may not in fact have any DoFs at all. we do the latter
3405  // by just identifying contiguous ranges of subdomain_ids
3406  // and filling IndexSets for those subdomains; subdomains
3407  // that don't appear will lead to IndexSets that are simply
3408  // never touched and remain empty as initialized above.
3409  unsigned int start_index = 0;
3410  unsigned int end_index = 0;
3411  while (start_index < n_dofs)
3412  {
3413  while ((end_index) < n_dofs &&
3414  (subdomain_association[end_index] ==
3415  subdomain_association[start_index]))
3416  ++end_index;
3417 
3418  // we've now identified a range of same indices. set that
3419  // range in the corresponding IndexSet
3420  if (end_index > start_index)
3421  {
3422  const unsigned int subdomain_owner =
3423  subdomain_association[start_index];
3424  locally_owned_dofs_per_processor[subdomain_owner].add_range(
3425  start_index, end_index);
3426  }
3427 
3428  // then move on to thinking about the next range
3429  start_index = end_index;
3430  }
3431  }
3432 
3433  // finally, restore current subdomain ids
3434  if (tr->with_artificial_cells())
3435  for (const auto &cell : tr->active_cell_iterators())
3436  cell->set_subdomain_id(
3437  saved_subdomain_ids[cell->active_cell_index()]);
3438 
3439  // return a NumberCache object made up from the sets of locally
3440  // owned DoFs
3441  return NumberCache(
3442  locally_owned_dofs_per_processor,
3443  this->dof_handler->get_triangulation().locally_owned_subdomain());
3444  }
3445 
3446 
3447 
3448  template <class DoFHandlerType>
3449  std::vector<NumberCache>
3450  ParallelShared<DoFHandlerType>::distribute_mg_dofs() const
3451  {
3452  const unsigned int dim = DoFHandlerType::dimension;
3453  const unsigned int spacedim = DoFHandlerType::space_dimension;
3454 
3455  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3456  (dynamic_cast<
3457  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3458  &this->dof_handler->get_triangulation()));
3459  Assert(tr != nullptr, ExcInternalError());
3460 
3461  const unsigned int n_procs =
3462  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3463  const unsigned int n_levels = tr->n_global_levels();
3464 
3465  std::vector<NumberCache> number_caches;
3466  number_caches.reserve(n_levels);
3467 
3468  // We create an index set for each level
3469  for (unsigned int lvl = 0; lvl < n_levels; ++lvl)
3470  {
3471  // If the underlying shared::Tria allows artificial cells,
3472  // then save the current set of level subdomain ids, and set
3473  // subdomain ids to the "true" owner of each cell. we later
3474  // restore these flags
3475  // Note: "allows_artificial_cells" is currently enforced for
3476  // MG computations.
3477  std::vector<types::subdomain_id> saved_level_subdomain_ids;
3478  saved_level_subdomain_ids.resize(tr->n_cells(lvl));
3479  {
3480  typename ::parallel::shared::Triangulation<dim, spacedim>::
3481  cell_iterator cell =
3482  this->dof_handler->get_triangulation().begin(
3483  lvl),
3484  endc =
3485  this->dof_handler->get_triangulation().end(lvl);
3486 
3487  const std::vector<types::subdomain_id> &true_level_subdomain_ids =
3488  tr->get_true_level_subdomain_ids_of_cells(lvl);
3489 
3490  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3491  {
3492  saved_level_subdomain_ids[index] = cell->level_subdomain_id();
3493  cell->set_level_subdomain_id(true_level_subdomain_ids[index]);
3494  }
3495  }
3496 
3497  // Next let the sequential algorithm do its magic. it is going to
3498  // enumerate DoFs on all cells on the given level, regardless of
3499  // owner
3500  const types::global_dof_index n_dofs_on_level =
3501  Implementation::distribute_dofs_on_level(
3502  numbers::invalid_subdomain_id, *this->dof_handler, lvl);
3503 
3504  // then re-enumerate them based on their level subdomain
3505  // association. for this, we first have to identify for each current
3506  // DoF index which subdomain they belong to. ideally, we would like
3507  // to call DoFRenumbering::subdomain_wise(), but because the
3508  // NumberCache of the current DoFHandler is not fully set up yet, we
3509  // can't quite do that. also, that function has to deal with other
3510  // kinds of triangulations as well, whereas we here know what kind
3511  // of triangulation we have and can simplify the code accordingly
3512  std::vector<types::global_dof_index> new_dof_indices(
3513  n_dofs_on_level, numbers::invalid_dof_index);
3514  {
3515  // first get the association of each dof with a subdomain and
3516  // determine the total number of subdomain ids used
3517  const std::vector<types::subdomain_id>
3518  level_subdomain_association =
3519  get_dof_level_subdomain_association(*this->dof_handler,
3520  n_dofs_on_level,
3521  n_procs,
3522  lvl);
3523 
3524  // then renumber the subdomains by first looking at those
3525  // belonging to subdomain 0, then those of subdomain 1, etc. note
3526  // that the algorithm is stable, i.e. if two dofs i,j have i<j and
3527  // belong to the same subdomain, then they will be in this order
3528  // also after reordering
3529  types::global_dof_index next_free_index = 0;
3530  for (types::subdomain_id level_subdomain = 0;
3531  level_subdomain < n_procs;
3532  ++level_subdomain)
3533  for (types::global_dof_index i = 0; i < n_dofs_on_level; ++i)
3534  if (level_subdomain_association[i] == level_subdomain)
3535  {
3536  Assert(new_dof_indices[i] == numbers::invalid_dof_index,
3537  ExcInternalError());
3538  new_dof_indices[i] = next_free_index;
3539  ++next_free_index;
3540  }
3541 
3542  // we should have numbered all dofs
3543  Assert(next_free_index == n_dofs_on_level, ExcInternalError());
3544  Assert(std::find(new_dof_indices.begin(),
3545  new_dof_indices.end(),
3546  numbers::invalid_dof_index) ==
3547  new_dof_indices.end(),
3548  ExcInternalError());
3549  }
3550 
3551  // finally do the renumbering. we can use the sequential
3552  // version of the function because we do things on all
3553  // cells and all cells have their subdomain ids and DoFs
3554  // correctly set
3555  Implementation::renumber_mg_dofs(
3556  new_dof_indices, IndexSet(0), *this->dof_handler, lvl, true);
3557 
3558  // update the number cache. for this, we first have to find the
3559  // level subdomain association for each DoF again following
3560  // renumbering, from which we can then compute the IndexSets of
3561  // locally owned DoFs for all processors. all other fields then
3562  // follow from this
3563  //
3564  // given the way we enumerate degrees of freedom, the locally owned
3565  // ranges must all be contiguous and consecutive. this makes filling
3566  // the IndexSets cheap. an assertion at the top verifies that this
3567  // assumption is true
3568  const std::vector<types::subdomain_id> level_subdomain_association =
3569  get_dof_level_subdomain_association(*this->dof_handler,
3570  n_dofs_on_level,
3571  n_procs,
3572  lvl);
3573 
3574  for (unsigned int i = 1; i < n_dofs_on_level; ++i)
3575  Assert(level_subdomain_association[i] >=
3576  level_subdomain_association[i - 1],
3577  ExcInternalError());
3578 
3579  std::vector<IndexSet> locally_owned_dofs_per_processor(
3580  n_procs, IndexSet(n_dofs_on_level));
3581  {
3582  // we know that the set of subdomain indices is contiguous from
3583  // the assertion above; find the start and end index for each
3584  // processor, taking into account that sometimes a processor
3585  // may not in fact have any DoFs at all. we do the latter
3586  // by just identifying contiguous ranges of level_subdomain_ids
3587  // and filling IndexSets for those subdomains; subdomains
3588  // that don't appear will lead to IndexSets that are simply
3589  // never touched and remain empty as initialized above.
3590  unsigned int start_index = 0;
3591  unsigned int end_index = 0;
3592  while (start_index < n_dofs_on_level)
3593  {
3594  while ((end_index) < n_dofs_on_level &&
3595  (level_subdomain_association[end_index] ==
3596  level_subdomain_association[start_index]))
3597  ++end_index;
3598 
3599  // we've now identified a range of same indices. set that
3600  // range in the corresponding IndexSet
3601  if (end_index > start_index)
3602  {
3603  const unsigned int level_subdomain_owner =
3604  level_subdomain_association[start_index];
3605  locally_owned_dofs_per_processor[level_subdomain_owner]
3606  .add_range(start_index, end_index);
3607  }
3608 
3609  // then move on to thinking about the next range
3610  start_index = end_index;
3611  }
3612  }
3613 
3614  // finally, restore current level subdomain ids
3615  {
3616  typename ::parallel::shared::Triangulation<dim, spacedim>::
3617  cell_iterator cell =
3618  this->dof_handler->get_triangulation().begin(
3619  lvl),
3620  endc =
3621  this->dof_handler->get_triangulation().end(lvl);
3622 
3623  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3624  cell->set_level_subdomain_id(saved_level_subdomain_ids[index]);
3625 
3626  // add NumberCache for current level
3627  number_caches.emplace_back(
3628  NumberCache(locally_owned_dofs_per_processor,
3629  this->dof_handler->get_triangulation()
3630  .locally_owned_subdomain()));
3631  }
3632  }
3633 
3634  return number_caches;
3635  }
3636 
3637 
3638 
3639  template <class DoFHandlerType>
3640  NumberCache
3641  ParallelShared<DoFHandlerType>::renumber_dofs(
3642  const std::vector<types::global_dof_index> &new_numbers) const
3643  {
3644 #ifndef DEAL_II_WITH_MPI
3645  (void)new_numbers;
3646  Assert(false, ExcNotImplemented());
3647  return NumberCache();
3648 #else
3649  const unsigned int dim = DoFHandlerType::dimension;
3650  const unsigned int spacedim = DoFHandlerType::space_dimension;
3651 
3652  // Similar to distribute_dofs() we need to have a special treatment in
3653  // case artificial cells are present.
3654  const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3655  (dynamic_cast<
3656  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3657  &this->dof_handler->get_triangulation()));
3658  Assert(tr != nullptr, ExcInternalError());
3659 
3660  typename ::parallel::shared::Triangulation<dim, spacedim>::
3661  active_cell_iterator
3662  cell = this->dof_handler->get_triangulation().begin_active(),
3663  endc = this->dof_handler->get_triangulation().end();
3664  std::vector<types::subdomain_id> current_subdomain_ids(
3665  tr->n_active_cells());
3666  const std::vector<types::subdomain_id> &true_subdomain_ids =
3667  tr->get_true_subdomain_ids_of_cells();
3668  if (tr->with_artificial_cells())
3669  for (unsigned int index = 0; cell != endc; cell++, index++)
3670  {
3671  current_subdomain_ids[index] = cell->subdomain_id();
3672  cell->set_subdomain_id(true_subdomain_ids[index]);
3673  }
3674 
3675  std::vector<types::global_dof_index> global_gathered_numbers(
3676  this->dof_handler->n_dofs(), 0);
3677  // as we call DoFRenumbering::subdomain_wise (*dof_handler) from
3678  // distribute_dofs(), we need to support sequential-like input.
3679  // Distributed-like input from, for example, component_wise renumbering
3680  // is also supported.
3681  if (new_numbers.size() == this->dof_handler->n_dofs())
3682  {
3683  global_gathered_numbers = new_numbers;
3684  }
3685  else
3686  {
3687  Assert(new_numbers.size() ==
3688  this->dof_handler->locally_owned_dofs().n_elements(),
3689  ExcInternalError());
3690  const unsigned int n_cpu =
3691  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3692  std::vector<types::global_dof_index> gathered_new_numbers(
3693  this->dof_handler->n_dofs(), 0);
3694  Assert(Utilities::MPI::this_mpi_process(tr->get_communicator()) ==
3695  this->dof_handler->get_triangulation()
3696  .locally_owned_subdomain(),
3697  ExcInternalError())
3698 
3699  // gather new numbers among processors into one vector
3700  {
3701  std::vector<types::global_dof_index> new_numbers_copy(
3702  new_numbers);
3703 
3704  // store the number of elements that are to be received from each
3705  // process
3706  std::vector<int> rcounts(n_cpu);
3707 
3708  types::global_dof_index shift = 0;
3709  // set rcounts based on new_numbers:
3710  int cur_count = new_numbers_copy.size();
3711  int ierr = MPI_Allgather(&cur_count,
3712  1,
3713  MPI_INT,
3714  rcounts.data(),
3715  1,
3716  MPI_INT,
3717  tr->get_communicator());
3718  AssertThrowMPI(ierr);
3719 
3720  // compute the displacements (relative to recvbuf)
3721  // at which to place the incoming data from process i
3722  std::vector<int> displacements(n_cpu);
3723  for (unsigned int i = 0; i < n_cpu; i++)
3724  {
3725  displacements[i] = shift;
3726  shift += rcounts[i];
3727  }
3728  Assert(new_numbers_copy.size() ==
3729  static_cast<unsigned int>(
3730  rcounts[Utilities::MPI::this_mpi_process(
3731  tr->get_communicator())]),
3732  ExcInternalError());
3733  ierr = MPI_Allgatherv(new_numbers_copy.data(),
3734  new_numbers_copy.size(),
3735  DEAL_II_DOF_INDEX_MPI_TYPE,
3736  gathered_new_numbers.data(),
3737  rcounts.data(),
3738  displacements.data(),
3739  DEAL_II_DOF_INDEX_MPI_TYPE,
3740  tr->get_communicator());
3741  AssertThrowMPI(ierr);
3742  }
3743 
3744  // put new numbers according to the current
3745  // locally_owned_dofs_per_processor IndexSets
3746  types::global_dof_index shift = 0;
3747  // flag_1 and flag_2 are
3748  // used to control that there is a
3749  // one-to-one relation between old and new DoFs.
3750  std::vector<unsigned int> flag_1(this->dof_handler->n_dofs(), 0);
3751  std::vector<unsigned int> flag_2(this->dof_handler->n_dofs(), 0);
3752  for (unsigned int i = 0; i < n_cpu; i++)
3753  {
3754  const IndexSet iset =
3755  this->dof_handler->locally_owned_dofs_per_processor()[i];
3756  for (types::global_dof_index ind = 0; ind < iset.n_elements();
3757  ind++)
3758  {
3759  const types::global_dof_index target =
3760  iset.nth_index_in_set(ind);
3761  const types::global_dof_index value =
3762  gathered_new_numbers[shift + ind];
3763  Assert(target < this->dof_handler->n_dofs(),
3764  ExcInternalError());
3765  Assert(value < this->dof_handler->n_dofs(),
3766  ExcInternalError());
3767  global_gathered_numbers[target] = value;
3768  flag_1[target]++;
3769  flag_2[value]++;
3770  }
3771  shift += iset.n_elements();
3772  }
3773 
3774  Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1,
3775  ExcInternalError());
3776  Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1,
3777  ExcInternalError());
3778  Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1,
3779  ExcInternalError());
3780  Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1,
3781  ExcInternalError());
3782  }
3783 
3784  // let the sequential algorithm do its magic; ignore the
3785  // return type, but reconstruct the number cache based on
3786  // which DoFs each process owns
3787  Implementation::renumber_dofs(global_gathered_numbers,
3788  IndexSet(0),
3789  *this->dof_handler,
3790  /*check_validity=*/true);
3791 
3792  const NumberCache number_cache(
3793  DoFTools::locally_owned_dofs_per_subdomain(*this->dof_handler),
3794  this->dof_handler->get_triangulation().locally_owned_subdomain());
3795 
3796  // restore artificial cells
3797  cell = tr->begin_active();
3798  if (tr->with_artificial_cells())
3799  for (unsigned int index = 0; cell != endc; cell++, index++)
3800  cell->set_subdomain_id(current_subdomain_ids[index]);
3801 
3802  return number_cache;
3803 #endif
3804  }
3805 
3806 
3807 
3808  template <class DoFHandlerType>
3809  NumberCache
3810  ParallelShared<DoFHandlerType>::renumber_mg_dofs(
3811  const unsigned int /*level*/,
3812  const std::vector<types::global_dof_index> & /*new_numbers*/) const
3813  {
3814  // multigrid is not currently implemented for shared triangulations
3815  Assert(false, ExcNotImplemented());
3816 
3817  return {};
3818  }
3819 
3820 
3821 
3822  /* --------------------- class ParallelDistributed ---------------- */
3823 
3824 #ifdef DEAL_II_WITH_MPI
3825 
3826  namespace
3827  {
3828  template <int dim, int spacedim>
3829  void
3830  get_mg_dofindices_recursively(
3831  const ::parallel::DistributedTriangulationBase<dim, spacedim>
3832  &tria,
3833  const typename DoFHandler<dim, spacedim>::level_cell_iterator
3834  & dealii_cell,
3835  const typename CellId::binary_type & quadrant,
3836  std::vector<::types::global_dof_index> &dof_numbers_and_indices)
3837  {
3838  if (dealii_cell->id() == CellId(quadrant))
3839  {
3840  // why would somebody request a cell that is not ours?
3841  Assert(dealii_cell->level_subdomain_id() ==
3842  tria.locally_owned_subdomain(),
3843  ExcInternalError());
3844 
3845  std::vector<::types::global_dof_index> local_dof_indices(
3846  dealii_cell->get_fe().dofs_per_cell);
3847  dealii_cell->get_mg_dof_indices(local_dof_indices);
3848 
3849  dof_numbers_and_indices.push_back(
3850  dealii_cell->get_fe().dofs_per_cell);
3851  dof_numbers_and_indices.insert(dof_numbers_and_indices.end(),
3852  local_dof_indices.begin(),
3853  local_dof_indices.end());
3854  return; // we are done
3855  }
3856 
3857  if (dealii_cell->is_active())
3858  return;
3859 
3860  if (!dealii_cell->id().is_ancestor_of(CellId(quadrant)))
3861  return;
3862 
3863  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
3864  ++c)
3865  get_mg_dofindices_recursively<dim, spacedim>(
3866  tria, dealii_cell->child(c), quadrant, dof_numbers_and_indices);
3867  }
3868 
3869 
3870 
3871  template <int dim, int spacedim>
3872  void
3873  find_marked_mg_ghost_cells_recursively(
3874  const typename ::parallel::
3875  DistributedTriangulationBase<dim, spacedim> &tria,
3876  const unsigned int tree_index,
3877  const typename DoFHandler<dim, spacedim>::level_cell_iterator
3878  &dealii_cell,
3879  std::map<
3880  ::types::subdomain_id,
3881  std::vector<std::pair<unsigned int, typename CellId::binary_type>>>
3882  &neighbor_cell_list)
3883  {
3884  // recurse...
3885  if (dealii_cell->has_children())
3886  {
3887  for (unsigned int c = 0;
3888  c < GeometryInfo<dim>::max_children_per_cell;
3889  ++c)
3890  find_marked_mg_ghost_cells_recursively<dim, spacedim>(
3891  tria, tree_index, dealii_cell->child(c), neighbor_cell_list);
3892  }
3893 
3894  if (dealii_cell->user_flag_set() &&
3895  dealii_cell->level_subdomain_id() !=
3896  tria.locally_owned_subdomain())
3897  {
3898  neighbor_cell_list[dealii_cell->level_subdomain_id()]
3899  .emplace_back(tree_index,
3900  dealii_cell->id().template to_binary<spacedim>());
3901  }
3902  }
3903 
3904 
3905 
3906  template <int dim, int spacedim>
3907  void
3908  set_mg_dofindices_recursively(
3909  const ::parallel::DistributedTriangulationBase<dim, spacedim>
3910  &tria,
3911  const typename DoFHandler<dim, spacedim>::level_cell_iterator
3912  & dealii_cell,
3913  const typename CellId::binary_type &quadrant,
3914  ::types::global_dof_index * dofs)
3915  {
3916  if (dealii_cell->id() == CellId(quadrant))
3917  {
3918  Assert(dealii_cell->level_subdomain_id() !=
3919  ::numbers::artificial_subdomain_id,
3920  ExcInternalError());
3921 
3922  // update dof indices of cell
3923  std::vector<::types::global_dof_index> dof_indices(
3924  dealii_cell->get_fe().dofs_per_cell);
3925  dealii_cell->get_mg_dof_indices(dof_indices);
3926 
3927  bool complete = true;
3928  for (unsigned int i = 0; i < dof_indices.size(); ++i)
3929  if (dofs[i] != numbers::invalid_dof_index)
3930  {
3931  Assert((dof_indices[i] == (numbers::invalid_dof_index)) ||
3932  (dof_indices[i] == dofs[i]),
3933  ExcInternalError());
3934  dof_indices[i] = dofs[i];
3935  }
3936  else
3937  complete = false;
3938 
3939  if (!complete)
3940  const_cast<
3941  typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
3942  dealii_cell)
3943  ->set_user_flag();
3944  else
3945  const_cast<
3946  typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
3947  dealii_cell)
3948  ->clear_user_flag();
3949 
3950  const_cast<
3951  typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
3952  dealii_cell)
3953  ->set_mg_dof_indices(dof_indices);
3954  return;
3955  }
3956 
3957  if (dealii_cell->is_active())
3958  return;
3959 
3960  if (!dealii_cell->id().is_ancestor_of(CellId(quadrant)))
3961  return;
3962 
3963  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
3964  ++c)
3965  set_mg_dofindices_recursively<dim, spacedim>(tria,
3966  dealii_cell->child(c),
3967  quadrant,
3968  dofs);
3969  }
3970 
3971 
3972 
3973  template <int dim, int spacedim, class DoFHandlerType>
3974  void
3975  communicate_mg_ghost_cells(
3976  const typename ::parallel::
3977  DistributedTriangulationBase<dim, spacedim> &tria,
3978  DoFHandlerType & dof_handler)
3979  {
3980  using QuadrantBufferType =
3981  std::vector<std::pair<unsigned int, typename CellId::binary_type>>;
3982  // build list of cells to request for each neighbor
3983  std::set<::types::subdomain_id> level_ghost_owners =
3984  tria.level_ghost_owners();
3985  std::map<::types::subdomain_id, QuadrantBufferType>
3986  neighbor_cell_list;
3987  for (const auto level_ghost_owner : level_ghost_owners)
3988  neighbor_cell_list[level_ghost_owner] = {};
3989 
3990  for (typename DoFHandlerType::level_cell_iterator cell =
3991  dof_handler.begin(0);
3992  cell != dof_handler.end(0);
3993  ++cell)
3994  {
3995  types::coarse_cell_id coarse_cell_id = 0;
3996  try
3997  {
3998  coarse_cell_id = cell->id().get_coarse_cell_id();
3999  }
4000  catch (...)
4001  {
4002  // In the case of parallel::fullydistributed::Triangulation,
4003  // a dummy cell throws an exception which is caught here.
4004  // We ignore this cell here.
4005  continue;
4006  };
4007 
4008  find_marked_mg_ghost_cells_recursively<dim, spacedim>(
4009  tria, coarse_cell_id, cell, neighbor_cell_list);
4010  }
4011  Assert(level_ghost_owners.size() == neighbor_cell_list.size(),
4012  ExcInternalError());
4013 
4014 
4015  // Before sending & receiving, make sure we protect this section with
4016  // a mutex:
4017  static Utilities::MPI::CollectiveMutex mutex;
4018  Utilities::MPI::CollectiveMutex::ScopedLock lock(
4019  mutex, tria.get_communicator());
4020 
4021  const int mpi_tag = Utilities::MPI::internal::Tags::
4022  dofhandler_communicate_mg_ghost_cells;
4023  const int mpi_tag_reply = Utilities::MPI::internal::Tags::
4024  dofhandler_communicate_mg_ghost_cells_reply;
4025 
4026  //* send our requests:
4027  std::vector<MPI_Request> requests(level_ghost_owners.size());
4028  {
4029  unsigned int idx = 0;
4030  for (const auto &it : neighbor_cell_list)
4031  {
4032  // send the data about the relevant cells
4033  const int ierr =
4034  MPI_Isend(it.second.data(),
4035  it.second.size() * sizeof(it.second[0]),
4036  MPI_BYTE,
4037  it.first,
4038  mpi_tag,
4039  tria.get_communicator(),
4040  &requests[idx]);
4041  AssertThrowMPI(ierr);
4042  ++idx;
4043  }
4044  }
4045 
4046  //* receive requests and reply with the ghost indices
4047  std::vector<QuadrantBufferType> quadrant_data_to_send(
4048  level_ghost_owners.size());
4049  std::vector<std::vector<types::global_dof_index>>
4050  send_dof_numbers_and_indices(level_ghost_owners.size());
4051  std::vector<MPI_Request> reply_requests(level_ghost_owners.size());
4052 
4053  for (unsigned int idx = 0; idx < level_ghost_owners.size(); ++idx)
4054  {
4055  MPI_Status status;
4056  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4057  mpi_tag,
4058  tria.get_communicator(),
4059  &status);
4060  AssertThrowMPI(ierr);
4061 
4062  int len;
4063  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4064  AssertThrowMPI(ierr);
4065  Assert(len % sizeof(quadrant_data_to_send[idx][0]) == 0,
4066  ExcInternalError());
4067 
4068  const unsigned int n_cells =
4069  len / sizeof(quadrant_data_to_send[idx][0]);
4070  quadrant_data_to_send[idx].resize(n_cells);
4071 
4072  ierr = MPI_Recv(quadrant_data_to_send[idx].data(),
4073  len,
4074  MPI_BYTE,
4075  status.MPI_SOURCE,
4076  status.MPI_TAG,
4077  tria.get_communicator(),
4078  &status);
4079  AssertThrowMPI(ierr);
4080 
4081  // store the dof indices for each cell
4082  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells);
4083  ++c)
4084  {
4085  const auto temp =
4086  CellId(quadrant_data_to_send[idx][c].first, 0, nullptr)
4087  .to_cell(tria);
4088 
4089  typename DoFHandlerType::level_cell_iterator cell(
4090  &dof_handler.get_triangulation(),
4091  0,
4092  temp->index(),
4093  &dof_handler);
4094 
4095  get_mg_dofindices_recursively<dim, spacedim>(
4096  tria,
4097  cell,
4098  quadrant_data_to_send[idx][c].second,
4099  send_dof_numbers_and_indices[idx]);
4100  }
4101 
4102  // send reply
4103  ierr = MPI_Isend(send_dof_numbers_and_indices[idx].data(),
4104  send_dof_numbers_and_indices[idx].size(),
4105  DEAL_II_DOF_INDEX_MPI_TYPE,
4106  status.MPI_SOURCE,
4107  mpi_tag_reply,
4108  tria.get_communicator(),
4109  &reply_requests[idx]);
4110  AssertThrowMPI(ierr);
4111  }
4112 
4113  //* finally receive the replies
4114  for (unsigned int idx = 0; idx < level_ghost_owners.size(); ++idx)
4115  {
4116  MPI_Status status;
4117  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4118  mpi_tag_reply,
4119  tria.get_communicator(),
4120  &status);
4121  AssertThrowMPI(ierr);
4122  int len;
4123  ierr = MPI_Get_count(&status, DEAL_II_DOF_INDEX_MPI_TYPE, &len);
4124  const QuadrantBufferType &quadrants =
4125  neighbor_cell_list[status.MPI_SOURCE];
4126  AssertThrowMPI(ierr);
4127  Assert((len > 0 && !quadrants.empty()) ||
4128  (len == 0 && quadrants.empty()),
4129  ExcInternalError());
4130  std::vector<types::global_dof_index>
4131  receive_dof_numbers_and_indices(len);
4132 
4133  ierr = MPI_Recv(receive_dof_numbers_and_indices.data(),
4134  len,
4135  DEAL_II_DOF_INDEX_MPI_TYPE,
4136  status.MPI_SOURCE,
4137  status.MPI_TAG,
4138  tria.get_communicator(),
4139  &status);
4140  AssertThrowMPI(ierr);
4141 
4142  // set the dof indices for each cell
4143  ::types::global_dof_index *dofs =
4144  receive_dof_numbers_and_indices.data();
4145  for (const auto &it : quadrants)
4146  {
4147  const auto temp = CellId(it.first, 0, nullptr).to_cell(tria);
4148 
4149  typename DoFHandlerType::level_cell_iterator cell(
4150  &tria, 0, temp->index(), &dof_handler);
4151 
4152  Assert(cell->get_fe().dofs_per_cell == dofs[0],
4153  ExcInternalError());
4154 
4155  set_mg_dofindices_recursively<dim, spacedim>(tria,
4156  cell,
4157  it.second,
4158  dofs + 1);
4159  dofs += 1 + dofs[0];
4160  }
4161  Assert(dofs == receive_dof_numbers_and_indices.data() +
4162  receive_dof_numbers_and_indices.size(),
4163  ExcInternalError());
4164  }
4165 
4166  // complete all sends, so that we can safely destroy the
4167  // buffers.
4168  if (requests.size() > 0)
4169  {
4170  const int ierr = MPI_Waitall(requests.size(),
4171  requests.data(),
4172  MPI_STATUSES_IGNORE);
4173  AssertThrowMPI(ierr);
4174  }
4175  if (reply_requests.size() > 0)
4176  {
4177  const int ierr = MPI_Waitall(reply_requests.size(),
4178  reply_requests.data(),
4179  MPI_STATUSES_IGNORE);
4180  AssertThrowMPI(ierr);
4181  }
4182  }
4183 
4184 
4185 
4186  template <int spacedim>
4187  void
4188  communicate_mg_ghost_cells(const typename ::parallel::
4189  distributed::Triangulation<1, spacedim> &,
4190  DoFHandler<1, spacedim> &)
4191  {
4192  Assert(false, ExcNotImplemented());
4193  }
4194 
4195 
4196 
4197  template <int spacedim>
4198  void
4199  communicate_mg_ghost_cells(const typename ::parallel::
4200  distributed::Triangulation<1, spacedim> &,
4201  hp::DoFHandler<1, spacedim> &)
4202  {
4203  Assert(false, ExcNotImplemented());
4204  }
4205 
4206 
4207 
4226  template <class DoFHandlerType>
4227  void
4228  communicate_dof_indices_on_marked_cells(
4229  const DoFHandlerType &dof_handler,
4230  const std::map<unsigned int, std::set<::types::subdomain_id>> &)
4231  {
4232 # ifndef DEAL_II_WITH_MPI
4233  (void)vertices_with_ghost_neighbors;
4234  Assert(false, ExcNotImplemented());
4235 # else
4236  const unsigned int dim = DoFHandlerType::dimension;
4237  const unsigned int spacedim = DoFHandlerType::space_dimension;
4238 
4239  // define functions that pack data on cells that are ghost cells
4240  // somewhere else, and unpack data on cells where we get information
4241  // from elsewhere
4242  auto pack =
4243  [](const typename DoFHandlerType::active_cell_iterator &cell)
4244  -> std_cxx17::optional<std::vector<types::global_dof_index>> {
4245  Assert(cell->is_locally_owned(), ExcInternalError());
4246 
4247  // first see whether we need to do anything at all on this cell.
4248  // this is determined by whether the user_flag is set on the
4249  // cell that indicates that the *complete* set of DoF indices
4250  // has not been sent
4251  if (cell->user_flag_set())
4252  {
4253  // get dof indices for the current cell
4254  std::vector<types::global_dof_index> local_dof_indices(
4255  cell->get_fe().dofs_per_cell);
4256  cell->get_dof_indices(local_dof_indices);
4257 
4258  // now see if there are dof indices that were previously
4259  // unknown. this can only happen in phase 1, and in
4260  // that case we know that the user flag must have been set
4261  //
4262  // in any case, if the cell *is* complete, we do not
4263  // need to send the data any more in the next phase. indicate
4264  // this by removing the user flag
4265  if (std::find(local_dof_indices.begin(),
4266  local_dof_indices.end(),
4267  numbers::invalid_dof_index) !=
4268  local_dof_indices.end())
4269  {
4270  Assert(cell->user_flag_set(), ExcInternalError());
4271  }
4272  else
4273  cell->clear_user_flag();
4274 
4275  return local_dof_indices;
4276  }
4277  else
4278  {
4279  // the fact that the user flag wasn't set means that there is
4280  // nothing we need to send that hasn't been sent so far.
4281  // so return an empty array, but also verify that indeed
4282  // the cell is complete
4283 # ifdef DEBUG
4284  std::vector<types::global_dof_index> local_dof_indices(
4285  cell->get_fe().dofs_per_cell);
4286  cell->get_dof_indices(local_dof_indices);
4287 
4288  const bool is_complete =
4289  (std::find(local_dof_indices.begin(),
4290  local_dof_indices.end(),
4291  numbers::invalid_dof_index) ==
4292  local_dof_indices.end());
4293  Assert(is_complete, ExcInternalError());
4294 # endif
4295  return std_cxx17::optional<
4296  std::vector<types::global_dof_index>>();
4297  }
4298  };
4299 
4300  auto unpack =
4301  [](const typename DoFHandlerType::active_cell_iterator &cell,
4302  const std::vector<types::global_dof_index> &received_dof_indices)
4303  -> void {
4304  // this function should only be called on ghost cells, and
4305  // on top of that, only on cells that have not been
4306  // completed -- which we indicate via the user flag.
4307  // check both
4308  Assert(cell->is_ghost(), ExcInternalError());
4309  Assert(cell->user_flag_set(), ExcInternalError());
4310 
4311  // if we just got an incomplete array of DoF indices, then we must
4312  // be in the first ghost exchange and the user flag must have been
4313  // set. we tested that already above.
4314  //
4315  // if we did get a complete array, then we may be in the first
4316  // or second ghost exchange, but in any case we need not exchange
4317  // another time. so delete the user flag
4318  const bool is_complete = (std::find(received_dof_indices.begin(),
4319  received_dof_indices.end(),
4320  numbers::invalid_dof_index) ==
4321  received_dof_indices.end());
4322  if (is_complete)
4323  cell->clear_user_flag();
4324 
4325  // in any case, set the DoF indices on this cell. some
4326  // of the ones we received may still be invalid because
4327  // the sending processor did not know them yet, so we
4328  // need to merge the ones we get with those that are
4329  // already set here and may have already been known. for
4330  // those that we already know *and* get, they must obviously
4331  // agree
4332  //
4333  // before getting the local dof indices, we need to update the
4334  // cell dof indices cache because we may have set dof indices
4335  // on a neighboring ghost cell before this one, which may have
4336  // affected the dof indices we know about the current cell
4337  std::vector<types::global_dof_index> local_dof_indices(
4338  cell->get_fe().dofs_per_cell);
4339  cell->update_cell_dof_indices_cache();
4340  cell->get_dof_indices(local_dof_indices);
4341 
4342  for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
4343  if (local_dof_indices[i] == numbers::invalid_dof_index)
4344  local_dof_indices[i] = received_dof_indices[i];
4345  else
4346  // we already know the dof index. check that there
4347  // is no conflict
4348  Assert((received_dof_indices[i] ==
4349  numbers::invalid_dof_index) ||
4350  (received_dof_indices[i] == local_dof_indices[i]),
4351  ExcInternalError());
4352 
4353  const_cast<typename DoFHandlerType::active_cell_iterator &>(cell)
4354  ->set_dof_indices(local_dof_indices);
4355  };
4356 
4357  GridTools::exchange_cell_data_to_ghosts<
4358  std::vector<types::global_dof_index>,
4359  DoFHandlerType>(dof_handler, pack, unpack);
4360 
4361  // finally update the cell DoF indices caches to make sure
4362  // our internal data structures are consistent
4363  update_all_active_cell_dof_indices_caches(dof_handler);
4364 
4365 
4366  // have a barrier so that sends between two calls to this
4367  // function are not mixed up.
4368  //
4369  // this is necessary because above we just see if there are
4370  // messages and then receive them, without discriminating
4371  // where they come from and whether they were sent in phase
4372  // 1 or 2 (the function is called twice in a row). the need
4373  // for a global communication step like this barrier could
4374  // be avoided by receiving messages specifically from those
4375  // processors from which we expect messages, and by using
4376  // different tags for phase 1 and 2, but the cost of a
4377  // barrier is negligible compared to everything else we do
4378  // here
4379  if (const auto *triangulation =
4380  dynamic_cast<const ::parallel::
4381  DistributedTriangulationBase<dim, spacedim> *>(
4382  &dof_handler.get_triangulation()))
4383  {
4384  const int ierr = MPI_Barrier(triangulation->get_communicator());
4385  AssertThrowMPI(ierr);
4386  }
4387  else
4388  {
4389  Assert(false,
4390  ExcMessage(
4391  "The function communicate_dof_indices_on_marked_cells() "
4392  "only works with parallel distributed triangulations."));
4393  }
4394 # endif
4395  }
4396 
4397 
4398 
4399  } // namespace
4400 
4401 #endif // DEAL_II_WITH_MPI
4402 
4403 
4404 
4405  template <class DoFHandlerType>
4406  ParallelDistributed<DoFHandlerType>::ParallelDistributed(
4407  DoFHandlerType &dof_handler)
4408  : dof_handler(&dof_handler)
4409  {}
4410 
4411 
4412 
4413  template <class DoFHandlerType>
4414  NumberCache
4415  ParallelDistributed<DoFHandlerType>::distribute_dofs() const
4416  {
4417 #ifndef DEAL_II_WITH_MPI
4418  Assert(false, ExcNotImplemented());
4419  return NumberCache();
4420 #else
4421  const unsigned int dim = DoFHandlerType::dimension;
4422  const unsigned int spacedim = DoFHandlerType::space_dimension;
4423 
4424  ::parallel::DistributedTriangulationBase<dim, spacedim>
4425  *triangulation =
4426  (dynamic_cast<
4427  ::parallel::DistributedTriangulationBase<dim, spacedim> *>(
4428  const_cast<::Triangulation<dim, spacedim> *>(
4429  &dof_handler->get_triangulation())));
4430  Assert(triangulation != nullptr, ExcInternalError());
4431 
4432  const types::subdomain_id subdomain_id =
4433  triangulation->locally_owned_subdomain();
4434 
4435 
4436  /*
4437  The following algorithm has a number of stages that are all
4438  documented in the paper that describes the parallel::distributed
4439  functionality:
4440 
4441  1/ locally enumerate dofs on locally owned cells
4442  2/ eliminate dof duplicates on all cells.
4443  un-numerate those that are on interfaces with ghost
4444  cells and that we don't own based on the tie-breaking
4445  criterion. unify dofs afterwards.
4446  3/ unify dofs and re-enumerate the remaining valid ones.
4447  the end result is that we only enumerate locally owned
4448  DoFs
4449  4/ shift indices so that each processor has a unique
4450  range of indices
4451  5/ for all locally owned cells that are ghost
4452  cells somewhere else, send our own DoF indices
4453  to the appropriate set of other processors.
4454  overwrite invalid DoF indices on ghost interfaces
4455  with the corresponding valid ones that we now know.
4456  6/ send DoF indices again to get the correct indices
4457  on ghost cells that we may not have known earlier
4458  */
4459 
4460  // --------- Phase 1: enumerate dofs on locally owned cells
4461  const types::global_dof_index n_initial_local_dofs =
4462  Implementation::distribute_dofs(subdomain_id, *dof_handler);
4463 
4464  // --------- Phase 2: eliminate dof duplicates on all cells:
4465  // - un-numerate dofs on interfaces to ghost cells
4466  // that we don't own
4467  // - in case of hp::DoFHandler, unify dofs
4468  std::vector<::types::global_dof_index> renumbering(
4469  n_initial_local_dofs, enumeration_dof_index);
4470 
4471  // first, we invalidate degrees of freedom that belong to processors
4472  // of a lower rank, from which we will receive the final (and lower)
4473  // degrees of freedom later.
4474  Implementation::
4475  invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering(
4476  renumbering, subdomain_id, *dof_handler);
4477 
4478  // then, we identify DoF duplicates if a hp::DoFHandler is used
4479  std::vector<std::map<types::global_dof_index, types::global_dof_index>>
4480  all_constrained_indices(dim);
4481  Implementation::compute_dof_identities(all_constrained_indices,
4482  *dof_handler);
4483 
4484  // --------- Phase 3: re-enumerate the valid degrees of freedom
4485  // consecutively. thus, we finally receive the
4486  // correct number of locally owned DoFs after
4487  // this step.
4488  //
4489  // the order in which we handle Phases 2 and 3 is important,
4490  // since we want to clarify ownership of degrees of freedom before
4491  // we actually unify and enumerate their indices. otherwise, we could
4492  // end up having a degee of freedom to which only invalid indices will
4493  // be assigned.
4494  const types::global_dof_index n_locally_owned_dofs =
4495  Implementation::enumerate_dof_indices_for_renumbering(
4496  renumbering, all_constrained_indices, *dof_handler);
4497 
4498  // --------- Phase 4: shift indices so that each processor has a unique
4499  // range of indices
4500  ::types::global_dof_index my_shift = 0;
4501  const int ierr =
4502  MPI_Exscan(DEAL_II_MPI_CONST_CAST(&n_locally_owned_dofs),
4503  &my_shift,
4504  1,
4505  DEAL_II_DOF_INDEX_MPI_TYPE,
4506  MPI_SUM,
4507  triangulation->get_communicator());
4508  AssertThrowMPI(ierr);
4509 
4510  // make dof indices globally consecutive
4511  for (auto &new_index : renumbering)
4512  if (new_index != numbers::invalid_dof_index)
4513  new_index += my_shift;
4514 
4515  // now re-enumerate all dofs to this shifted and condensed
4516  // numbering form. we renumber some dofs as invalid, so
4517  // choose the nocheck-version.
4518  Implementation::renumber_dofs(renumbering,
4519  IndexSet(0),
4520  *dof_handler,
4521  /*check_validity=*/false);
4522 
4523  // now a little bit of housekeeping
4524  const ::types::global_dof_index n_global_dofs =
4525  Utilities::MPI::sum(n_locally_owned_dofs,
4526  triangulation->get_communicator());
4527 
4528  NumberCache number_cache;
4529  number_cache.n_global_dofs = n_global_dofs;
4530  number_cache.n_locally_owned_dofs = n_locally_owned_dofs;
4531  number_cache.locally_owned_dofs = IndexSet(n_global_dofs);
4532  number_cache.locally_owned_dofs.add_range(my_shift,
4533  my_shift +
4534  n_locally_owned_dofs);
4535  number_cache.locally_owned_dofs.compress();
4536 
4537  // this ends the phase where we enumerate degrees of freedom on
4538  // each processor. what is missing is communicating DoF indices
4539  // on ghost cells
4540 
4541  // --------- Phase 5: for all locally owned cells that are ghost
4542  // cells somewhere else, send our own DoF indices
4543  // to the appropriate set of other processors
4544  {
4545  std::vector<bool> user_flags;
4546  triangulation->save_user_flags(user_flags);
4547  triangulation->clear_user_flags();
4548 
4549  // figure out which cells are ghost cells on which we have
4550  // to exchange DoF indices
4551  const std::map<unsigned int, std::set<::types::subdomain_id>>
4552  vertices_with_ghost_neighbors =
4553  GridTools::compute_vertices_with_ghost_neighbors(*triangulation);
4554 
4555  // mark all cells that either have to send data (locally
4556  // owned cells that are adjacent to ghost neighbors in some
4557  // way) or receive data (all ghost cells) via the user flags
4558  for (const auto &cell : dof_handler->active_cell_iterators())
4559  if (cell->is_locally_owned())
4560  {
4561  for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
4562  if (vertices_with_ghost_neighbors.find(cell->vertex_index(
4563  v)) != vertices_with_ghost_neighbors.end())
4564  {
4565  cell->set_user_flag();
4566  break;
4567  }
4568  }
4569  else if (cell->is_ghost())
4570  cell->set_user_flag();
4571 
4572 
4573 
4574  // Send and receive cells. After this, only the local cells
4575  // are marked, that received new data. This has to be
4576  // communicated in a second communication step.
4577  //
4578  // as explained in the 'distributed' paper, this has to be
4579  // done twice
4580  communicate_dof_indices_on_marked_cells(
4581  *dof_handler, vertices_with_ghost_neighbors);
4582 
4583  // in case of hp::DoFHandlers, we may have received valid
4584  // indices of degrees of freedom that are dominated by a fe
4585  // object adjacent to a ghost interface.
4586  // thus, we overwrite the remaining invalid indices with
4587  // the valid ones in this step.
4588  Implementation::merge_invalid_dof_indices_on_ghost_interfaces(
4589  *dof_handler);
4590 
4591  // --------- Phase 6: all locally owned cells have their correct
4592  // DoF indices set. however, some ghost cells
4593  // may still have invalid ones. thus, exchange
4594  // one more time.
4595  communicate_dof_indices_on_marked_cells(
4596  *dof_handler, vertices_with_ghost_neighbors);
4597 
4598  // at this point, we must have taken care of the data transfer
4599  // on all cells we had previously marked. verify this
4600 # ifdef DEBUG
4601  for (const auto &cell : dof_handler->active_cell_iterators())
4602  Assert(cell->user_flag_set() == false, ExcInternalError());
4603 # endif
4604 
4605  triangulation->load_user_flags(user_flags);
4606  }
4607 
4608 # ifdef DEBUG
4609  // check that we are really done
4610  {
4611  std::vector<::types::global_dof_index> local_dof_indices;
4612 
4613  for (const auto &cell : dof_handler->active_cell_iterators())
4614  if (!cell->is_artificial())
4615  {
4616  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
4617  cell->get_dof_indices(local_dof_indices);
4618  if (local_dof_indices.end() !=
4619  std::find(local_dof_indices.begin(),
4620  local_dof_indices.end(),
4621  numbers::invalid_dof_index))
4622  {
4623  if (cell->is_ghost())
4624  {
4625  Assert(false,
4626  ExcMessage(
4627  "A ghost cell ended up with incomplete "
4628  "DoF index information. This should not "
4629  "have happened!"));
4630  }
4631  else
4632  {
4633  Assert(
4634  false,
4635  ExcMessage(
4636  "A locally owned cell ended up with incomplete "
4637  "DoF index information. This should not "
4638  "have happened!"));
4639  }
4640  }
4641  }
4642  }
4643 # endif // DEBUG
4644  return number_cache;
4645 #endif // DEAL_II_WITH_MPI
4646  }
4647 
4648 
4649 
4650  template <class DoFHandlerType>
4651  std::vector<NumberCache>
4652  ParallelDistributed<DoFHandlerType>::distribute_mg_dofs() const
4653  {
4654 #ifndef DEAL_II_WITH_MPI
4655  Assert(false, ExcNotImplemented());
4656  return std::vector<NumberCache>();
4657 #else
4658  const unsigned int dim = DoFHandlerType::dimension;
4659  const unsigned int spacedim = DoFHandlerType::space_dimension;
4660 
4661  ::parallel::DistributedTriangulationBase<dim, spacedim>
4662  *triangulation =
4663  (dynamic_cast<
4664  ::parallel::DistributedTriangulationBase<dim, spacedim> *>(
4665  const_cast<::Triangulation<dim, spacedim> *>(
4666  &dof_handler->get_triangulation())));
4667  Assert(triangulation != nullptr, ExcInternalError());
4668 
4669  AssertThrow((triangulation->is_multilevel_hierarchy_constructed()),
4670  ExcMessage(
4671  "Multigrid DoFs can only be distributed on a parallel "
4672  "Triangulation if the flag construct_multigrid_hierarchy "
4673  "is set in the constructor."));
4674 
4675  // loop over all levels that exist globally (across all
4676  // processors), even if the current processor does not in fact
4677  // have any cells on that level or if the local part of the
4678  // Triangulation has fewer levels. we need to do this because
4679  // we need to communicate across all processors on all levels
4680  const unsigned int n_levels = triangulation->n_global_levels();
4681  std::vector<NumberCache> number_caches;
4682  number_caches.reserve(n_levels);
4683  for (unsigned int level = 0; level < n_levels; ++level)
4684  {
4685  NumberCache level_number_cache;
4686 
4687  //* 1. distribute on own subdomain
4688  const unsigned int n_initial_local_dofs =
4689  Implementation::distribute_dofs_on_level(
4690  triangulation->locally_owned_subdomain(), *dof_handler, level);
4691 
4692  //* 2. iterate over ghostcells and kill dofs that are not
4693  // owned by us
4694  std::vector<::types::global_dof_index> renumbering(
4695  n_initial_local_dofs);
4696  for (::types::global_dof_index i = 0; i < renumbering.size();
4697  ++i)
4698  renumbering[i] = i;
4699 
4700  if (level < triangulation->n_levels())
4701  {
4702  std::vector<::types::global_dof_index> local_dof_indices;
4703 
4704  typename DoFHandlerType::level_cell_iterator
4705  cell = dof_handler->begin(level),
4706  endc = dof_handler->end(level);
4707 
4708  for (; cell != endc; ++cell)
4709  if (cell->level_subdomain_id() !=
4710  numbers::artificial_subdomain_id &&
4711  (cell->level_subdomain_id() <
4712  triangulation->locally_owned_subdomain()))
4713  {
4714  // we found a neighboring ghost cell whose
4715  // subdomain is "stronger" than our own
4716  // subdomain
4717 
4718  // delete all dofs that live there and that we
4719  // have previously assigned a number to
4720  // (i.e. the ones on the interface)
4721  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
4722  cell->get_mg_dof_indices(local_dof_indices);
4723  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
4724  ++i)
4725  if (local_dof_indices[i] != numbers::invalid_dof_index)
4726  renumbering[local_dof_indices[i]] =
4727  numbers::invalid_dof_index;
4728  }
4729  }
4730 
4731  level_number_cache.n_locally_owned_dofs = 0;
4732  for (types::global_dof_index &index : renumbering)
4733  if (index != numbers::invalid_dof_index)
4734  index = level_number_cache.n_locally_owned_dofs++;
4735 
4736  //* 3. communicate local dofcount and shift ids to make
4737  // them unique
4738  ::types::global_dof_index my_shift = 0;
4739  int ierr = MPI_Exscan(DEAL_II_MPI_CONST_CAST(
4740  &level_number_cache.n_locally_owned_dofs),
4741  &my_shift,
4742  1,
4743  DEAL_II_DOF_INDEX_MPI_TYPE,
4744  MPI_SUM,
4745  triangulation->get_communicator());
4746  AssertThrowMPI(ierr);
4747 
4748  // The last processor knows about the total number of dofs, so we
4749  // can use a cheaper broadcast rather than an MPI_Allreduce via
4750  // MPI::sum().
4751  level_number_cache.n_global_dofs =
4752  my_shift + level_number_cache.n_locally_owned_dofs;
4753  ierr = MPI_Bcast(&level_number_cache.n_global_dofs,
4754  1,
4755  DEAL_II_DOF_INDEX_MPI_TYPE,
4756  Utilities::MPI::n_mpi_processes(
4757  triangulation->get_communicator()) -
4758  1,
4759  triangulation->get_communicator());
4760 
4761  // shift indices
4762  for (types::global_dof_index &index : renumbering)
4763  if (index != numbers::invalid_dof_index)
4764  index += my_shift;
4765 
4766  // now re-enumerate all dofs to this shifted and condensed
4767  // numbering form. we renumber some dofs as invalid, so
4768  // choose the nocheck-version of the function
4769  //
4770  // of course there is nothing for us to renumber if the
4771  // level we are currently dealing with doesn't even exist
4772  // within the current triangulation, so skip renumbering
4773  // in that case
4774  if (level < triangulation->n_levels())
4775  Implementation::renumber_mg_dofs(
4776  renumbering, IndexSet(0), *dof_handler, level, false);
4777 
4778  // now a little bit of housekeeping
4779  level_number_cache.locally_owned_dofs =
4780  IndexSet(level_number_cache.n_global_dofs);
4781  level_number_cache.locally_owned_dofs.add_range(
4782  my_shift, my_shift + level_number_cache.n_locally_owned_dofs);
4783  level_number_cache.locally_owned_dofs.compress();
4784 
4785  number_caches.emplace_back(level_number_cache);
4786  }
4787 
4788 
4789  //* communicate ghost DoFs
4790  // We mark all ghost cells by setting the user_flag and then request
4791  // these cells from the corresponding owners. As this information
4792  // can be incomplete,
4793  {
4794  std::vector<bool> user_flags;
4795  triangulation->save_user_flags(user_flags);
4796  triangulation->clear_user_flags();
4797 
4798  // mark all ghost cells for transfer
4799  {
4800  typename DoFHandlerType::level_cell_iterator cell,
4801  endc = dof_handler->end();
4802  for (cell = dof_handler->begin(); cell != endc; ++cell)
4803  if (cell->level_subdomain_id() !=
4804  ::numbers::artificial_subdomain_id &&
4805  !cell->is_locally_owned_on_level())
4806  cell->set_user_flag();
4807  }
4808 
4809  // Phase 1. Request all marked cells from corresponding owners. If we
4810  // managed to get every DoF, remove the user_flag, otherwise we
4811  // will request them again in the step below.
4812  communicate_mg_ghost_cells(*triangulation, *dof_handler);
4813 
4814  // have a barrier so that sends from above and below this
4815  // place are not mixed up.
4816  //
4817  // this is necessary because above we just see if there are
4818  // messages and then receive them, without discriminating
4819  // where they come from and whether they were sent in phase
4820  // 1 or 2 in communicate_mg_ghost_cells() on another
4821  // processor. the need for a global communication step like
4822  // this barrier could be avoided by receiving messages
4823  // specifically from those processors from which we expect
4824  // messages, and by using different tags for phase 1 and 2,
4825  // but the cost of a barrier is negligible compared to
4826  // everything else we do here
4827  const int ierr = MPI_Barrier(triangulation->get_communicator());
4828  AssertThrowMPI(ierr);
4829 
4830  // Phase 2, only request the cells that were not completed
4831  // in Phase 1.
4832  communicate_mg_ghost_cells(*triangulation, *dof_handler);
4833 
4834 # ifdef DEBUG
4835  // make sure we have removed all flags:
4836  {
4837  typename DoFHandlerType::level_cell_iterator cell,
4838  endc = dof_handler->end();
4839  for (cell = dof_handler->begin(); cell != endc; ++cell)
4840  if (cell->level_subdomain_id() !=
4841  ::numbers::artificial_subdomain_id &&
4842  !cell->is_locally_owned_on_level())
4843  Assert(cell->user_flag_set() == false, ExcInternalError());
4844  }
4845 # endif
4846 
4847  triangulation->load_user_flags(user_flags);
4848  }
4849 
4850 
4851 
4852 # ifdef DEBUG
4853  // check that we are really done
4854  {
4855  std::vector<::types::global_dof_index> local_dof_indices;
4856  typename DoFHandlerType::level_cell_iterator cell,
4857  endc = dof_handler->end();
4858 
4859  for (cell = dof_handler->begin(); cell != endc; ++cell)
4860  if (cell->level_subdomain_id() !=
4861  ::numbers::artificial_subdomain_id)
4862  {
4863  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
4864  cell->get_mg_dof_indices(local_dof_indices);
4865  if (local_dof_indices.end() !=
4866  std::find(local_dof_indices.begin(),
4867  local_dof_indices.end(),
4868  numbers::invalid_dof_index))
4869  {
4870  Assert(false, ExcMessage("not all DoFs got distributed!"));
4871  }
4872  }
4873  }
4874 # endif // DEBUG
4875 
4876  return number_caches;
4877 
4878 #endif // DEAL_II_WITH_MPI
4879  }
4880 
4881 
4882  template <class DoFHandlerType>
4883  NumberCache
4884  ParallelDistributed<DoFHandlerType>::renumber_dofs(
4885  const std::vector<::types::global_dof_index> &new_numbers) const
4886  {
4887  (void)new_numbers;
4888 
4889  Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
4890  ExcInternalError());
4891 
4892 #ifndef DEAL_II_WITH_MPI
4893  Assert(false, ExcNotImplemented());
4894  return NumberCache();
4895 #else
4896  const unsigned int dim = DoFHandlerType::dimension;
4897  const unsigned int spacedim = DoFHandlerType::space_dimension;
4898 
4899  ::parallel::DistributedTriangulationBase<dim, spacedim>
4900  *triangulation =
4901  (dynamic_cast<
4902  ::parallel::DistributedTriangulationBase<dim, spacedim> *>(
4903  const_cast<::Triangulation<dim, spacedim> *>(
4904  &dof_handler->get_triangulation())));
4905  Assert(triangulation != nullptr, ExcInternalError());
4906 
4907 
4908  // We start by checking whether only the numbering within the MPI
4909  // ranks changed. In that case, we can apply the renumbering with some
4910  // local renumbering only (this is similar to the renumber_mg_dofs()
4911  // function below)
4912  const bool locally_owned_set_changes =
4913  std::any_of(new_numbers.cbegin(),
4914  new_numbers.cend(),
4915  [this](const types::global_dof_index i) {
4916  return dof_handler->locally_owned_dofs().is_element(
4917  i) == false;
4918  });
4919 
4920  if (Utilities::MPI::sum(static_cast<unsigned int>(
4921  locally_owned_set_changes),
4922  triangulation->get_communicator()) == 0)
4923  {
4924  // Since only the order within the local subdomains has changed,
4925  // all we need to do is to propagate the knowledge about the
4926  // numbers from the locally owned dofs (given by the new_numbers
4927  // array) to all ghosted dofs on neighboring processors. We can do
4928  // this by ghost layer exchange routines as in parallel vectors:
4929  // We create an IndexSet for the relevant dofs and then export
4930  // into an array of those values via Utilities::MPI::Partitioner.
4931  IndexSet relevant_dofs;
4932  DoFTools::extract_locally_relevant_dofs(*dof_handler,
4933  relevant_dofs);
4934  std::vector<types::global_dof_index> ghosted_new_numbers(
4935  relevant_dofs.n_elements());
4936  {
4937  Utilities::MPI::Partitioner partitioner(
4938  dof_handler->locally_owned_dofs(),
4939  relevant_dofs,
4940  triangulation->get_communicator());
4941 
4942  // choose some number that makes it unlikely to get conflicts
4943  // with other ongoing non-blocking communication (there
4944  // shouldn't be any at this place in most programs).
4945  const unsigned int communication_channel = 19;
4946  std::vector<types::global_dof_index> temp_array(
4947  partitioner.n_import_indices());
4948  std::vector<MPI_Request> requests;
4949  partitioner.export_to_ghosted_array_start(
4950  communication_channel,
4951  make_array_view(new_numbers),
4952  make_array_view(temp_array),
4953  ArrayView<types::global_dof_index>(
4954  ghosted_new_numbers.data() + new_numbers.size(),
4955  partitioner.n_ghost_indices()),
4956  requests);
4957  partitioner.export_to_ghosted_array_finish(
4958  ArrayView<types::global_dof_index>(
4959  ghosted_new_numbers.data() + new_numbers.size(),
4960  partitioner.n_ghost_indices()),
4961  requests);
4962 
4963  // we need to fill the indices of the locally owned part into
4964  // the new numbers array, which is not provided by the parallel
4965  // partitioner. their right position is somewhere in the middle
4966  // of the array, so we first copy the ghosted part from smaller
4967  // ranks to the front, then insert the data in the middle.
4968  unsigned int n_ghosts_on_smaller_ranks = 0;
4969  for (std::pair<unsigned int, unsigned int> t :
4970  partitioner.ghost_targets())
4971  {
4972  if (t.first > partitioner.this_mpi_process())
4973  break;
4974  n_ghosts_on_smaller_ranks += t.second;
4975  }
4976  if (n_ghosts_on_smaller_ranks > 0)
4977  {
4978  Assert(ghosted_new_numbers.data() != nullptr,
4979  ExcInternalError());
4980  std::memmove(ghosted_new_numbers.data(),
4981  ghosted_new_numbers.data() + new_numbers.size(),
4982  sizeof(types::global_dof_index) *
4983  n_ghosts_on_smaller_ranks);
4984  }
4985  if (new_numbers.size() > 0)
4986  {
4987  Assert(new_numbers.data() != nullptr, ExcInternalError());
4988  std::memcpy(ghosted_new_numbers.data() +
4989  n_ghosts_on_smaller_ranks,
4990  new_numbers.data(),
4991  sizeof(types::global_dof_index) *
4992  new_numbers.size());
4993  }
4994  }
4995 
4996  // In case we do not carry any relevant dof (but only some remote
4997  // processor), we do not need to call the renumbering. We call the
4998  // version without validity check because vertex dofs will be
4999  // set already in the artificial region.
5000  if (relevant_dofs.n_elements() > 0)
5001  Implementation::renumber_dofs(ghosted_new_numbers,
5002  relevant_dofs,
5003  *dof_handler,
5004  /*check_validity=*/false);
5005 
5006  NumberCache number_cache;
5007  number_cache.locally_owned_dofs = dof_handler->locally_owned_dofs();
5008  number_cache.n_global_dofs = dof_handler->n_dofs();
5009  number_cache.n_locally_owned_dofs =
5010  number_cache.locally_owned_dofs.n_elements();
5011  return number_cache;
5012  }
5013  else
5014  {
5015  // Now back to the more complicated case
5016  //
5017  // First figure out the new set of locally owned DoF indices.
5018  // If we own no DoFs, we still need to go through this function,
5019  // but we can skip this calculation.
5020  //
5021  // The IndexSet::add_indices() function is substantially more
5022  // efficient if the set of indices is already sorted because
5023  // it can then insert ranges instead of individual elements.
5024  // consequently, pre-sort the array of new indices
5025  IndexSet my_locally_owned_new_dof_indices(dof_handler->n_dofs());
5026  if (dof_handler->n_locally_owned_dofs() > 0)
5027  {
5028  std::vector<::types::global_dof_index>
5029  new_numbers_sorted = new_numbers;
5030  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
5031 
5032  my_locally_owned_new_dof_indices.add_indices(
5033  new_numbers_sorted.begin(), new_numbers_sorted.end());
5034  my_locally_owned_new_dof_indices.compress();
5035 
5036  Assert(my_locally_owned_new_dof_indices.n_elements() ==
5037  new_numbers.size(),
5038  ExcInternalError());
5039  }
5040 
5041  // delete all knowledge of DoF indices that are not locally
5042  // owned. we do so by getting DoF indices on cells, checking
5043  // whether they are locally owned, if not, setting them to
5044  // an invalid value, and then setting them again on the current
5045  // cell
5046  //
5047  // DoFs we (i) know about, and (ii) don't own locally must be
5048  // located either on ghost cells, or on the interface between a
5049  // locally owned cell and a ghost cell. In any case, it is
5050  // sufficient to kill them only from the ghost side cell, so loop
5051  // only over ghost cells
5052  {
5053  std::vector<::types::global_dof_index> local_dof_indices;
5054 
5055  for (auto cell : dof_handler->active_cell_iterators())
5056  if (cell->is_ghost())
5057  {
5058  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
5059  cell->get_dof_indices(local_dof_indices);
5060 
5061  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
5062  ++i)
5063  // delete a DoF index if it has not already been deleted
5064  // (e.g., by visiting a neighboring cell, if it is on the
5065  // boundary), and if we don't own it
5066  if ((local_dof_indices[i] !=
5067  numbers::invalid_dof_index) &&
5068  (!dof_handler->locally_owned_dofs().is_element(
5069  local_dof_indices[i])))
5070  local_dof_indices[i] = numbers::invalid_dof_index;
5071 
5072  cell->set_dof_indices(local_dof_indices);
5073  }
5074  }
5075 
5076 
5077  // renumber. Skip when there is nothing to do because we own no DoF.
5078  if (dof_handler->locally_owned_dofs().n_elements() > 0)
5079  Implementation::renumber_dofs(new_numbers,
5080  dof_handler->locally_owned_dofs(),
5081  *dof_handler,
5082  /*check_validity=*/false);
5083 
5084  // Communicate newly assigned DoF indices to other processors
5085  // and get the same information for our own ghost cells.
5086  //
5087  // This is the same as phase 5+6 in the distribute_dofs() algorithm,
5088  // taking into account that we have to unify a few DoFs in between
5089  // then communication phases if we do hp numbering
5090  {
5091  std::vector<bool> user_flags;
5092  triangulation->save_user_flags(user_flags);
5093  triangulation->clear_user_flags();
5094 
5095  // mark all own cells for transfer
5096  for (const auto &cell : dof_handler->active_cell_iterators())
5097  if (!cell->is_artificial())
5098  cell->set_user_flag();
5099 
5100  // figure out which cells are ghost cells on which we have
5101  // to exchange DoF indices
5102  const std::map<unsigned int,
5103  std::set<::types::subdomain_id>>
5104  vertices_with_ghost_neighbors =
5105  GridTools::compute_vertices_with_ghost_neighbors(
5106  *triangulation);
5107 
5108 
5109  // Send and receive cells. After this, only the local cells
5110  // are marked, that received new data. This has to be
5111  // communicated in a second communication step.
5112  //
5113  // as explained in the 'distributed' paper, this has to be
5114  // done twice
5115  communicate_dof_indices_on_marked_cells(
5116  *dof_handler, vertices_with_ghost_neighbors);
5117 
5118  // in case of hp::DoFHandlers, we may have received valid
5119  // indices of degrees of freedom that are dominated by a fe
5120  // object adjacent to a ghost interface.
5121  // thus, we overwrite the remaining invalid indices with
5122  // the valid ones in this step.
5123  Implementation::merge_invalid_dof_indices_on_ghost_interfaces(
5124  *dof_handler);
5125 
5126  communicate_dof_indices_on_marked_cells(
5127  *dof_handler, vertices_with_ghost_neighbors);
5128 
5129  triangulation->load_user_flags(user_flags);
5130  }
5131 
5132  NumberCache number_cache;
5133  number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
5134  number_cache.n_global_dofs = dof_handler->n_dofs();
5135  number_cache.n_locally_owned_dofs =
5136  number_cache.locally_owned_dofs.n_elements();
5137  return number_cache;
5138  }
5139 #endif
5140  }
5141 
5142 
5143 
5144  template <class DoFHandlerType>
5145  NumberCache
5146  ParallelDistributed<DoFHandlerType>::renumber_mg_dofs(
5147  const unsigned int level,
5148  const std::vector<types::global_dof_index> &new_numbers) const
5149  {
5150  // we only implement the case where the multigrid numbers are
5151  // renumbered within the processor's partition, rather than the most
5152  // general case
5153  const IndexSet index_set = dof_handler->locally_owned_mg_dofs(level);
5154 
5155 #ifdef DEAL_II_WITH_MPI
5156 
5157  constexpr int dim = DoFHandlerType::dimension;
5158  constexpr int spacedim = DoFHandlerType::space_dimension;
5159  const ::parallel::TriangulationBase<dim, spacedim> *tr =
5160  (dynamic_cast<const ::parallel::TriangulationBase<dim, spacedim>
5161  *>(&this->dof_handler->get_triangulation()));
5162  Assert(tr != nullptr, ExcInternalError());
5163 
5164  const unsigned int my_rank =
5165  Utilities::MPI::this_mpi_process(tr->get_communicator());
5166 
5167 # ifdef DEBUG
5168  for (types::global_dof_index i : new_numbers)
5169  {
5170  Assert(index_set.is_element(i),
5171  ExcNotImplemented(
5172  "Renumberings that change the locally owned mg dofs "
5173  "partitioning are currently not implemented for "
5174  "the multigrid levels"));
5175  }
5176 # endif
5177 
5178  // we need to access all locally relevant degrees of freedom. we
5179  // use Utilities::MPI::Partitioner for handling the data exchange
5180  // of the new numbers, which is simply the extraction of ghost data
5181  IndexSet relevant_dofs;
5182  DoFTools::extract_locally_relevant_level_dofs(*dof_handler,
5183  level,
5184  relevant_dofs);
5185  std::vector<types::global_dof_index> ghosted_new_numbers(
5186  relevant_dofs.n_elements());
5187  {
5188  Utilities::MPI::Partitioner partitioner(index_set,
5189  relevant_dofs,
5190  tr->get_communicator());
5191  std::vector<types::global_dof_index> temp_array(
5192  partitioner.n_import_indices());
5193  const unsigned int communication_channel = 17;
5194  std::vector<MPI_Request> requests;
5195  partitioner.export_to_ghosted_array_start(
5196  communication_channel,
5197  make_array_view(new_numbers),
5198  make_array_view(temp_array),
5199  ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
5200  new_numbers.size(),
5201  partitioner.n_ghost_indices()),
5202  requests);
5203  partitioner.export_to_ghosted_array_finish(
5204  ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
5205  new_numbers.size(),
5206  partitioner.n_ghost_indices()),
5207  requests);
5208 
5209  // we need to fill the indices of the locally owned part into the
5210  // new numbers array. their right position is somewhere in the
5211  // middle of the array, so we first copy the ghosted part from
5212  // smaller ranks to the front, then insert the data in the middle.
5213  unsigned int n_ghosts_on_smaller_ranks = 0;
5214  for (std::pair<unsigned int, unsigned int> t :
5215  partitioner.ghost_targets())
5216  {
5217  if (t.first > my_rank)
5218  break;
5219  n_ghosts_on_smaller_ranks += t.second;
5220  }
5221  if (n_ghosts_on_smaller_ranks > 0)
5222  {
5223  Assert(ghosted_new_numbers.data() != nullptr, ExcInternalError());
5224  std::memmove(ghosted_new_numbers.data(),
5225  ghosted_new_numbers.data() + new_numbers.size(),
5226  sizeof(types::global_dof_index) *
5227  n_ghosts_on_smaller_ranks);
5228  }
5229  if (new_numbers.size() > 0)
5230  {
5231  Assert(new_numbers.data() != nullptr, ExcInternalError());
5232  std::memcpy(ghosted_new_numbers.data() +
5233  n_ghosts_on_smaller_ranks,
5234  new_numbers.data(),
5235  sizeof(types::global_dof_index) * new_numbers.size());
5236  }
5237  }
5238 
5239  // in case we do not own any of the given level (but only some remote
5240  // processor), we do not need to call the renumbering
5241  if (level < this->dof_handler->get_triangulation().n_levels() &&
5242  relevant_dofs.n_elements() > 0)
5243  Implementation::renumber_mg_dofs(
5244  ghosted_new_numbers, relevant_dofs, *dof_handler, level, true);
5245 #else
5246  (void)new_numbers;
5247  Assert(false, ExcNotImplemented());
5248 #endif
5249 
5250  NumberCache number_cache;
5251  number_cache.locally_owned_dofs = index_set;
5252  number_cache.n_global_dofs = dof_handler->n_dofs();
5253  number_cache.n_locally_owned_dofs =
5254  number_cache.locally_owned_dofs.n_elements();
5255  return number_cache;
5256  }
5257  } // namespace Policy
5258  } // namespace DoFHandlerImplementation
5259 } // namespace internal
5260 
5261 
5262 
5263 /*-------------- Explicit Instantiations -------------------------------*/
5264 #include "dof_handler_policy.inst"
5265 
5266 
5267 DEAL_II_NAMESPACE_CLOSE
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_mg_dofs
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)
Definition: dof_handler_policy.cc:2870
IndexSet::compress
void compress() const
Definition: index_set.h:1643
IndexSet
Definition: index_set.h:74
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
static void renumber_face_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, hp::DoFHandler< 1, spacedim > &)
Definition: dof_handler_policy.cc:2350
FiniteElement::hp_line_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:922
internal::DoFHandlerImplementation::Policy::Implementation::unify_dof_indices
static types::global_dof_index unify_dof_indices(const DoFHandler< dim, spacedim > &, const unsigned int n_dofs_before_identification, const bool)
Definition: dof_handler_policy.cc:1010
internal::DoFHandlerImplementation::Policy::ParallelDistributed::distribute_dofs
virtual NumberCache distribute_dofs() const override
Definition: dof_handler_policy.cc:4415
numbers::invalid_dof_index
const types::global_dof_index invalid_dof_index
Definition: types.h:206
internal::DoFHandlerImplementation::Policy::ParallelDistributed::ParallelDistributed
ParallelDistributed(DoFHandlerType &dof_handler)
Definition: dof_handler_policy.cc:4406
internal::DoFHandlerImplementation::Policy::Implementation::renumber_dofs
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)
Definition: dof_handler_policy.cc:2651
DoFHandler::faces
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > faces
Definition: dof_handler.h:1349
DoFHandler::cell_iterator
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:357
internal::DoFHandlerImplementation::Policy::ParallelShared::renumber_dofs
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
Definition: dof_handler_policy.cc:3641
parallel::DistributedTriangulationBase
Definition: tria_base.h:352
DoFHandler::levels
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > levels
Definition: dof_handler.h:1337
Utilities::MPI::Partitioner
Definition: partitioner.h:195
internal::DoFHandlerImplementation::Policy::ParallelShared::distribute_dofs
virtual NumberCache distribute_dofs() const override
Definition: dof_handler_policy.cc:3287
MultithreadInfo::n_threads
static unsigned int n_threads()
Definition: multithread_info.cc:169
DoFHandler::level_cell_iterator
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:393
DEAL_II_DOF_INDEX_MPI_TYPE
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition: types.h:86
hp::DoFHandler::set_dof_index
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
Definition: dof_handler.cc:2348
internal::DoFHandlerImplementation::Policy::Implementation::invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering
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)
Definition: dof_handler_policy.cc:1706
FiniteElement::hp_quad_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:933
internal::DoFHandlerImplementation::Policy::Sequential::distribute_dofs
virtual NumberCache distribute_dofs() const override
Definition: dof_handler_policy.cc:3049
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Triangulation
Definition: tria.h:1109
tria.h
internal::DoFHandlerImplementation::Policy::ParallelDistributed::renumber_mg_dofs
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
Definition: dof_handler_policy.cc:5146
hp::DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:312
hp::DoFHandler::get_fe_collection
const hp::FECollection< dim, spacedim > & get_fe_collection() const
ArrayView
Definition: array_view.h:77
memory_consumption.h
IndexSet::size
size_type size() const
Definition: index_set.h:1635
utilities.h
internal::DoFHandlerImplementation::Policy::Implementation::renumber_cell_dofs
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)
Definition: dof_handler_policy.cc:2291
hp::DoFHandler::line_iterator
typename ActiveSelector::line_iterator line_iterator
Definition: dof_handler.h:246
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_quad_dofs_on_ghost_interfaces
static void merge_invalid_quad_dofs_on_ghost_interfaces(hp::DoFHandler< dim, spacedim > &)
Definition: dof_handler_policy.cc:1462
tria_iterator.h
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_dof_indices_on_ghost_interfaces
static void merge_invalid_dof_indices_on_ghost_interfaces(hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler_policy.cc:1613
hp::DoFHandler::fe_collection
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1051
internal::DoFHandlerImplementation::Policy::ParallelDistributed::distribute_mg_dofs
virtual std::vector< NumberCache > distribute_mg_dofs() const override
Definition: dof_handler_policy.cc:4652
GeometryInfo
Definition: geometry_info.h:1224
DoFTools::max_dofs_per_cell
unsigned int max_dofs_per_cell(const DoFHandler< dim, spacedim > &dh)
partitioner.h
hp::DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
Utilities::MPI::internal::Tags::dofhandler_communicate_mg_ghost_cells
@ dofhandler_communicate_mg_ghost_cells
dof_handler_policy.cc: communicate_mg_ghost_cells()
Definition: mpi_tags.h:66
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition: tria_iterator_base.h:38
internal::DoFHandlerImplementation::Policy::Sequential::Sequential
Sequential(DoFHandlerType &dof_handler)
Definition: dof_handler_policy.cc:3041
make_array_view
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
hp::DoFHandler::quad_iterator
typename ActiveSelector::quad_iterator quad_iterator
Definition: dof_handler.h:270
internal::DoFHandlerImplementation::Policy::Implementation::enumerate_dof_indices_for_renumbering
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 &)
Definition: dof_handler_policy.cc:944
internal::DoFHandlerImplementation::Policy::Implementation::compute_dof_identities
static void compute_dof_identities(const std::vector< std::map< types::global_dof_index, types::global_dof_index >> &, const DoFHandler< dim, spacedim > &)
Definition: dof_handler_policy.cc:876
identity
Definition: template_constraints.h:268
Utilities::pack
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1209
second
Point< 2 > second
Definition: grid_out.cc:4353
internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs_on_cell
static types::global_dof_index distribute_mg_dofs_on_cell(const typename DoFHandler< dim, spacedim >::level_cell_iterator &cell, types::global_dof_index next_free_dof, const std::integral_constant< int, 1 > &)
Definition: dof_handler_policy.cc:1744
hp::DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index) const
DoFTools::locally_owned_dofs_per_subdomain
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandlerType &dof_handler)
Definition: dof_tools.cc:1375
Threads::new_task
Task< RT > new_task(const std::function< RT()> &function)
Definition: thread_management.h:1647
hp::DoFHandler
Definition: dof_handler.h:203
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_mg_dofs
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)
Definition: dof_handler_policy.cc:2801
DoFHandler
Definition: dof_handler.h:205
internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs_on_cell
static types::global_dof_index distribute_mg_dofs_on_cell(const typename DoFHandler< dim, spacedim >::level_cell_iterator &cell, types::global_dof_index next_free_dof, const std::integral_constant< int, 2 > &)
Definition: dof_handler_policy.cc:1816
Utilities::MPI::internal::Tags::dofhandler_communicate_mg_ghost_cells_reply
@ dofhandler_communicate_mg_ghost_cells_reply
dof_handler_policy.cc: communicate_mg_ghost_cells()
Definition: mpi_tags.h:69
Table
Definition: table.h:699
internal::DoFHandlerImplementation::NumberCache::locally_owned_dofs
IndexSet locally_owned_dofs
Definition: number_cache.h:158
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_line_dofs_on_ghost_interfaces
static void merge_invalid_line_dofs_on_ghost_interfaces(hp::DoFHandler< 1, spacedim > &)
Definition: dof_handler_policy.cc:1193
IndexSet::n_elements
size_type n_elements() const
Definition: index_set.h:1833
internal::DoFHandlerImplementation::Policy::Implementation::compute_quad_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const hp::DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler_policy.cc:743
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
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)
Definition: dof_handler_policy.cc:2147
WorkStream::run
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:1185
thread_management.h
level
unsigned int level
Definition: grid_out.cc:4355
work_stream.h
internal::DoFHandlerImplementation::Policy::Implementation::distribute_dofs
static types::global_dof_index distribute_dofs(const types::subdomain_id subdomain_id, DoFHandlerType &dof_handler)
Definition: dof_handler_policy.cc:1653
hp::DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: dof_handler.cc:1379
internal::DoFHandlerImplementation::Policy::Implementation::renumber_cell_dofs
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)
Definition: dof_handler_policy.cc:2086
Utilities::MPI::Partitioner::export_to_ghosted_array_finish
void export_to_ghosted_array_finish(const ArrayView< Number, MemorySpaceType > &ghost_array, std::vector< MPI_Request > &requests) const
DoFHandler::begin
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:922
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
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)
Definition: dof_handler_policy.cc:2127
FiniteElement::hp_vertex_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:911
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
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)
Definition: dof_handler_policy.cc:2462
internal::DoFHandlerImplementation::Policy::Implementation::renumber_vertex_dofs
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)
Definition: dof_handler_policy.cc:2177
internal::DoFHandlerImplementation::NumberCache
Definition: number_cache.h:37
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
tree_index
std::vector< unsigned int > tree_index
Definition: tria.cc:2237
internal::DoFHandlerImplementation::Policy::ParallelShared::distribute_mg_dofs
virtual std::vector< NumberCache > distribute_mg_dofs() const override
Definition: dof_handler_policy.cc:3450
types::global_dof_index
unsigned int global_dof_index
Definition: types.h:76
internal::DoFHandlerImplementation::Policy::Implementation::unify_dof_indices
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)
Definition: dof_handler_policy.cc:1020
hp::DoFHandler::get_active_fe_indices
void get_active_fe_indices(std::vector< unsigned int > &active_fe_indices) const
Definition: dof_handler.cc:1554
geometry_info.h
internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs_on_cell
static types::global_dof_index distribute_mg_dofs_on_cell(const hp::DoFHandler< 1, spacedim > &dof_handler, const typename hp::DoFHandler< 1, spacedim >::active_cell_iterator &cell, types::global_dof_index next_free_dof)
Definition: dof_handler_policy.cc:1943
internal::DoFHandlerImplementation::Policy::Implementation::compute_vertex_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_vertex_dof_identities(const hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler_policy.cc:228
Utilities::MPI::Partitioner::this_mpi_process
unsigned int this_mpi_process() const
DoFHandler::mg_vertex_dofs
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1329
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_mg_dofs
static void renumber_face_mg_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &, const unsigned int, const bool)
Definition: dof_handler_policy.cc:2787
CellId
Definition: cell_id.h:69
grid_tools.h
GridTools::compute_vertices_with_ghost_neighbors
std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
Definition: grid_tools.cc:5517
Utilities::MPI::Partitioner::n_ghost_indices
unsigned int n_ghost_indices() const
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
Utilities::unpack
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1353
internal::DoFHandlerImplementation::Policy::Implementation::compute_line_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const hp::DoFHandler< 1, spacedim > &)
Definition: dof_handler_policy.cc:370
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs_on_cell
static types::global_dof_index distribute_mg_dofs_on_cell(const hp::DoFHandler< 2, spacedim > &dof_handler, const typename hp::DoFHandler< 2, spacedim >::active_cell_iterator &cell, types::global_dof_index next_free_dof)
Definition: dof_handler_policy.cc:1959
DoFHandler::vertex_dofs
std::vector< types::global_dof_index > vertex_dofs
Definition: dof_handler.h:1323
fe.h
shared_tria.h
Utilities::MPI::CollectiveMutex
Definition: mpi.h:322
IndexSet::is_element
bool is_element(const size_type index) const
Definition: index_set.h:1766
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_line_dofs_on_ghost_interfaces
static void merge_invalid_line_dofs_on_ghost_interfaces(hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler_policy.cc:1200
internal::DoFHandlerImplementation::NumberCache::n_locally_owned_dofs
types::global_dof_index n_locally_owned_dofs
Definition: number_cache.h:147
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_dof_indices_on_ghost_interfaces
static void merge_invalid_dof_indices_on_ghost_interfaces(const DoFHandler< dim, spacedim > &)
Definition: dof_handler_policy.cc:1606
internal::DoFHandlerImplementation::Policy::Implementation::compute_quad_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const hp::DoFHandler< dim, spacedim > &)
Definition: dof_handler_policy.cc:730
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1707
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
static void renumber_face_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &)
Definition: dof_handler_policy.cc:2115
Utilities::MPI::Partitioner::n_import_indices
unsigned int n_import_indices() const
DoFHandler::MGVertexDoFs
Definition: dof_handler.h:1221
internal::DoFHandlerImplementation::Policy::Implementation::renumber_cell_mg_dofs
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)
Definition: dof_handler_policy.cc:2751
internal::DoFHandlerImplementation::Policy::Implementation::distribute_dofs_on_level
static types::global_dof_index distribute_dofs_on_level(const types::subdomain_id level_subdomain_id, DoFHandlerType &dof_handler, const unsigned int level)
Definition: dof_handler_policy.cc:1991
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_quad_dofs_on_ghost_interfaces
static void merge_invalid_quad_dofs_on_ghost_interfaces(hp::DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler_policy.cc:1473
FiniteElement
Definition: fe.h:648
Threads::TaskGroup::join_all
void join_all() const
Definition: thread_management.h:1833
internal::DoFHandlerImplementation::Policy::Implementation::renumber_mg_dofs
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)
Definition: dof_handler_policy.cc:2982
internal::DoFHandlerImplementation::Policy::Sequential::renumber_dofs
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
Definition: dof_handler_policy.cc:3105
DoFHandler::end
cell_iterator end() const
Definition: dof_handler.cc:951
types::subdomain_id
unsigned int subdomain_id
Definition: types.h:43
vertices_with_ghost_neighbors
std::map< unsigned int, std::set<::types::subdomain_id > > * vertices_with_ghost_neighbors
Definition: p4est_wrappers.cc:72
internal::DoFHandlerImplementation::Policy::ParallelDistributed::renumber_dofs
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
Definition: dof_handler_policy.cc:4884
unsigned int
internal::DoFHandlerImplementation::Policy::Sequential::renumber_mg_dofs
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
Definition: dof_handler_policy.cc:3129
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
tria.h
internal::DoFHandlerImplementation::Policy::ParallelShared::ParallelShared
ParallelShared(DoFHandlerType &dof_handler)
Definition: dof_handler_policy.cc:3145
internal::DoFHandlerImplementation::NumberCache::n_global_dofs
types::global_dof_index n_global_dofs
Definition: number_cache.h:137
types::coarse_cell_id
global_cell_index coarse_cell_id
Definition: types.h:114
internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs_on_cell
static types::global_dof_index distribute_mg_dofs_on_cell(const hp::DoFHandler< 3, spacedim > &dof_handler, const typename hp::DoFHandler< 3, spacedim >::active_cell_iterator &cell, types::global_dof_index next_free_dof)
Definition: dof_handler_policy.cc:1975
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
IndexSet::nth_index_in_set
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1881
DEAL_II_MPI_CONST_CAST
#define DEAL_II_MPI_CONST_CAST(expr)
Definition: mpi.h:76
DoFTools::extract_locally_relevant_dofs
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1173
internal::DoFHandlerImplementation::Policy::Implementation
Definition: dof_handler_policy.cc:218
Utilities::MPI::Partitioner::ghost_targets
const std::vector< std::pair< unsigned int, unsigned int > > & ghost_targets() const
GridTools::exchange_cell_data_to_ghosts
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)
dof_handler.h
dof_accessor.h
Threads::TaskGroup
Definition: thread_management.h:1798
GridTools::shift
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:817
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
internal::DoFHandlerImplementation::Policy::Implementation::renumber_vertex_mg_dofs
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)
Definition: dof_handler_policy.cc:2695
hp
Definition: hp.h:117
Utilities::MPI::Partitioner::export_to_ghosted_array_start
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
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:191
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1922
CellId::binary_type
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:79
internal::DoFHandlerImplementation::Policy::Implementation::compute_dof_identities
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)
Definition: dof_handler_policy.cc:885
dof_handler_policy.h
internal::DoFHandlerImplementation::Policy::Implementation::compute_line_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler_policy.cc:378
DoFTools::extract_locally_relevant_level_dofs
void extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1215
numbers::invalid_subdomain_id
const types::subdomain_id invalid_subdomain_id
Definition: types.h:285
numbers::artificial_subdomain_id
const types::subdomain_id artificial_subdomain_id
Definition: types.h:302
internal::DoFHandlerImplementation::Policy::ParallelShared::renumber_mg_dofs
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
Definition: dof_handler_policy.cc:3810
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
internal
Definition: aligned_vector.h:369
internal::DoFHandlerImplementation::Policy::Implementation::renumber_vertex_dofs
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)
Definition: dof_handler_policy.cc:2046
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12618
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_vertex_dofs_on_ghost_interfaces
static void merge_invalid_vertex_dofs_on_ghost_interfaces(hp::DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler_policy.cc:1051
IndexSet::add_indices
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1704
DoFHandler::quad_iterator
typename ActiveSelector::quad_iterator quad_iterator
Definition: dof_handler.h:272
internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs_on_cell
static types::global_dof_index distribute_mg_dofs_on_cell(const typename DoFHandler< dim, spacedim >::level_cell_iterator &cell, types::global_dof_index next_free_dof, const std::integral_constant< int, 3 > &)
Definition: dof_handler_policy.cc:1870
first
Point< 2 > first
Definition: grid_out.cc:4352
CellId::to_cell
Triangulation< dim, spacedim >::cell_iterator to_cell(const Triangulation< dim, spacedim > &tria) const
Definition: cell_id.cc:159
internal::DoFHandlerImplementation::Policy::Implementation::renumber_mg_dofs
static void renumber_mg_dofs(const std::vector<::types::global_dof_index > &, const IndexSet &, hp::DoFHandler< dim, spacedim > &, const unsigned int, const bool)
Definition: dof_handler_policy.cc:3023
DoFHandler::line_iterator
typename ActiveSelector::line_iterator line_iterator
Definition: dof_handler.h:248
DoFHandler::mg_levels
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition: dof_handler.h:1341
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
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)
Definition: dof_handler_policy.cc:2363
quadrants
std::vector< typename ::internal::p4est::types< dim >::quadrant > quadrants
Definition: tria.cc:2241
IndexSet::add_range
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1674
internal::DoFHandlerImplementation::Policy::Sequential::distribute_mg_dofs
virtual std::vector< NumberCache > distribute_mg_dofs() const override
Definition: dof_handler_policy.cc:3068
Utilities::MPI::CollectiveMutex::ScopedLock
Definition: mpi.h:330
int
Generated by   doxygen 1.8.17