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