Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
dof_handler_policy.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/geometry_info.h>
18 #include <deal.II/base/memory_consumption.h>
19 #include <deal.II/base/partitioner.h>
20 #include <deal.II/base/thread_management.h>
21 #include <deal.II/base/utilities.h>
22 #include <deal.II/base/work_stream.h>
23 
24 #include <deal.II/distributed/shared_tria.h>
25 #include <deal.II/distributed/tria.h>
26 
27 #include <deal.II/dofs/dof_accessor.h>
28 #include <deal.II/dofs/dof_handler.h>
29 #include <deal.II/dofs/dof_handler_policy.h>
30 
31 #include <deal.II/fe/fe.h>
32 
33 #include <deal.II/grid/grid_tools.h>
34 #include <deal.II/grid/tria.h>
35 #include <deal.II/grid/tria_iterator.h>
36 
37 #include <algorithm>
38 #include <memory>
39 #include <numeric>
40 #include <set>
41 
42 DEAL_II_NAMESPACE_OPEN
43 
44 
45 namespace internal
46 {
47  namespace DoFHandlerImplementation
48  {
49  namespace Policy
50  {
51  // use class ::DoFHandler instead
52  // of namespace internal::DoFHandler in
53  // the following
54  using ::DoFHandler;
55 
56  namespace hp
57  {
58  using ::hp::DoFHandler;
59  }
60 
61 
62  namespace
63  {
70  const types::global_dof_index enumeration_dof_index =
72 
77  template <class DoFHandlerType>
78  void
79  update_all_active_cell_dof_indices_caches(
80  const DoFHandlerType &dof_handler)
81  {
82  typename DoFHandlerType::active_cell_iterator
83  beginc = dof_handler.begin_active(),
84  endc = dof_handler.end();
85 
86  auto worker =
87  [](const typename DoFHandlerType::active_cell_iterator &cell,
88  void *,
89  void *) {
90  if (!cell->is_artificial())
91  cell->update_cell_dof_indices_cache();
92  };
93 
94  // parallelize filling all of the cell caches. by using
95  // WorkStream, we make sure that we only run through the
96  // range of iterators once, whereas a parallel_for loop
97  // for example has to split the range multiple times,
98  // which is expensive because cell iterators are not
99  // random access iterators with a cheap operator-
100  WorkStream::run(beginc,
101  endc,
102  worker,
103  /* copier */ std::function<void(void *)>(),
104  /* scratch_data */ nullptr,
105  /* copy_data */ nullptr,
107  /* chunk_size = */ 32);
108  }
109 
110 
115  template <class DoFHandlerType>
116  void
117  update_all_level_cell_dof_indices_caches(
118  const DoFHandlerType &dof_handler)
119  {
120  typename DoFHandlerType::level_cell_iterator beginc =
121  dof_handler.begin(),
122  endc = dof_handler.end();
123 
124  auto worker =
125  [](const typename DoFHandlerType::level_cell_iterator &cell,
126  void *,
127  void *) {
128  if (cell->has_children() || !cell->is_artificial())
129  cell->update_cell_dof_indices_cache();
130  };
131 
132  // parallelize filling all of the cell caches. by using
133  // WorkStream, we make sure that we only run through the
134  // range of iterators once, whereas a parallel_for loop
135  // for example has to split the range multiple times,
136  // which is expensive because cell iterators are not
137  // random access iterators with a cheap operator-
138  WorkStream::run(beginc,
139  endc,
140  worker,
141  /* copier */ std::function<void(void *)>(),
142  /* scratch_data */ nullptr,
143  /* copy_data */ nullptr,
145  /* chunk_size = */ 32);
146  }
147 
148 
149  using DoFIdentities =
150  std::vector<std::pair<unsigned int, unsigned int>>;
151 
152 
163  template <int structdim, int dim, int spacedim>
164  void
165  ensure_existence_of_dof_identities(
166  const FiniteElement<dim, spacedim> &fe1,
167  const FiniteElement<dim, spacedim> &fe2,
168  std::unique_ptr<DoFIdentities> & identities)
169  {
170  // see if we need to fill this entry, or whether it already
171  // exists
172  if (identities.get() == nullptr)
173  {
174  switch (structdim)
175  {
176  case 0:
177  {
178  identities = std_cxx14::make_unique<DoFIdentities>(
179  fe1.hp_vertex_dof_identities(fe2));
180  break;
181  }
182 
183  case 1:
184  {
185  identities = std_cxx14::make_unique<DoFIdentities>(
186  fe1.hp_line_dof_identities(fe2));
187  break;
188  }
189 
190  case 2:
191  {
192  identities = std_cxx14::make_unique<DoFIdentities>(
193  fe1.hp_quad_dof_identities(fe2));
194  break;
195  }
196 
197  default:
198  Assert(false, ExcNotImplemented());
199  }
200 
201  // double check whether the newly created entries make
202  // any sense at all
203  for (unsigned int i = 0; i < identities->size(); ++i)
204  {
205  Assert((*identities)[i].first <
206  fe1.template n_dofs_per_object<structdim>(),
207  ExcInternalError());
208  Assert((*identities)[i].second <
209  fe2.template n_dofs_per_object<structdim>(),
210  ExcInternalError());
211  }
212  }
213  }
214  } // namespace
215 
216 
217 
218  struct Implementation
219  {
220  /* -------------- distribute_dofs functionality ------------- */
221 
229  template <int spacedim>
231  distribute_dofs_on_cell(
232  const DoFHandler<1, spacedim> &dof_handler,
234  types::global_dof_index next_free_dof)
235  {
236  // distribute dofs of vertices
237  if (dof_handler.get_fe().dofs_per_vertex > 0)
238  for (unsigned int v = 0; v < GeometryInfo<1>::vertices_per_cell;
239  ++v)
240  {
241  if (cell->vertex_dof_index(v, 0) == numbers::invalid_dof_index)
242  for (unsigned int d = 0;
243  d < dof_handler.get_fe().dofs_per_vertex;
244  ++d)
245  {
246  Assert((cell->vertex_dof_index(v, d) ==
248  ExcInternalError());
249  cell->set_vertex_dof_index(v, d, next_free_dof++);
250  }
251  else
252  for (unsigned int d = 0;
253  d < dof_handler.get_fe().dofs_per_vertex;
254  ++d)
255  Assert((cell->vertex_dof_index(v, d) !=
257  ExcInternalError());
258  }
259 
260  // dofs of line
261  for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_line; ++d)
262  cell->set_dof_index(d, next_free_dof++);
263 
264  return next_free_dof;
265  }
266 
267 
268 
269  template <int spacedim>
271  distribute_dofs_on_cell(
272  const DoFHandler<2, spacedim> &dof_handler,
274  types::global_dof_index next_free_dof)
275  {
276  if (dof_handler.get_fe().dofs_per_vertex > 0)
277  // number dofs on vertices
278  for (unsigned int vertex = 0;
279  vertex < GeometryInfo<2>::vertices_per_cell;
280  ++vertex)
281  // check whether dofs for this vertex have been distributed
282  // (checking the first dof should be good enough)
283  if (cell->vertex_dof_index(vertex, 0) ==
285  for (unsigned int d = 0;
286  d < dof_handler.get_fe().dofs_per_vertex;
287  ++d)
288  cell->set_vertex_dof_index(vertex, d, next_free_dof++);
289 
290  // for the four sides
291  if (dof_handler.get_fe().dofs_per_line > 0)
292  for (unsigned int side = 0; side < GeometryInfo<2>::faces_per_cell;
293  ++side)
294  {
295  const typename DoFHandler<2, spacedim>::line_iterator line =
296  cell->line(side);
297 
298  // distribute dofs if necessary: check whether line dof is
299  // already numbered (checking the first dof should be good
300  // enough)
301  if (line->dof_index(0) == numbers::invalid_dof_index)
302  // if not: distribute dofs
303  for (unsigned int d = 0;
304  d < dof_handler.get_fe().dofs_per_line;
305  ++d)
306  line->set_dof_index(d, next_free_dof++);
307  }
308 
309 
310  // dofs of quad
311  if (dof_handler.get_fe().dofs_per_quad > 0)
312  for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_quad;
313  ++d)
314  cell->set_dof_index(d, next_free_dof++);
315 
316  return next_free_dof;
317  }
318 
319 
320 
321  template <int spacedim>
323  distribute_dofs_on_cell(
324  const DoFHandler<3, spacedim> &dof_handler,
326  types::global_dof_index next_free_dof)
327  {
328  if (dof_handler.get_fe().dofs_per_vertex > 0)
329  // number dofs on vertices
330  for (unsigned int vertex = 0;
331  vertex < GeometryInfo<3>::vertices_per_cell;
332  ++vertex)
333  // check whether dofs for this vertex have been distributed
334  // (checking the first dof should be good enough)
335  if (cell->vertex_dof_index(vertex, 0) ==
337  for (unsigned int d = 0;
338  d < dof_handler.get_fe().dofs_per_vertex;
339  ++d)
340  cell->set_vertex_dof_index(vertex, d, next_free_dof++);
341 
342  // for the lines
343  if (dof_handler.get_fe().dofs_per_line > 0)
344  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
345  {
346  const typename DoFHandler<3, spacedim>::line_iterator line =
347  cell->line(l);
348 
349  // distribute dofs if necessary: check whether line dof is
350  // already numbered (checking the first dof should be good
351  // enough)
352  if (line->dof_index(0) == numbers::invalid_dof_index)
353  // if not: distribute dofs
354  for (unsigned int d = 0;
355  d < dof_handler.get_fe().dofs_per_line;
356  ++d)
357  line->set_dof_index(d, next_free_dof++);
358  }
359 
360  // for the quads
361  if (dof_handler.get_fe().dofs_per_quad > 0)
362  for (unsigned int q = 0; q < GeometryInfo<3>::quads_per_cell; ++q)
363  {
364  const typename DoFHandler<3, spacedim>::quad_iterator quad =
365  cell->quad(q);
366 
367  // distribute dofs if necessary: check whether line dof is
368  // already numbered (checking the first dof should be good
369  // enough)
370  if (quad->dof_index(0) == numbers::invalid_dof_index)
371  // if not: distribute dofs
372  for (unsigned int d = 0;
373  d < dof_handler.get_fe().dofs_per_quad;
374  ++d)
375  quad->set_dof_index(d, next_free_dof++);
376  }
377 
378 
379  // dofs of hex
380  if (dof_handler.get_fe().dofs_per_hex > 0)
381  for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_hex; ++d)
382  cell->set_dof_index(d, next_free_dof++);
383 
384  return next_free_dof;
385  }
386 
387 
388 
389  // same for the hp::DoFHandler
390  template <int spacedim>
392  distribute_dofs_on_cell(
395  & cell,
396  types::global_dof_index next_free_dof)
397  {
398  const unsigned int dim = 1;
399 
400  const FiniteElement<dim, spacedim> &fe = cell->get_fe();
401  const unsigned int fe_index = cell->active_fe_index();
402 
403  // number dofs on vertices. to do so, check whether dofs for
404  // this vertex have been distributed and for the present fe
405  // (only check the first dof), and if this isn't the case
406  // distribute new ones there
407  if (fe.dofs_per_vertex > 0)
408  for (unsigned int vertex = 0;
409  vertex < GeometryInfo<1>::vertices_per_cell;
410  ++vertex)
411  if (cell->vertex_dof_index(vertex, 0, fe_index) ==
413  for (unsigned int d = 0; d < fe.dofs_per_vertex;
414  ++d, ++next_free_dof)
415  cell->set_vertex_dof_index(vertex,
416  d,
417  next_free_dof,
418  fe_index);
419 
420  // finally for the line. this one shouldn't be numbered yet
421  if (fe.dofs_per_line > 0)
422  {
423  Assert((cell->dof_index(0, fe_index) ==
425  ExcInternalError());
426 
427  for (unsigned int d = 0; d < fe.dofs_per_line;
428  ++d, ++next_free_dof)
429  cell->set_dof_index(d, next_free_dof, fe_index);
430  }
431 
432  // note that this cell has been processed
433  cell->set_user_flag();
434 
435  return next_free_dof;
436  }
437 
438 
439 
440  template <int spacedim>
442  distribute_dofs_on_cell(
445  & cell,
446  types::global_dof_index next_free_dof)
447  {
448  const unsigned int dim = 2;
449 
450  const FiniteElement<dim, spacedim> &fe = cell->get_fe();
451  const unsigned int fe_index = cell->active_fe_index();
452 
453  // number dofs on vertices. to do so, check whether dofs for
454  // this vertex have been distributed and for the present fe
455  // (only check the first dof), and if this isn't the case
456  // distribute new ones there
457  if (fe.dofs_per_vertex > 0)
458  for (unsigned int vertex = 0;
459  vertex < GeometryInfo<2>::vertices_per_cell;
460  ++vertex)
461  if (cell->vertex_dof_index(vertex, 0, fe_index) ==
463  for (unsigned int d = 0; d < fe.dofs_per_vertex;
464  ++d, ++next_free_dof)
465  cell->set_vertex_dof_index(vertex,
466  d,
467  next_free_dof,
468  fe_index);
469 
470  // next the sides. do the same as above: check whether the
471  // line is already numbered for the present fe_index, and if
472  // not do it
473  if (fe.dofs_per_line > 0)
474  for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell; ++l)
475  {
477  line = cell->line(l);
478 
479  if (line->dof_index(0, fe_index) == numbers::invalid_dof_index)
480  for (unsigned int d = 0; d < fe.dofs_per_line;
481  ++d, ++next_free_dof)
482  line->set_dof_index(d, next_free_dof, fe_index);
483  }
484 
485 
486  // finally for the quad. this one shouldn't be numbered yet
487  if (fe.dofs_per_quad > 0)
488  {
489  Assert((cell->dof_index(0, fe_index) ==
491  ExcInternalError());
492 
493  for (unsigned int d = 0; d < fe.dofs_per_quad;
494  ++d, ++next_free_dof)
495  cell->set_dof_index(d, next_free_dof, fe_index);
496  }
497 
498  // note that this cell has been processed
499  cell->set_user_flag();
500 
501  return next_free_dof;
502  }
503 
504 
505 
506  template <int spacedim>
508  distribute_dofs_on_cell(
511  & cell,
512  types::global_dof_index next_free_dof)
513  {
514  const unsigned int dim = 3;
515 
516  const FiniteElement<dim, spacedim> &fe = cell->get_fe();
517  const unsigned int fe_index = cell->active_fe_index();
518 
519  // number dofs on vertices. to do so, check whether dofs for
520  // this vertex have been distributed and for the present fe
521  // (only check the first dof), and if this isn't the case
522  // distribute new ones there
523  if (fe.dofs_per_vertex > 0)
524  for (unsigned int vertex = 0;
525  vertex < GeometryInfo<3>::vertices_per_cell;
526  ++vertex)
527  if (cell->vertex_dof_index(vertex, 0, fe_index) ==
529  for (unsigned int d = 0; d < fe.dofs_per_vertex;
530  ++d, ++next_free_dof)
531  cell->set_vertex_dof_index(vertex,
532  d,
533  next_free_dof,
534  fe_index);
535 
536  // next the eight lines. do the same as above: check whether
537  // the line is already numbered for the present fe_index,
538  // and if not do it
539  if (fe.dofs_per_line > 0)
540  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
541  {
543  line = cell->line(l);
544 
545  if (line->dof_index(0, fe_index) == numbers::invalid_dof_index)
546  for (unsigned int d = 0; d < fe.dofs_per_line;
547  ++d, ++next_free_dof)
548  line->set_dof_index(d, next_free_dof, fe_index);
549  }
550 
551  // same for quads
552  if (fe.dofs_per_quad > 0)
553  for (unsigned int q = 0; q < GeometryInfo<3>::quads_per_cell; ++q)
554  {
556  quad = cell->quad(q);
557 
558  if (quad->dof_index(0, fe_index) == numbers::invalid_dof_index)
559  for (unsigned int d = 0; d < fe.dofs_per_quad;
560  ++d, ++next_free_dof)
561  quad->set_dof_index(d, next_free_dof, fe_index);
562  }
563 
564 
565  // finally for the hex. this one shouldn't be numbered yet
566  // because there is no other cell from which we could already
567  // have gotten to numbering anything on this hex (=cell)
568  if (fe.dofs_per_hex > 0)
569  {
570  Assert((cell->dof_index(0, fe_index) ==
572  ExcInternalError());
573 
574  for (unsigned int d = 0; d < fe.dofs_per_hex;
575  ++d, ++next_free_dof)
576  cell->set_dof_index(d, next_free_dof, fe_index);
577  }
578 
579  // note that this cell has been processed
580  cell->set_user_flag();
581 
582  return next_free_dof;
583  }
584 
585 
586 
591  template <int dim, int spacedim>
592  static std::map<types::global_dof_index, types::global_dof_index>
593  compute_vertex_dof_identities(
594  const hp::DoFHandler<dim, spacedim> &dof_handler)
595  {
596  std::map<types::global_dof_index, types::global_dof_index>
597  dof_identities;
598 
599  // Note: we may wish to have something here similar to what
600  // we do for lines and quads, namely that we only identify
601  // dofs for any fe towards the most dominating one. however,
602  // it is not clear whether this is actually necessary for
603  // vertices at all, I can't think of a finite element that
604  // would make that necessary...
606  vertex_dof_identities(dof_handler.get_fe_collection().size(),
607  dof_handler.get_fe_collection().size());
608 
609  // loop over all vertices and see which one we need to work on
610  for (unsigned int vertex_index = 0;
611  vertex_index < dof_handler.get_triangulation().n_vertices();
612  ++vertex_index)
613  if (dof_handler.get_triangulation()
614  .get_used_vertices()[vertex_index] == true)
615  {
616  const unsigned int n_active_fe_indices =
617  ::internal::DoFAccessorImplementation::Implementation::
618  n_active_vertex_fe_indices(dof_handler, vertex_index);
619 
620  if (n_active_fe_indices > 1)
621  {
622  const std::set<unsigned int> fe_indices =
623  ::internal::DoFAccessorImplementation::
624  Implementation::get_active_vertex_fe_indices(
625  dof_handler, vertex_index);
626 
627  // find out which is the most dominating finite
628  // element of the ones that are used on this vertex
629  unsigned int most_dominating_fe_index =
630  dof_handler.get_fe_collection().find_dominating_fe(
631  fe_indices,
632  /*codim*/ dim);
633 
634  // if we haven't found a dominating finite element,
635  // choose the very first one to be dominant
636  if (most_dominating_fe_index ==
638  most_dominating_fe_index =
639  ::internal::DoFAccessorImplementation::
640  Implementation::nth_active_vertex_fe_index(
641  dof_handler, vertex_index, 0);
642 
643  // loop over the indices of all the finite
644  // elements that are not dominating, and
645  // identify their dofs to the most dominating
646  // one
647  for (const auto &other_fe_index : fe_indices)
648  if (other_fe_index != most_dominating_fe_index)
649  {
650  // make sure the entry in the equivalence
651  // table exists
652  ensure_existence_of_dof_identities<0>(
653  dof_handler.get_fe(most_dominating_fe_index),
654  dof_handler.get_fe(other_fe_index),
655  vertex_dof_identities[most_dominating_fe_index]
656  [other_fe_index]);
657 
658  // then loop through the identities we
659  // have. first get the global numbers of the
660  // dofs we want to identify and make sure they
661  // are not yet constrained to anything else,
662  // except for to each other. use the rule that
663  // we will always constrain the dof with the
664  // higher fe index to the one with the lower,
665  // to avoid circular reasoning.
666  DoFIdentities &identities =
667  *vertex_dof_identities[most_dominating_fe_index]
668  [other_fe_index];
669  for (const auto &identity : identities)
670  {
671  const types::global_dof_index master_dof_index =
672  ::internal::DoFAccessorImplementation::
673  Implementation::get_vertex_dof_index(
674  dof_handler,
675  vertex_index,
676  most_dominating_fe_index,
677  identity.first);
678  const types::global_dof_index slave_dof_index =
679  ::internal::DoFAccessorImplementation::
680  Implementation::get_vertex_dof_index(
681  dof_handler,
682  vertex_index,
683  other_fe_index,
684  identity.second);
685 
686  // on subdomain boundaries, we will
687  // encounter invalid DoFs on ghost cells,
688  // for which we have not yet distributed
689  // valid indices. depending on which finte
690  // element is dominating the other on this
691  // interface, we either have to constrain
692  // the valid to the invalid indices, or vice
693  // versa.
694  //
695  // we only store an identity if we are about
696  // to overwrite a valid DoF. we will skip
697  // constraining invalid DoFs for now, and
698  // consider them later in Phase 5.
699  if (slave_dof_index != numbers::invalid_dof_index)
700  {
701  // if the DoF indices of both elements
702  // are already distributed, i.e., both
703  // of these 'fe_indices' are associated
704  // with a locally owned cell, then we
705  // should either not have a dof_identity
706  // yet, or it must come out here to be
707  // exactly as we had computed before
708  if (master_dof_index !=
710  Assert((dof_identities.find(
711  master_dof_index) ==
712  dof_identities.end()) ||
713  (dof_identities[slave_dof_index] ==
714  master_dof_index),
715  ExcInternalError());
716 
717  dof_identities[slave_dof_index] =
718  master_dof_index;
719  }
720  }
721  }
722  }
723  }
724 
725  return dof_identities;
726  }
727 
728 
733  template <int spacedim>
734  static std::map<types::global_dof_index, types::global_dof_index>
735  compute_line_dof_identities(const hp::DoFHandler<1, spacedim> &)
736  {
737  return std::map<types::global_dof_index, types::global_dof_index>();
738  }
739 
740 
741  template <int dim, int spacedim>
742  static std::map<types::global_dof_index, types::global_dof_index>
743  compute_line_dof_identities(
744  const hp::DoFHandler<dim, spacedim> &dof_handler)
745  {
746  std::map<types::global_dof_index, types::global_dof_index>
747  dof_identities;
748 
749  // we will mark lines that we have already treated, so first save and
750  // clear the user flags on lines and later restore them
751  std::vector<bool> user_flags;
752  dof_handler.get_triangulation().save_user_flags_line(user_flags);
753  const_cast<::Triangulation<dim, spacedim> &>(
754  dof_handler.get_triangulation())
755  .clear_user_flags_line();
756 
757  // An implementation of the algorithm described in the hp paper,
758  // including the modification mentioned later in the "complications in
759  // 3-d" subsections
760  //
761  // as explained there, we do something only if there are exactly 2
762  // finite elements associated with an object. if there is only one,
763  // then there is nothing to do anyway, and if there are 3 or more,
764  // then we can get into trouble. note that this only happens for lines
765  // in 3d and higher, and for quads only in 4d and higher, so this
766  // isn't a particularly frequent case
767  //
768  // there is one case, however, that we would like to handle (see, for
769  // example, the hp/crash_15 testcase): if we have
770  // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
771  // then we should be able to handle this because we can simply unify
772  // *all* dofs, not only a some. so what we do is to first treat all
773  // pairs of finite elements that have *identical* dofs, and then only
774  // deal with those that are not identical of which we can handle at
775  // most 2
776  ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
777  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
778 
780  cell = dof_handler.begin_active();
781  cell != dof_handler.end();
782  ++cell)
783  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
784  if (cell->line(l)->user_flag_set() == false)
785  {
787  line = cell->line(l);
788  line->set_user_flag();
789 
790  unsigned int unique_sets_of_dofs =
791  line->n_active_fe_indices();
792 
793  // do a first loop over all sets of dofs and do identity
794  // uniquification
795  const unsigned int n_active_fe_indices =
796  line->n_active_fe_indices();
797  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
798  for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
799  {
800  const unsigned int fe_index_1 =
801  line->nth_active_fe_index(f),
802  fe_index_2 =
803  line->nth_active_fe_index(g);
804 
805  // as described in the hp paper, we only unify on lines
806  // when there are at most two different FE objects
807  // assigned on it.
808  // however, more than two 'active_fe_indices' can be
809  // attached that still fulfill the above criterion,
810  // i.e. when two different FiniteElement objects are
811  // assigned to neighboring cells that map their degrees
812  // of freedom one-to-one.
813  // we cannot verify with certainty if two dofs each of
814  // separate FiniteElement objects actually map
815  // one-to-one. however, checking for the number of
816  // 'dofs_per_line' turned out to be a reasonable
817  // approach, that also works for e.g. two different
818  // FE_Q objects of the same order, from which one is
819  // enhanced by a bubble function that is zero on the
820  // boundary.
821  if ((dof_handler.get_fe(fe_index_1).dofs_per_line ==
822  dof_handler.get_fe(fe_index_2).dofs_per_line) &&
823  (dof_handler.get_fe(fe_index_1).dofs_per_line > 0))
824  {
825  // the number of dofs per line is identical
826  const unsigned int dofs_per_line =
827  dof_handler.get_fe(fe_index_1).dofs_per_line;
828 
829  ensure_existence_of_dof_identities<1>(
830  dof_handler.get_fe(fe_index_1),
831  dof_handler.get_fe(fe_index_2),
832  line_dof_identities[fe_index_1][fe_index_2]);
833  // see if these sets of dofs are identical. the
834  // first condition for this is that indeed there are
835  // n identities
836  if (line_dof_identities[fe_index_1][fe_index_2]
837  ->size() == dofs_per_line)
838  {
839  unsigned int i = 0;
840  for (; i < dofs_per_line; ++i)
841  if (((*(line_dof_identities[fe_index_1]
842  [fe_index_2]))[i]
843  .first != i) &&
844  ((*(line_dof_identities[fe_index_1]
845  [fe_index_2]))[i]
846  .second != i))
847  // not an identity
848  break;
849 
850  if (i == dofs_per_line)
851  {
852  // The line dofs (i.e., the ones interior to
853  // a line) of these two finite elements are
854  // identical. Note that there could be
855  // situations when one element still
856  // dominates another, e.g.: FE_Q(2) x
857  // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
858 
859  --unique_sets_of_dofs;
860 
861  // determine which one of both finite
862  // elements is the dominating one.
863  const std::set<unsigned int> fe_indices{
864  fe_index_1, fe_index_2};
865 
866  unsigned int dominating_fe_index =
867  dof_handler.get_fe_collection()
868  .find_dominating_fe(fe_indices,
869  /*codim=*/dim - 1);
870  unsigned int other_fe_index =
872 
873  if (dominating_fe_index !=
875  other_fe_index =
876  (dominating_fe_index == fe_index_1) ?
877  fe_index_2 :
878  fe_index_1;
879  else
880  {
881  // if we haven't found a dominating
882  // finite element, choose the one with
883  // the lower index to be dominating
884  dominating_fe_index = fe_index_1;
885  other_fe_index = fe_index_2;
886  }
887 
888  for (unsigned int j = 0; j < dofs_per_line;
889  ++j)
890  {
892  master_dof_index = line->dof_index(
893  j, dominating_fe_index);
895  slave_dof_index =
896  line->dof_index(j, other_fe_index);
897 
898  // on subdomain boundaries, we will
899  // encounter invalid DoFs on ghost
900  // cells, for which we have not yet
901  // distributed valid indices. depending
902  // on which finte element is dominating
903  // the other on this interface, we
904  // either have to constrain the valid to
905  // the invalid indices, or vice versa.
906  //
907  // we only store an identity if we are
908  // about to overwrite a valid DoF. we
909  // will skip constraining invalid DoFs
910  // for now, and consider them later in
911  // Phase 5.
912  if (slave_dof_index !=
914  {
915  if (master_dof_index !=
917  {
918  // if master dof was already
919  // constrained, constrain to
920  // that one, otherwise constrain
921  // slave to master
922  if (dof_identities.find(
923  master_dof_index) !=
924  dof_identities.end())
925  {
926  // if the DoF indices of
927  // both elements are already
928  // distributed, i.e., both
929  // of these 'fe_indices' are
930  // associated with a locally
931  // owned cell, then we
932  // should either not have a
933  // dof_identity yet, or it
934  // must come out here to be
935  // exactly as we had
936  // computed before
937  Assert(
938  dof_identities.find(
939  dof_identities
940  [master_dof_index]) ==
941  dof_identities.end(),
942  ExcInternalError());
943 
944  dof_identities
945  [slave_dof_index] =
946  dof_identities
947  [master_dof_index];
948  }
949  else
950  {
951  // see comment above for an
952  // explanation of this
953  // assertion
954  Assert(
955  (dof_identities.find(
956  master_dof_index) ==
957  dof_identities.end()) ||
958  (dof_identities
959  [slave_dof_index] ==
960  master_dof_index),
961  ExcInternalError());
962 
963  dof_identities
964  [slave_dof_index] =
965  master_dof_index;
966  }
967  }
968  else
969  {
970  // set slave_dof to
971  // master_dof_index, which is
972  // invalid
973  dof_identities
974  [slave_dof_index] =
976  }
977  }
978  }
979  }
980  }
981  }
982  }
983 
984  // if at this point, there is only one unique set of dofs
985  // left, then we have taken care of everything above. if there
986  // are two, then we need to deal with them here. if there are
987  // more, then we punt, as described in the paper (and
988  // mentioned above)
989  // TODO: The check for 'dim==2' was inserted by intuition. It
990  // fixes
991  // the previous problems with @ref step_27 "step-27" in 3D. But an
992  // explanation for this is still required, and what we do here
993  // is not what we describe in the paper!.
994  if ((unique_sets_of_dofs == 2) && (dim == 2))
995  {
996  const std::set<unsigned int> fe_indices =
997  line->get_active_fe_indices();
998 
999  // find out which is the most dominating finite element of
1000  // the ones that are used on this line
1001  const unsigned int most_dominating_fe_index =
1002  dof_handler.get_fe_collection().find_dominating_fe(
1003  fe_indices,
1004  /*codim=*/dim - 1);
1005 
1006  // if we found the most dominating element, then use this
1007  // to eliminate some of the degrees of freedom by
1008  // identification. otherwise, the code that computes
1009  // hanging node constraints will have to deal with it by
1010  // computing appropriate constraints along this face/edge
1011  if (most_dominating_fe_index !=
1013  {
1014  // loop over the indices of all the finite elements
1015  // that are not dominating, and identify their dofs to
1016  // the most dominating one
1017  for (const auto &other_fe_index : fe_indices)
1018  if (other_fe_index != most_dominating_fe_index)
1019  {
1020  ensure_existence_of_dof_identities<1>(
1021  dof_handler.get_fe(most_dominating_fe_index),
1022  dof_handler.get_fe(other_fe_index),
1023  line_dof_identities[most_dominating_fe_index]
1024  [other_fe_index]);
1025 
1026  DoFIdentities &identities =
1027  *line_dof_identities[most_dominating_fe_index]
1028  [other_fe_index];
1029  for (const auto &identity : identities)
1030  {
1032  master_dof_index = line->dof_index(
1033  identity.first,
1034  most_dominating_fe_index);
1036  slave_dof_index =
1037  line->dof_index(identity.second,
1038  other_fe_index);
1039 
1040  // on subdomain boundaries, we will
1041  // encounter invalid DoFs on ghost cells,
1042  // for which we have not yet distributed
1043  // valid indices. depending on which finte
1044  // element is dominating the other on this
1045  // interface, we either have to constrain
1046  // the valid to the invalid indices, or vice
1047  // versa.
1048  //
1049  // we only store an identity if we are about
1050  // to overwrite a valid DoF. we will skip
1051  // constraining invalid DoFs for now, and
1052  // consider them later in Phase 5.
1053  if (slave_dof_index !=
1055  {
1056  // if the DoF indices of both elements
1057  // are already distributed, i.e., both
1058  // of these 'fe_indices' are associated
1059  // with a locally owned cell, then we
1060  // should either not have a dof_identity
1061  // yet, or it must come out here to be
1062  // exactly as we had computed before
1063  if (master_dof_index !=
1065  Assert((dof_identities.find(
1066  master_dof_index) ==
1067  dof_identities.end()) ||
1068  (dof_identities
1069  [slave_dof_index] ==
1070  master_dof_index),
1071  ExcInternalError());
1072 
1073  dof_identities[slave_dof_index] =
1074  master_dof_index;
1075  }
1076  }
1077  }
1078  }
1079  }
1080  }
1081 
1082  // finally restore the user flags
1083  const_cast<::Triangulation<dim, spacedim> &>(
1084  dof_handler.get_triangulation())
1085  .load_user_flags_line(user_flags);
1086 
1087  return dof_identities;
1088  }
1089 
1090 
1091 
1096  template <int dim, int spacedim>
1097  static std::map<types::global_dof_index, types::global_dof_index>
1098  compute_quad_dof_identities(const hp::DoFHandler<dim, spacedim> &)
1099  {
1100  // this function should only be called for dim<3 where there are
1101  // no quad dof identies. for dim==3, the specialization below should
1102  // take care of it
1103  Assert(dim < 3, ExcInternalError());
1104 
1105  return std::map<types::global_dof_index, types::global_dof_index>();
1106  }
1107 
1108 
1109  template <int spacedim>
1110  static std::map<types::global_dof_index, types::global_dof_index>
1111  compute_quad_dof_identities(
1112  const hp::DoFHandler<3, spacedim> &dof_handler)
1113  {
1114  const int dim = 3;
1115 
1116  std::map<types::global_dof_index, types::global_dof_index>
1117  dof_identities;
1118 
1119 
1120  // we will mark quads that we have already treated, so first
1121  // save and clear the user flags on quads and later restore
1122  // them
1123  std::vector<bool> user_flags;
1124  dof_handler.get_triangulation().save_user_flags_quad(user_flags);
1125  const_cast<::Triangulation<dim, spacedim> &>(
1126  dof_handler.get_triangulation())
1127  .clear_user_flags_quad();
1128 
1129  // An implementation of the algorithm described in the hp
1130  // paper, including the modification mentioned later in the
1131  // "complications in 3-d" subsections
1132  //
1133  // as explained there, we do something only if there are
1134  // exactly 2 finite elements associated with an object. if
1135  // there is only one, then there is nothing to do anyway,
1136  // and if there are 3 or more, then we can get into
1137  // trouble. note that this only happens for lines in 3d and
1138  // higher, and for quads only in 4d and higher, so this
1139  // isn't a particularly frequent case
1140  ::Table<2, std::unique_ptr<DoFIdentities>> quad_dof_identities(
1141  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
1142 
1144  cell = dof_handler.begin_active();
1145  cell != dof_handler.end();
1146  ++cell)
1147  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
1148  if ((cell->quad(q)->user_flag_set() == false) &&
1149  (cell->quad(q)->n_active_fe_indices() == 2))
1150  {
1152  quad = cell->quad(q);
1153  quad->set_user_flag();
1154 
1155  const std::set<unsigned int> fe_indices =
1156  quad->get_active_fe_indices();
1157 
1158  // find out which is the most dominating finite
1159  // element of the ones that are used on this quad
1160  const unsigned int most_dominating_fe_index =
1161  dof_handler.get_fe_collection().find_dominating_fe(
1162  fe_indices,
1163  /*codim=*/dim - 2);
1164 
1165  // if we found the most dominating element, then use
1166  // this to eliminate some of the degrees of freedom
1167  // by identification. otherwise, the code that
1168  // computes hanging node constraints will have to
1169  // deal with it by computing appropriate constraints
1170  // along this face/edge
1171  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
1172  {
1173  // loop over the indices of all the finite
1174  // elements that are not dominating, and
1175  // identify their dofs to the most dominating
1176  // one
1177  for (const auto &other_fe_index : fe_indices)
1178  if (other_fe_index != most_dominating_fe_index)
1179  {
1180  ensure_existence_of_dof_identities<2>(
1181  dof_handler.get_fe(most_dominating_fe_index),
1182  dof_handler.get_fe(other_fe_index),
1183  quad_dof_identities[most_dominating_fe_index]
1184  [other_fe_index]);
1185 
1186  DoFIdentities &identities =
1187  *quad_dof_identities[most_dominating_fe_index]
1188  [other_fe_index];
1189  for (const auto &identity : identities)
1190  {
1191  const types::global_dof_index master_dof_index =
1192  quad->dof_index(identity.first,
1193  most_dominating_fe_index);
1194  const types::global_dof_index slave_dof_index =
1195  quad->dof_index(identity.second,
1196  other_fe_index);
1197 
1198  // we only store an identity if we are about to
1199  // overwrite a valid degree of freedom. we will
1200  // skip invalid degrees of freedom (that are
1201  // associated with ghost cells) for now, and
1202  // consider them later in phase 5.
1203  if (slave_dof_index !=
1205  {
1206  // if the DoF indices of both elements are
1207  // already distributed, i.e., both of these
1208  // 'fe_indices' are associated with a
1209  // locally owned cell, then we should either
1210  // not have a dof_identity yet, or it must
1211  // come out here to be exactly as we had
1212  // computed before
1213  if (master_dof_index !=
1215  Assert(
1216  (dof_identities.find(
1217  master_dof_index) ==
1218  dof_identities.end()) ||
1219  (dof_identities[slave_dof_index] ==
1220  master_dof_index),
1221  ExcInternalError());
1222 
1223  dof_identities[slave_dof_index] =
1224  master_dof_index;
1225  }
1226  }
1227  }
1228  }
1229  }
1230 
1231  // finally restore the user flags
1232  const_cast<::Triangulation<dim, spacedim> &>(
1233  dof_handler.get_triangulation())
1234  .load_user_flags_quad(user_flags);
1235 
1236  return dof_identities;
1237  }
1238 
1239 
1240 
1245  template <int dim, int spacedim>
1246  static void
1247  compute_dof_identities(
1248  const std::vector<
1249  std::map<types::global_dof_index, types::global_dof_index>> &,
1250  const DoFHandler<dim, spacedim> &)
1251  {}
1252 
1253 
1254  template <int dim, int spacedim>
1255  static void
1256  compute_dof_identities(std::vector<std::map<types::global_dof_index,
1258  &all_constrained_indices,
1259  const hp::DoFHandler<dim, spacedim> &dof_handler)
1260  {
1261  Assert(all_constrained_indices.size() == dim, ExcInternalError());
1262 
1263  Threads::TaskGroup<> tasks;
1264 
1265  unsigned int i = 0;
1266  tasks += Threads::new_task([&, i]() {
1267  all_constrained_indices[i] =
1268  compute_vertex_dof_identities(dof_handler);
1269  });
1270 
1271  if (dim > 1)
1272  {
1273  ++i;
1274  tasks += Threads::new_task([&, i]() {
1275  all_constrained_indices[i] =
1276  compute_line_dof_identities(dof_handler);
1277  });
1278  }
1279 
1280  if (dim > 2)
1281  {
1282  ++i;
1283  tasks += Threads::new_task([&, i]() {
1284  all_constrained_indices[i] =
1285  compute_quad_dof_identities(dof_handler);
1286  });
1287  }
1288 
1289  tasks.join_all();
1290  }
1291 
1292 
1293 
1313  template <class DoFHandlerType>
1315  enumerate_dof_indices_for_renumbering(
1316  std::vector<types::global_dof_index> &new_dof_indices,
1317  const std::vector<
1318  std::map<types::global_dof_index, types::global_dof_index>>
1319  &all_constrained_indices,
1320  const DoFHandlerType &)
1321  {
1322  Assert(all_constrained_indices.size() == DoFHandlerType::dimension,
1323  ExcInternalError());
1324 
1325  // first preset the new DoF indices that are identities
1326  for (const auto &constrained_dof_indices : all_constrained_indices)
1327  for (const auto &p : constrained_dof_indices)
1328  if (new_dof_indices[p.first] != numbers::invalid_dof_index)
1329  {
1330  Assert(new_dof_indices[p.first] == enumeration_dof_index,
1331  ExcInternalError());
1332 
1333  new_dof_indices[p.first] = p.second;
1334  }
1335 
1336  // then enumerate the rest
1337  types::global_dof_index next_free_dof = 0;
1338  for (auto &new_dof_index : new_dof_indices)
1339  if (new_dof_index == enumeration_dof_index)
1340  new_dof_index = next_free_dof++;
1341 
1342  // then loop over all those that are constrained and record the
1343  // new dof number for those
1344  for (const auto &constrained_dof_indices : all_constrained_indices)
1345  for (const auto &p : constrained_dof_indices)
1346  if (new_dof_indices[p.first] != numbers::invalid_dof_index)
1347  {
1348  Assert(new_dof_indices[p.first] != enumeration_dof_index,
1349  ExcInternalError());
1350 
1351  if (p.second != numbers::invalid_dof_index)
1352  new_dof_indices[p.first] = new_dof_indices[p.second];
1353  }
1354 
1355  for (const types::global_dof_index new_dof_index : new_dof_indices)
1356  {
1357  (void)new_dof_index;
1358  Assert(new_dof_index != enumeration_dof_index,
1359  ExcInternalError());
1360  Assert(new_dof_index < next_free_dof ||
1361  new_dof_index == numbers::invalid_dof_index,
1362  ExcInternalError());
1363  }
1364 
1365  return next_free_dof;
1366  }
1367 
1368 
1369 
1379  template <int dim, int spacedim>
1381  unify_dof_indices(const DoFHandler<dim, spacedim> &,
1382  const unsigned int n_dofs_before_identification,
1383  const bool)
1384  {
1385  return n_dofs_before_identification;
1386  }
1387 
1388 
1389  template <int dim, int spacedim>
1391  unify_dof_indices(hp::DoFHandler<dim, spacedim> &dof_handler,
1392  const unsigned int n_dofs_before_identification,
1393  const bool check_validity)
1394  {
1395  std::vector<
1396  std::map<types::global_dof_index, types::global_dof_index>>
1397  all_constrained_indices(dim);
1398  compute_dof_identities(all_constrained_indices, dof_handler);
1399 
1400  std::vector<::types::global_dof_index> renumbering(
1401  n_dofs_before_identification, enumeration_dof_index);
1402  const types::global_dof_index n_dofs =
1403  enumerate_dof_indices_for_renumbering(renumbering,
1404  all_constrained_indices,
1405  dof_handler);
1406 
1407  renumber_dofs(renumbering, IndexSet(0), dof_handler, check_validity);
1408 
1409  update_all_active_cell_dof_indices_caches(dof_handler);
1410 
1411  return n_dofs;
1412  }
1413 
1414 
1415 
1420  template <int dim, int spacedim>
1421  static void
1422  merge_invalid_vertex_dofs_on_ghost_interfaces(
1423  hp::DoFHandler<dim, spacedim> &dof_handler)
1424  {
1425  // Note: we may wish to have something here similar to what
1426  // we do for lines and quads, namely that we only identify
1427  // dofs for any fe towards the most dominating one. however,
1428  // it is not clear whether this is actually necessary for
1429  // vertices at all, I can't think of a finite element that
1430  // would make that necessary...
1432  vertex_dof_identities(dof_handler.get_fe_collection().size(),
1433  dof_handler.get_fe_collection().size());
1434 
1435  // mark all vertices on ghost cells
1436  std::vector<bool> include_vertex(
1437  dof_handler.get_triangulation().n_vertices(), false);
1438  if (dynamic_cast<
1440  &dof_handler.get_triangulation()) != nullptr)
1441  for (const auto &cell : dof_handler.active_cell_iterators())
1442  if (cell->is_ghost())
1443  for (unsigned int v = 0;
1444  v < GeometryInfo<dim>::vertices_per_cell;
1445  ++v)
1446  include_vertex[cell->vertex_index(v)] = true;
1447 
1448  // loop over all vertices and see which one we need to work on
1449  for (unsigned int vertex_index = 0;
1450  vertex_index < dof_handler.get_triangulation().n_vertices();
1451  ++vertex_index)
1452  if ((dof_handler.get_triangulation()
1453  .get_used_vertices()[vertex_index] == true) &&
1454  (include_vertex[vertex_index] == true))
1455  {
1456  const unsigned int n_active_fe_indices =
1457  ::internal::DoFAccessorImplementation::Implementation::
1458  n_active_vertex_fe_indices(dof_handler, vertex_index);
1459 
1460  if (n_active_fe_indices > 1)
1461  {
1462  const std::set<unsigned int> fe_indices =
1463  ::internal::DoFAccessorImplementation::
1464  Implementation::get_active_vertex_fe_indices(
1465  dof_handler, vertex_index);
1466 
1467  // find out which is the most dominating finite
1468  // element of the ones that are used on this vertex
1469  const unsigned int most_dominating_fe_index =
1470  dof_handler.get_fe_collection().find_dominating_fe(
1471  fe_indices,
1472  /*codim=*/dim);
1473 
1474  // if we found the most dominating element, then use
1475  // this to eliminate some of the degrees of freedom
1476  // by identification. otherwise, the code that
1477  // computes hanging node constraints will have to
1478  // deal with it by computing appropriate constraints
1479  // along this face/edge
1480  if (most_dominating_fe_index !=
1482  {
1483  // loop over the indices of all the finite
1484  // elements that are not dominating, and
1485  // identify their dofs to the most dominating
1486  // one
1487  for (const auto &other_fe_index : fe_indices)
1488  if (other_fe_index != most_dominating_fe_index)
1489  {
1490  // make sure the entry in the equivalence
1491  // table exists
1492  ensure_existence_of_dof_identities<0>(
1493  dof_handler.get_fe(most_dominating_fe_index),
1494  dof_handler.get_fe(other_fe_index),
1495  vertex_dof_identities[most_dominating_fe_index]
1496  [other_fe_index]);
1497 
1498  // then loop through the identities we
1499  // have. first get the global numbers of the
1500  // dofs we want to identify and make sure they
1501  // are not yet constrained to anything else,
1502  // except for to each other. use the rule that
1503  // we will always constrain the dof with the
1504  // higher fe index to the one with the lower,
1505  // to avoid circular reasoning.
1506  DoFIdentities &identities =
1507  *vertex_dof_identities[most_dominating_fe_index]
1508  [other_fe_index];
1509  for (const auto &identity : identities)
1510  {
1512  master_dof_index = ::internal::
1513  DoFAccessorImplementation::
1514  Implementation::get_vertex_dof_index(
1515  dof_handler,
1516  vertex_index,
1517  most_dominating_fe_index,
1518  identity.first);
1520  slave_dof_index = ::internal::
1521  DoFAccessorImplementation::
1522  Implementation::get_vertex_dof_index(
1523  dof_handler,
1524  vertex_index,
1525  other_fe_index,
1526  identity.second);
1527 
1528  // check if we are on an interface between
1529  // a locally owned and a ghost cell on which
1530  // we need to work on.
1531  //
1532  // all degrees of freedom belonging to
1533  // dominating fe indices or to a processor
1534  // with a higher rank have been set at this
1535  // point (either in Phase 2, or after the
1536  // first ghost exchange in Phase 5). thus,
1537  // we only have to set the indices of
1538  // degrees of freedom that have been
1539  // previously flagged invalid.
1540  if ((slave_dof_index ==
1542  (master_dof_index !=
1544  ::internal::
1545  DoFAccessorImplementation::
1546  Implementation::set_vertex_dof_index(
1547  dof_handler,
1548  vertex_index,
1549  other_fe_index,
1550  identity.second,
1551  master_dof_index);
1552  }
1553  }
1554  }
1555  }
1556  }
1557  }
1558 
1559 
1560 
1565  template <int spacedim>
1566  static void merge_invalid_line_dofs_on_ghost_interfaces(
1568  {}
1569 
1570 
1571  template <int dim, int spacedim>
1572  static void
1573  merge_invalid_line_dofs_on_ghost_interfaces(
1574  hp::DoFHandler<dim, spacedim> &dof_handler)
1575  {
1576  // we will mark lines that we have already treated, so first save and
1577  // clear the user flags on lines and later restore them
1578  std::vector<bool> user_flags;
1579  dof_handler.get_triangulation().save_user_flags_line(user_flags);
1580  const_cast<::Triangulation<dim, spacedim> &>(
1581  dof_handler.get_triangulation())
1582  .clear_user_flags_line();
1583 
1584  // mark all lines on ghost cells
1585  for (const auto &cell : dof_handler.active_cell_iterators())
1586  if (cell->is_ghost())
1587  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
1588  ++l)
1589  cell->line(l)->set_user_flag();
1590 
1591  // An implementation of the algorithm described in the hp paper,
1592  // including the modification mentioned later in the "complications in
1593  // 3-d" subsections
1594  //
1595  // as explained there, we do something only if there are exactly 2
1596  // finite elements associated with an object. if there is only one,
1597  // then there is nothing to do anyway, and if there are 3 or more,
1598  // then we can get into trouble. note that this only happens for lines
1599  // in 3d and higher, and for quads only in 4d and higher, so this
1600  // isn't a particularly frequent case
1601  //
1602  // there is one case, however, that we would like to handle (see, for
1603  // example, the hp/crash_15 testcase): if we have
1604  // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
1605  // then we should be able to handle this because we can simply unify
1606  // *all* dofs, not only a some. so what we do is to first treat all
1607  // pairs of finite elements that have *identical* dofs, and then only
1608  // deal with those that are not identical of which we can handle at
1609  // most 2
1610  ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
1611  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
1612 
1614  cell = dof_handler.begin_active();
1615  cell != dof_handler.end();
1616  ++cell)
1617  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1618  if ((cell->is_locally_owned()) &&
1619  (cell->line(l)->user_flag_set() == true))
1620  {
1622  line = cell->line(l);
1623  line->clear_user_flag();
1624 
1625  unsigned int unique_sets_of_dofs =
1626  line->n_active_fe_indices();
1627 
1628  // do a first loop over all sets of dofs and do identity
1629  // uniquification
1630  const unsigned int n_active_fe_indices =
1631  line->n_active_fe_indices();
1632  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1633  for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
1634  {
1635  const unsigned int fe_index_1 =
1636  line->nth_active_fe_index(f),
1637  fe_index_2 =
1638  line->nth_active_fe_index(g);
1639 
1640  if ((dof_handler.get_fe(fe_index_1).dofs_per_line ==
1641  dof_handler.get_fe(fe_index_2).dofs_per_line) &&
1642  (dof_handler.get_fe(fe_index_1).dofs_per_line > 0))
1643  {
1644  // the number of dofs per line is identical
1645  const unsigned int dofs_per_line =
1646  dof_handler.get_fe(fe_index_1).dofs_per_line;
1647 
1648  ensure_existence_of_dof_identities<1>(
1649  dof_handler.get_fe(fe_index_1),
1650  dof_handler.get_fe(fe_index_2),
1651  line_dof_identities[fe_index_1][fe_index_2]);
1652  // see if these sets of dofs are identical. the
1653  // first condition for this is that indeed there are
1654  // n identities
1655  if (line_dof_identities[fe_index_1][fe_index_2]
1656  ->size() == dofs_per_line)
1657  {
1658  unsigned int i = 0;
1659  for (; i < dofs_per_line; ++i)
1660  if (((*(line_dof_identities[fe_index_1]
1661  [fe_index_2]))[i]
1662  .first != i) &&
1663  ((*(line_dof_identities[fe_index_1]
1664  [fe_index_2]))[i]
1665  .second != i))
1666  // not an identity
1667  break;
1668 
1669  if (i == dofs_per_line)
1670  {
1671  // The line dofs (i.e., the ones interior to
1672  // a line) of these two finite elements are
1673  // identical. Note that there could be
1674  // situations when one element still
1675  // dominates another, e.g.: FE_Q(2) x
1676  // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
1677 
1678  --unique_sets_of_dofs;
1679 
1680  // determine which one of both finite
1681  // elements is the dominating one.
1682  const std::set<unsigned int> fe_indices{
1683  fe_index_1, fe_index_2};
1684 
1685  unsigned int dominating_fe_index =
1686  dof_handler.get_fe_collection()
1687  .find_dominating_fe(fe_indices,
1688  /*codim*/ dim - 1);
1689  unsigned int other_fe_index =
1691 
1692  if (dominating_fe_index !=
1694  other_fe_index =
1695  (dominating_fe_index == fe_index_1) ?
1696  fe_index_2 :
1697  fe_index_1;
1698  else
1699  {
1700  // if we haven't found a dominating
1701  // finite element, choose the one with
1702  // the lower index to be dominating
1703  dominating_fe_index = fe_index_1;
1704  other_fe_index = fe_index_2;
1705  }
1706 
1707  for (unsigned int j = 0; j < dofs_per_line;
1708  ++j)
1709  {
1711  master_dof_index = line->dof_index(
1712  j, dominating_fe_index);
1714  slave_dof_index =
1715  line->dof_index(j, other_fe_index);
1716 
1717  // check if we are on an interface
1718  // between a locally owned and a ghost
1719  // cell on which we need to work on.
1720  //
1721  // all degrees of freedom belonging to
1722  // dominating fe_indices or to a
1723  // processor with a higher rank have
1724  // been set at this point (either in
1725  // Phase 2, or after the first ghost
1726  // exchange in Phase 5). thus, we only
1727  // have to set the indices of degrees
1728  // of freedom that have been previously
1729  // flagged invalid.
1730  if ((slave_dof_index ==
1732  (master_dof_index !=
1734  line->set_dof_index(j,
1735  master_dof_index,
1736  fe_index_2);
1737  }
1738  }
1739  }
1740  }
1741  }
1742 
1743  // if at this point, there is only one unique set of dofs
1744  // left, then we have taken care of everything above. if there
1745  // are two, then we need to deal with them here. if there are
1746  // more, then we punt, as described in the paper (and
1747  // mentioned above)
1748  // TODO: The check for 'dim==2' was inserted by intuition. It
1749  // fixes
1750  // the previous problems with @ref step_27 "step-27" in 3D. But an
1751  // explanation for this is still required, and what we do here
1752  // is not what we describe in the paper!.
1753  if ((unique_sets_of_dofs == 2) && (dim == 2))
1754  {
1755  const std::set<unsigned int> fe_indices =
1756  line->get_active_fe_indices();
1757 
1758  // find out which is the most dominating finite element of
1759  // the ones that are used on this line
1760  const unsigned int most_dominating_fe_index =
1761  dof_handler.get_fe_collection().find_dominating_fe(
1762  fe_indices,
1763  /*codim=*/dim - 1);
1764 
1765  // if we found the most dominating element, then use this
1766  // to eliminate some of the degrees of freedom by
1767  // identification. otherwise, the code that computes
1768  // hanging node constraints will have to deal with it by
1769  // computing appropriate constraints along this face/edge
1770  if (most_dominating_fe_index !=
1772  {
1773  // loop over the indices of all the finite elements
1774  // that are not dominating, and identify their dofs to
1775  // the most dominating one
1776  for (const auto &other_fe_index : fe_indices)
1777  if (other_fe_index != most_dominating_fe_index)
1778  {
1779  ensure_existence_of_dof_identities<1>(
1780  dof_handler.get_fe(most_dominating_fe_index),
1781  dof_handler.get_fe(other_fe_index),
1782  line_dof_identities[most_dominating_fe_index]
1783  [other_fe_index]);
1784 
1785  DoFIdentities &identities =
1786  *line_dof_identities[most_dominating_fe_index]
1787  [other_fe_index];
1788  for (const auto &identity : identities)
1789  {
1791  master_dof_index = line->dof_index(
1792  identity.first,
1793  most_dominating_fe_index);
1795  slave_dof_index =
1796  line->dof_index(identity.second,
1797  other_fe_index);
1798 
1799  // check if we are on an interface between
1800  // a locally owned and a ghost cell on which
1801  // we need to work on.
1802  //
1803  // all degrees of freedom belonging to
1804  // dominating fe indices or to a processor
1805  // with a higher rank have been set at this
1806  // point (either in Phase 2, or after the
1807  // first ghost exchange in Phase 5). thus,
1808  // we only have to set the indices of
1809  // degrees of freedom that have been
1810  // previously flagged invalid.
1811  if ((slave_dof_index ==
1813  (master_dof_index !=
1815  line->set_dof_index(identity.second,
1816  master_dof_index,
1817  other_fe_index);
1818  }
1819  }
1820  }
1821  }
1822  }
1823 
1824  // finally restore the user flags
1825  const_cast<::Triangulation<dim, spacedim> &>(
1826  dof_handler.get_triangulation())
1827  .load_user_flags_line(user_flags);
1828  }
1829 
1830 
1831 
1836  template <int dim, int spacedim>
1837  static void
1838  merge_invalid_quad_dofs_on_ghost_interfaces(
1840  {
1841  // this function should only be called for dim<3 where there are
1842  // no quad dof identies. for dim>=3, the specialization below should
1843  // take care of it
1844  Assert(dim < 3, ExcInternalError());
1845  }
1846 
1847 
1848  template <int spacedim>
1849  static void merge_invalid_quad_dofs_on_ghost_interfaces(
1850  hp::DoFHandler<3, spacedim> &dof_handler)
1851  {
1852  const int dim = 3;
1853 
1854  // we will mark quads that we have already treated, so first
1855  // save and clear the user flags on quads and later restore
1856  // them
1857  std::vector<bool> user_flags;
1858  dof_handler.get_triangulation().save_user_flags_quad(user_flags);
1859  const_cast<::Triangulation<dim, spacedim> &>(
1860  dof_handler.get_triangulation())
1861  .clear_user_flags_quad();
1862 
1863  // mark all quads on ghost cells
1864  for (const auto &cell : dof_handler.active_cell_iterators())
1865  if (cell->is_ghost())
1866  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell;
1867  ++q)
1868  cell->quad(q)->set_user_flag();
1869 
1870  // An implementation of the algorithm described in the hp
1871  // paper, including the modification mentioned later in the
1872  // "complications in 3-d" subsections
1873  //
1874  // as explained there, we do something only if there are
1875  // exactly 2 finite elements associated with an object. if
1876  // there is only one, then there is nothing to do anyway,
1877  // and if there are 3 or more, then we can get into
1878  // trouble. note that this only happens for lines in 3d and
1879  // higher, and for quads only in 4d and higher, so this
1880  // isn't a particularly frequent case
1881  ::Table<2, std::unique_ptr<DoFIdentities>> quad_dof_identities(
1882  dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
1883 
1885  cell = dof_handler.begin_active();
1886  cell != dof_handler.end();
1887  ++cell)
1888  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
1889  if ((cell->is_locally_owned()) &&
1890  (cell->quad(q)->user_flag_set() == true) &&
1891  (cell->quad(q)->n_active_fe_indices() == 2))
1892  {
1894  quad = cell->quad(q);
1895  quad->clear_user_flag();
1896 
1897  const std::set<unsigned int> fe_indices =
1898  quad->get_active_fe_indices();
1899 
1900  // find out which is the most dominating finite
1901  // element of the ones that are used on this quad
1902  const unsigned int most_dominating_fe_index =
1903  dof_handler.get_fe_collection().find_dominating_fe(
1904  fe_indices,
1905  /*codim=*/dim - 2);
1906 
1907  // if we found the most dominating element, then use
1908  // this to eliminate some of the degrees of freedom
1909  // by identification. otherwise, the code that
1910  // computes hanging node constraints will have to
1911  // deal with it by computing appropriate constraints
1912  // along this face/edge
1913  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
1914  {
1915  // loop over the indices of all the finite
1916  // elements that are not dominating, and
1917  // identify their dofs to the most dominating
1918  // one
1919  for (const auto &other_fe_index : fe_indices)
1920  if (other_fe_index != most_dominating_fe_index)
1921  {
1922  ensure_existence_of_dof_identities<2>(
1923  dof_handler.get_fe(most_dominating_fe_index),
1924  dof_handler.get_fe(other_fe_index),
1925  quad_dof_identities[most_dominating_fe_index]
1926  [other_fe_index]);
1927 
1928  DoFIdentities &identities =
1929  *quad_dof_identities[most_dominating_fe_index]
1930  [other_fe_index];
1931  for (const auto &identity : identities)
1932  {
1933  const types::global_dof_index master_dof_index =
1934  quad->dof_index(identity.first,
1935  most_dominating_fe_index);
1936  const types::global_dof_index slave_dof_index =
1937  quad->dof_index(identity.second,
1938  other_fe_index);
1939 
1940  // check if we are on an interface between
1941  // a locally owned and a ghost cell on which
1942  // we need to work on.
1943  //
1944  // all degrees of freedom belonging to
1945  // dominating fe indices or to a processor with
1946  // a higher rank have been set at this point
1947  // (either in Phase 2, or after the first ghost
1948  // exchange in Phase 5). thus, we only have to
1949  // set the indices of degrees of freedom that
1950  // have been previously flagged invalid.
1951  if ((slave_dof_index ==
1953  (master_dof_index !=
1955  quad->set_dof_index(identity.second,
1956  master_dof_index,
1957  other_fe_index);
1958  }
1959  }
1960  }
1961  }
1962 
1963  // finally restore the user flags
1964  const_cast<::Triangulation<dim, spacedim> &>(
1965  dof_handler.get_triangulation())
1966  .load_user_flags_quad(user_flags);
1967  }
1968 
1969 
1970 
1983  template <int dim, int spacedim>
1984  static void
1985  merge_invalid_dof_indices_on_ghost_interfaces(
1986  const DoFHandler<dim, spacedim> &)
1987  {}
1988 
1989 
1990  template <int dim, int spacedim>
1991  static void
1992  merge_invalid_dof_indices_on_ghost_interfaces(
1993  hp::DoFHandler<dim, spacedim> &dof_handler)
1994  {
1995  {
1996  Threads::TaskGroup<> tasks;
1997 
1998  tasks += Threads::new_task([&]() {
1999  merge_invalid_vertex_dofs_on_ghost_interfaces(dof_handler);
2000  });
2001 
2002  if (dim > 1)
2003  {
2004  tasks += Threads::new_task([&]() {
2005  merge_invalid_line_dofs_on_ghost_interfaces(dof_handler);
2006  });
2007  }
2008 
2009  if (dim > 2)
2010  {
2011  tasks += Threads::new_task([&]() {
2012  merge_invalid_quad_dofs_on_ghost_interfaces(dof_handler);
2013  });
2014  }
2015 
2016  tasks.join_all();
2017  }
2018 
2019  update_all_active_cell_dof_indices_caches(dof_handler);
2020  }
2021 
2022 
2023 
2030  template <class DoFHandlerType>
2032  distribute_dofs(const types::subdomain_id subdomain_id,
2033  DoFHandlerType & dof_handler)
2034  {
2035  Assert(dof_handler.get_triangulation().n_levels() > 0,
2036  ExcMessage("Empty triangulation"));
2037 
2038  // Step 1: distribute dofs on all cells, but definitely
2039  // exclude artificial cells
2040  types::global_dof_index next_free_dof = 0;
2041  typename DoFHandlerType::active_cell_iterator
2042  cell = dof_handler.begin_active(),
2043  endc = dof_handler.end();
2044 
2045  for (; cell != endc; ++cell)
2046  if (!cell->is_artificial())
2047  if ((subdomain_id == numbers::invalid_subdomain_id) ||
2048  (cell->subdomain_id() == subdomain_id))
2049  next_free_dof =
2050  Implementation::distribute_dofs_on_cell(dof_handler,
2051  cell,
2052  next_free_dof);
2053 
2054  update_all_active_cell_dof_indices_caches(dof_handler);
2055 
2056  return next_free_dof;
2057  }
2058 
2059 
2060 
2074  template <class DoFHandlerType>
2075  static void
2076  invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering(
2077  std::vector<types::global_dof_index> &renumbering,
2078  const types::subdomain_id subdomain_id,
2079  const DoFHandlerType & dof_handler)
2080  {
2081  std::vector<types::global_dof_index> local_dof_indices;
2082 
2083  for (const auto &cell : dof_handler.active_cell_iterators())
2084  if (cell->is_ghost() && (cell->subdomain_id() < subdomain_id))
2085  {
2086  // we found a neighboring ghost cell whose subdomain
2087  // is "stronger" than our own subdomain
2088 
2089  // delete all dofs that live there and that we have
2090  // previously assigned a number to (i.e. the ones on
2091  // the interface)
2092  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
2093  cell->get_dof_indices(local_dof_indices);
2094  for (const auto &local_dof_index : local_dof_indices)
2095  if (local_dof_index != numbers::invalid_dof_index)
2096  renumbering[local_dof_index] = numbers::invalid_dof_index;
2097  }
2098  }
2099 
2100 
2101 
2102  /* -------------- distribute_mg_dofs functionality ------------- */
2103 
2104 
2112  template <int dim, int spacedim>
2114  distribute_mg_dofs_on_cell(
2115  const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
2116  types::global_dof_index next_free_dof,
2117  const std::integral_constant<int, 1> &)
2118  {
2119  // distribute dofs of vertices
2120  if (cell->get_fe().dofs_per_vertex > 0)
2121  for (unsigned int v = 0; v < GeometryInfo<1>::vertices_per_cell;
2122  ++v)
2123  {
2124  typename DoFHandler<dim, spacedim>::level_cell_iterator
2125  neighbor = cell->neighbor(v);
2126 
2127  if (neighbor.state() == IteratorState::valid)
2128  {
2129  // has neighbor already been processed?
2130  if (neighbor->user_flag_set() &&
2131  (neighbor->level() == cell->level()))
2132  // copy dofs if the neighbor is on the same level (only
2133  // then are mg dofs the same)
2134  {
2135  if (v == 0)
2136  for (unsigned int d = 0;
2137  d < cell->get_fe().dofs_per_vertex;
2138  ++d)
2139  cell->set_mg_vertex_dof_index(
2140  cell->level(),
2141  0,
2142  d,
2143  neighbor->mg_vertex_dof_index(cell->level(),
2144  1,
2145  d));
2146  else
2147  for (unsigned int d = 0;
2148  d < cell->get_fe().dofs_per_vertex;
2149  ++d)
2150  cell->set_mg_vertex_dof_index(
2151  cell->level(),
2152  1,
2153  d,
2154  neighbor->mg_vertex_dof_index(cell->level(),
2155  0,
2156  d));
2157 
2158  // next neighbor
2159  continue;
2160  }
2161  }
2162 
2163  // otherwise: create dofs newly
2164  for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
2165  ++d)
2166  cell->set_mg_vertex_dof_index(cell->level(),
2167  v,
2168  d,
2169  next_free_dof++);
2170  }
2171 
2172  // dofs of line
2173  if (cell->get_fe().dofs_per_line > 0)
2174  for (unsigned int d = 0; d < cell->get_fe().dofs_per_line; ++d)
2175  cell->set_mg_dof_index(cell->level(), d, next_free_dof++);
2176 
2177  // note that this cell has been processed
2178  cell->set_user_flag();
2179 
2180  return next_free_dof;
2181  }
2182 
2183 
2184 
2185  template <int dim, int spacedim>
2187  distribute_mg_dofs_on_cell(
2188  const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
2189  types::global_dof_index next_free_dof,
2190  const std::integral_constant<int, 2> &)
2191  {
2192  if (cell->get_fe().dofs_per_vertex > 0)
2193  // number dofs on vertices
2194  for (unsigned int vertex = 0;
2195  vertex < GeometryInfo<2>::vertices_per_cell;
2196  ++vertex)
2197  // check whether dofs for this
2198  // vertex have been distributed
2199  // (only check the first dof)
2200  if (cell->mg_vertex_dof_index(cell->level(), vertex, 0) ==
2202  for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
2203  ++d)
2204  cell->set_mg_vertex_dof_index(cell->level(),
2205  vertex,
2206  d,
2207  next_free_dof++);
2208 
2209  // for the four sides
2210  if (cell->get_fe().dofs_per_line > 0)
2211  for (unsigned int side = 0; side < GeometryInfo<2>::faces_per_cell;
2212  ++side)
2213  {
2215  cell->line(side);
2216 
2217  // distribute dofs if necessary: check whether line dof is
2218  // already numbered (check only first dof)
2219  if (line->mg_dof_index(cell->level(), 0) ==
2221  // if not: distribute dofs
2222  for (unsigned int d = 0; d < cell->get_fe().dofs_per_line;
2223  ++d)
2224  line->set_mg_dof_index(cell->level(), d, next_free_dof++);
2225  }
2226 
2227 
2228  // dofs of quad
2229  if (cell->get_fe().dofs_per_quad > 0)
2230  for (unsigned int d = 0; d < cell->get_fe().dofs_per_quad; ++d)
2231  cell->set_mg_dof_index(cell->level(), d, next_free_dof++);
2232 
2233 
2234  // note that this cell has been processed
2235  cell->set_user_flag();
2236 
2237  return next_free_dof;
2238  }
2239 
2240 
2241 
2242  template <int dim, int spacedim>
2244  distribute_mg_dofs_on_cell(
2245  const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
2246  types::global_dof_index next_free_dof,
2247  const std::integral_constant<int, 3> &)
2248  {
2249  if (cell->get_fe().dofs_per_vertex > 0)
2250  // number dofs on vertices
2251  for (unsigned int vertex = 0;
2252  vertex < GeometryInfo<3>::vertices_per_cell;
2253  ++vertex)
2254  // check whether dofs for this vertex have been distributed
2255  // (only check the first dof)
2256  if (cell->mg_vertex_dof_index(cell->level(), vertex, 0) ==
2258  for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
2259  ++d)
2260  cell->set_mg_vertex_dof_index(cell->level(),
2261  vertex,
2262  d,
2263  next_free_dof++);
2264 
2265  // for the lines
2266  if (cell->get_fe().dofs_per_line > 0)
2267  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
2268  {
2270  cell->line(l);
2271 
2272  // distribute dofs if necessary:
2273  // check whether line dof is already
2274  // numbered (check only first dof)
2275  if (line->mg_dof_index(cell->level(), 0) ==
2277  // if not: distribute dofs
2278  for (unsigned int d = 0; d < cell->get_fe().dofs_per_line;
2279  ++d)
2280  line->set_mg_dof_index(cell->level(), d, next_free_dof++);
2281  }
2282 
2283  // for the quads
2284  if (cell->get_fe().dofs_per_quad > 0)
2285  for (unsigned int q = 0; q < GeometryInfo<3>::quads_per_cell; ++q)
2286  {
2288  cell->quad(q);
2289 
2290  // distribute dofs if necessary:
2291  // check whether line dof is already
2292  // numbered (check only first dof)
2293  if (quad->mg_dof_index(cell->level(), 0) ==
2295  // if not: distribute dofs
2296  for (unsigned int d = 0; d < cell->get_fe().dofs_per_quad;
2297  ++d)
2298  quad->set_mg_dof_index(cell->level(), d, next_free_dof++);
2299  }
2300 
2301 
2302  // dofs of cell
2303  if (cell->get_fe().dofs_per_hex > 0)
2304  for (unsigned int d = 0; d < cell->get_fe().dofs_per_hex; ++d)
2305  cell->set_mg_dof_index(cell->level(), d, next_free_dof++);
2306 
2307 
2308  // note that this cell has been processed
2309  cell->set_user_flag();
2310 
2311  return next_free_dof;
2312  }
2313 
2314 
2315 
2316  // same for the hp::DoFHandler
2317  template <int spacedim>
2319  distribute_mg_dofs_on_cell(
2320  const hp::DoFHandler<1, spacedim> &dof_handler,
2322  & cell,
2323  types::global_dof_index next_free_dof)
2324  {
2325  (void)dof_handler;
2326  (void)cell;
2327  (void)next_free_dof;
2328  return 0;
2329  }
2330 
2331 
2332 
2333  template <int spacedim>
2335  distribute_mg_dofs_on_cell(
2336  const hp::DoFHandler<2, spacedim> &dof_handler,
2338  & cell,
2339  types::global_dof_index next_free_dof)
2340  {
2341  (void)dof_handler;
2342  (void)cell;
2343  (void)next_free_dof;
2344  return 0;
2345  }
2346 
2347 
2348 
2349  template <int spacedim>
2351  distribute_mg_dofs_on_cell(
2352  const hp::DoFHandler<3, spacedim> &dof_handler,
2354  & cell,
2355  types::global_dof_index next_free_dof)
2356  {
2357  (void)dof_handler;
2358  (void)cell;
2359  (void)next_free_dof;
2360  return 0;
2361  }
2362 
2363 
2364 
2365  template <class DoFHandlerType>
2367  distribute_dofs_on_level(const types::subdomain_id level_subdomain_id,
2368  DoFHandlerType & dof_handler,
2369  const unsigned int level)
2370  {
2371  const unsigned int dim = DoFHandlerType::dimension;
2372  const unsigned int spacedim = DoFHandlerType::space_dimension;
2373 
2374  const ::Triangulation<dim, spacedim> &tria =
2375  dof_handler.get_triangulation();
2376  Assert(tria.n_levels() > 0, ExcMessage("Empty triangulation"));
2377  if (level >= tria.n_levels())
2378  return 0; // this is allowed for multigrid
2379 
2380  // Clear user flags because we will need them. But first we save
2381  // them and make sure that we restore them later such that at
2382  // the end of this function the Triangulation will be in the
2383  // same state as it was at the beginning of this function.
2384  std::vector<bool> user_flags;
2385  tria.save_user_flags(user_flags);
2386  const_cast<::Triangulation<dim, spacedim> &>(tria)
2387  .clear_user_flags();
2388 
2389  types::global_dof_index next_free_dof = 0;
2390  typename DoFHandler<dim, spacedim>::level_cell_iterator
2391  cell = dof_handler.begin(level),
2392  endc = dof_handler.end(level);
2393 
2394  for (; cell != endc; ++cell)
2395  if ((level_subdomain_id == numbers::invalid_subdomain_id) ||
2396  (cell->level_subdomain_id() == level_subdomain_id))
2397  next_free_dof =
2398  Implementation::distribute_mg_dofs_on_cell<dim, spacedim>(
2399  cell, next_free_dof, std::integral_constant<int, dim>());
2400 
2401  // finally restore the user flags
2402  const_cast<::Triangulation<dim, spacedim> &>(tria)
2403  .load_user_flags(user_flags);
2404 
2405  return next_free_dof;
2406  }
2407 
2408 
2409 
2410  /* --------------------- renumber_dofs functionality ---------------- */
2411 
2412 
2420  template <int dim, int spacedim>
2421  static void
2422  renumber_vertex_dofs(
2423  const std::vector<types::global_dof_index> &new_numbers,
2424  const IndexSet & indices_we_care_about,
2425  DoFHandler<dim, spacedim> & dof_handler,
2426  const bool check_validity)
2427  {
2428  // we can not use cell iterators in this function since then
2429  // we would renumber the dofs on the interface of two cells
2430  // more than once. Anyway, this way it's not only more
2431  // correct but also faster; note, however, that dof numbers
2432  // may be invalid_dof_index, namely when the appropriate
2433  // vertex/line/etc is unused
2434  for (std::vector<types::global_dof_index>::iterator i =
2435  dof_handler.vertex_dofs.begin();
2436  i != dof_handler.vertex_dofs.end();
2437  ++i)
2438  if (*i != numbers::invalid_dof_index)
2439  *i = (indices_we_care_about.size() == 0) ?
2440  (new_numbers[*i]) :
2441  (new_numbers[indices_we_care_about.index_within_set(*i)]);
2442  else if (check_validity)
2443  // if index is invalid_dof_index: check if this one
2444  // really is unused
2445  Assert(dof_handler.get_triangulation().vertex_used(
2446  (i - dof_handler.vertex_dofs.begin()) /
2447  dof_handler.get_fe().dofs_per_vertex) == false,
2448  ExcInternalError());
2449  }
2450 
2451 
2452 
2460  template <int dim, int spacedim>
2461  static void
2462  renumber_cell_dofs(
2463  const std::vector<types::global_dof_index> &new_numbers,
2464  const IndexSet & indices_we_care_about,
2465  DoFHandler<dim, spacedim> & dof_handler)
2466  {
2467  for (unsigned int level = 0; level < dof_handler.levels.size();
2468  ++level)
2469  for (std::vector<types::global_dof_index>::iterator i =
2470  dof_handler.levels[level]->dof_object.dofs.begin();
2471  i != dof_handler.levels[level]->dof_object.dofs.end();
2472  ++i)
2473  if (*i != numbers::invalid_dof_index)
2474  *i =
2475  ((indices_we_care_about.size() == 0) ?
2476  new_numbers[*i] :
2477  new_numbers[indices_we_care_about.index_within_set(*i)]);
2478  }
2479 
2480 
2481 
2489  template <int spacedim>
2490  static void
2491  renumber_face_dofs(
2492  const std::vector<types::global_dof_index> & /*new_numbers*/,
2493  const IndexSet & /*indices_we_care_about*/,
2494  DoFHandler<1, spacedim> & /*dof_handler*/)
2495  {
2496  // nothing to do in 1d since there are no separate faces
2497  }
2498 
2499 
2500 
2501  template <int spacedim>
2502  static void
2503  renumber_face_dofs(
2504  const std::vector<types::global_dof_index> &new_numbers,
2505  const IndexSet & indices_we_care_about,
2506  DoFHandler<2, spacedim> & dof_handler)
2507  {
2508  // treat dofs on lines
2509  for (std::vector<types::global_dof_index>::iterator i =
2510  dof_handler.faces->lines.dofs.begin();
2511  i != dof_handler.faces->lines.dofs.end();
2512  ++i)
2513  if (*i != numbers::invalid_dof_index)
2514  *i = ((indices_we_care_about.size() == 0) ?
2515  new_numbers[*i] :
2516  new_numbers[indices_we_care_about.index_within_set(*i)]);
2517  }
2518 
2519 
2520 
2521  template <int spacedim>
2522  static void
2523  renumber_face_dofs(
2524  const std::vector<types::global_dof_index> &new_numbers,
2525  const IndexSet & indices_we_care_about,
2526  DoFHandler<3, spacedim> & dof_handler)
2527  {
2528  // treat dofs on lines
2529  for (std::vector<types::global_dof_index>::iterator i =
2530  dof_handler.faces->lines.dofs.begin();
2531  i != dof_handler.faces->lines.dofs.end();
2532  ++i)
2533  if (*i != numbers::invalid_dof_index)
2534  *i = ((indices_we_care_about.size() == 0) ?
2535  new_numbers[*i] :
2536  new_numbers[indices_we_care_about.index_within_set(*i)]);
2537 
2538  // treat dofs on quads
2539  for (std::vector<types::global_dof_index>::iterator i =
2540  dof_handler.faces->quads.dofs.begin();
2541  i != dof_handler.faces->quads.dofs.end();
2542  ++i)
2543  if (*i != numbers::invalid_dof_index)
2544  *i = ((indices_we_care_about.size() == 0) ?
2545  new_numbers[*i] :
2546  new_numbers[indices_we_care_about.index_within_set(*i)]);
2547  }
2548 
2549 
2550 
2551  template <int dim, int spacedim>
2552  static void
2553  renumber_vertex_dofs(
2554  const std::vector<types::global_dof_index> &new_numbers,
2555  const IndexSet & indices_we_care_about,
2556  hp::DoFHandler<dim, spacedim> & dof_handler,
2557  const bool check_validity)
2558  {
2559  for (unsigned int vertex_index = 0;
2560  vertex_index < dof_handler.get_triangulation().n_vertices();
2561  ++vertex_index)
2562  {
2563  const unsigned int n_active_fe_indices =
2564  ::internal::DoFAccessorImplementation::Implementation::
2565  n_active_vertex_fe_indices(dof_handler, vertex_index);
2566 
2567  // if this vertex is unused, then we really ought not to have
2568  // allocated any space for it, i.e., n_active_fe_indices should be
2569  // zero, and there is no space to actually store dof indices for
2570  // this vertex
2571  if (dof_handler.get_triangulation().vertex_used(vertex_index) ==
2572  false)
2573  Assert(n_active_fe_indices == 0, ExcInternalError());
2574 
2575  // otherwise the vertex is used; it may still not hold any dof
2576  // indices if it is located on an artificial cell and not adjacent
2577  // to a ghost cell, but in that case there is simply nothing for
2578  // us to do
2579  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2580  {
2581  const unsigned int fe_index =
2582  ::internal::DoFAccessorImplementation::
2583  Implementation::nth_active_vertex_fe_index(dof_handler,
2584  vertex_index,
2585  f);
2586 
2587  for (unsigned int d = 0;
2588  d < dof_handler.get_fe(fe_index).dofs_per_vertex;
2589  ++d)
2590  {
2591  const types::global_dof_index old_dof_index =
2592  ::internal::DoFAccessorImplementation::
2593  Implementation::get_vertex_dof_index(dof_handler,
2594  vertex_index,
2595  fe_index,
2596  d);
2597 
2598  // if check_validity was set, then we are to verify that
2599  // the previous indices were all valid. this really should
2600  // be the case: we allocated space for these vertex dofs,
2601  // i.e., at least one adjacent cell has a valid
2602  // active_fe_index, so there are DoFs that really live
2603  // on this vertex. if check_validity is set, then we
2604  // must make sure that they have been set to something
2605  // useful
2606  if (check_validity)
2607  Assert(old_dof_index != numbers::invalid_dof_index,
2608  ExcInternalError());
2609 
2610  if (old_dof_index != numbers::invalid_dof_index)
2611  {
2612  // In the following blocks, we first check whether
2613  // we were given an IndexSet of DoFs to touch. If not
2614  // (the first 'if' case here), then we are in the
2615  // sequential case and are allowed to touch all DoFs.
2616  //
2617  // If yes (the 'else' case), then we need to
2618  // distinguish whether the DoF whose number we want to
2619  // touch is in fact locally owned (i.e., is in the
2620  // index set) and then we can actually assign it a new
2621  // number; otherwise, we have encountered a
2622  // non-locally owned DoF for which we don't know the
2623  // new number yet and so set it to an invalid index.
2624  // This will later be fixed up after the first ghost
2625  // exchange phase when we unify hp DoFs on neighboring
2626  // cells.
2627  if (indices_we_care_about.size() == 0)
2629  Implementation::set_vertex_dof_index(
2630  dof_handler,
2631  vertex_index,
2632  fe_index,
2633  d,
2634  new_numbers[old_dof_index]);
2635  else
2636  {
2637  if (indices_we_care_about.is_element(
2638  old_dof_index))
2640  Implementation::set_vertex_dof_index(
2641  dof_handler,
2642  vertex_index,
2643  fe_index,
2644  d,
2645  new_numbers[indices_we_care_about
2646  .index_within_set(
2647  old_dof_index)]);
2648  else
2649  ::internal::DoFAccessorImplementation::
2650  Implementation::set_vertex_dof_index(
2651  dof_handler,
2652  vertex_index,
2653  fe_index,
2654  d,
2656  }
2657  }
2658  }
2659  }
2660  }
2661  }
2662 
2663 
2664 
2665  template <int dim, int spacedim>
2666  static void
2667  renumber_cell_dofs(
2668  const std::vector<types::global_dof_index> &new_numbers,
2669  const IndexSet & indices_we_care_about,
2670  hp::DoFHandler<dim, spacedim> & dof_handler)
2671  {
2673  cell = dof_handler.begin_active();
2674  cell != dof_handler.end();
2675  ++cell)
2676  if (!cell->is_artificial())
2677  {
2678  const unsigned int fe_index = cell->active_fe_index();
2679 
2680  for (unsigned int d = 0;
2681  d < dof_handler.get_fe(fe_index)
2682  .template n_dofs_per_object<dim>();
2683  ++d)
2684  {
2685  const types::global_dof_index old_dof_index =
2686  cell->dof_index(d, fe_index);
2687  if (old_dof_index != numbers::invalid_dof_index)
2688  {
2689  // In the following blocks, we first check whether
2690  // we were given an IndexSet of DoFs to touch. If not
2691  // (the first 'if' case here), then we are in the
2692  // sequential case and are allowed to touch all DoFs.
2693  //
2694  // If yes (the 'else' case), then we need to distinguish
2695  // whether the DoF whose number we want to touch is in
2696  // fact locally owned (i.e., is in the index set) and
2697  // then we can actually assign it a new number;
2698  // otherwise, we have encountered a non-locally owned
2699  // DoF for which we don't know the new number yet and so
2700  // set it to an invalid index. This will later be fixed
2701  // up after the first ghost exchange phase when we unify
2702  // hp DoFs on neighboring cells.
2703  if (indices_we_care_about.size() == 0)
2704  cell->set_dof_index(d,
2705  new_numbers[old_dof_index],
2706  fe_index);
2707  else
2708  {
2709  if (indices_we_care_about.is_element(old_dof_index))
2710  cell->set_dof_index(
2711  d,
2712  new_numbers[indices_we_care_about
2713  .index_within_set(old_dof_index)],
2714  fe_index);
2715  else
2716  cell->set_dof_index(d,
2718  fe_index);
2719  }
2720  }
2721  }
2722  }
2723  }
2724 
2725 
2726 
2727  template <int spacedim>
2728  static void
2729  renumber_face_dofs(
2730  const std::vector<types::global_dof_index> & /*new_numbers*/,
2731  const IndexSet & /*indices_we_care_about*/,
2732  hp::DoFHandler<1, spacedim> & /*dof_handler*/)
2733  {
2734  // nothing to do in 1d since there are no separate faces -- we've
2735  // already taken care of this when dealing with the vertices
2736  }
2737 
2738 
2739 
2740  template <int spacedim>
2741  static void
2742  renumber_face_dofs(
2743  const std::vector<types::global_dof_index> &new_numbers,
2744  const IndexSet & indices_we_care_about,
2745  hp::DoFHandler<2, spacedim> & dof_handler)
2746  {
2747  const unsigned int dim = 2;
2748 
2749  // deal with DoFs on lines
2750  {
2751  // save user flags on lines so we can use them to mark lines
2752  // we've already treated
2753  std::vector<bool> saved_line_user_flags;
2754  const_cast<::Triangulation<dim, spacedim> &>(
2755  dof_handler.get_triangulation())
2756  .save_user_flags_line(saved_line_user_flags);
2757  const_cast<::Triangulation<dim, spacedim> &>(
2758  dof_handler.get_triangulation())
2759  .clear_user_flags_line();
2760 
2762  cell = dof_handler.begin_active();
2763  cell != dof_handler.end();
2764  ++cell)
2765  if (!cell->is_artificial())
2766  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
2767  ++l)
2768  if (cell->line(l)->user_flag_set() == false)
2769  {
2770  const typename hp::DoFHandler<dim,
2771  spacedim>::line_iterator
2772  line = cell->line(l);
2773  line->set_user_flag();
2774 
2775  const unsigned int n_active_fe_indices =
2776  line->n_active_fe_indices();
2777 
2778  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2779  {
2780  const unsigned int fe_index =
2781  line->nth_active_fe_index(f);
2782 
2783  for (unsigned int d = 0;
2784  d < dof_handler.get_fe(fe_index).dofs_per_line;
2785  ++d)
2786  {
2787  const types::global_dof_index old_dof_index =
2788  line->dof_index(d, fe_index);
2789  if (old_dof_index != numbers::invalid_dof_index)
2790  {
2791  // In the following blocks, we first check
2792  // whether we were given an IndexSet of DoFs
2793  // to touch. If not (the first 'if' case
2794  // here), then we are in the sequential case
2795  // and are allowed to touch all DoFs.
2796  //
2797  // If yes (the 'else' case), then we need to
2798  // distinguish whether the DoF whose number we
2799  // want to touch is in fact locally owned
2800  // (i.e., is in the index set) and then we can
2801  // actually assign it a new number; otherwise,
2802  // we have encountered a non-locally owned DoF
2803  // for which we don't know the new number yet
2804  // and so set it to an invalid index. This
2805  // will later be fixed up after the first
2806  // ghost exchange phase when we unify hp DoFs
2807  // on neighboring cells.
2808  if (indices_we_care_about.size() == 0)
2809  line->set_dof_index(
2810  d, new_numbers[old_dof_index], fe_index);
2811  else
2812  {
2813  if (indices_we_care_about.is_element(
2814  old_dof_index))
2815  line->set_dof_index(
2816  d,
2817  new_numbers[indices_we_care_about
2819  old_dof_index)],
2820  fe_index);
2821  else
2822  line->set_dof_index(
2823  d,
2825  fe_index);
2826  }
2827  }
2828  }
2829  }
2830  }
2831 
2832  // at the end, restore the user
2833  // flags for the lines
2834  const_cast<::Triangulation<dim, spacedim> &>(
2835  dof_handler.get_triangulation())
2836  .load_user_flags_line(saved_line_user_flags);
2837  }
2838  }
2839 
2840 
2841 
2842  template <int spacedim>
2843  static void
2844  renumber_face_dofs(
2845  const std::vector<types::global_dof_index> &new_numbers,
2846  const IndexSet & indices_we_care_about,
2847  hp::DoFHandler<3, spacedim> & dof_handler)
2848  {
2849  const unsigned int dim = 3;
2850 
2851  // deal with DoFs on lines
2852  {
2853  // save user flags on lines so we can use them to mark lines
2854  // we've already treated
2855  std::vector<bool> saved_line_user_flags;
2856  const_cast<::Triangulation<dim, spacedim> &>(
2857  dof_handler.get_triangulation())
2858  .save_user_flags_line(saved_line_user_flags);
2859  const_cast<::Triangulation<dim, spacedim> &>(
2860  dof_handler.get_triangulation())
2861  .clear_user_flags_line();
2862 
2864  cell = dof_handler.begin_active();
2865  cell != dof_handler.end();
2866  ++cell)
2867  if (!cell->is_artificial())
2868  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
2869  ++l)
2870  if (cell->line(l)->user_flag_set() == false)
2871  {
2872  const typename hp::DoFHandler<dim,
2873  spacedim>::line_iterator
2874  line = cell->line(l);
2875  line->set_user_flag();
2876 
2877  const unsigned int n_active_fe_indices =
2878  line->n_active_fe_indices();
2879 
2880  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2881  {
2882  const unsigned int fe_index =
2883  line->nth_active_fe_index(f);
2884 
2885  for (unsigned int d = 0;
2886  d < dof_handler.get_fe(fe_index).dofs_per_line;
2887  ++d)
2888  {
2889  const types::global_dof_index old_dof_index =
2890  line->dof_index(d, fe_index);
2891  if (old_dof_index != numbers::invalid_dof_index)
2892  {
2893  // In the following blocks, we first check
2894  // whether we were given an IndexSet of DoFs
2895  // to touch. If not (the first 'if' case
2896  // here), then we are in the sequential case
2897  // and are allowed to touch all DoFs.
2898  //
2899  // If yes (the 'else' case), then we need to
2900  // distinguish whether the DoF whose number we
2901  // want to touch is in fact locally owned
2902  // (i.e., is in the index set) and then we can
2903  // actually assign it a new number; otherwise,
2904  // we have encountered a non-locally owned DoF
2905  // for which we don't know the new number yet
2906  // and so set it to an invalid index. This
2907  // will later be fixed up after the first
2908  // ghost exchange phase when we unify hp DoFs
2909  // on neighboring cells.
2910  if (indices_we_care_about.size() == 0)
2911  line->set_dof_index(
2912  d, new_numbers[old_dof_index], fe_index);
2913  else if (indices_we_care_about.is_element(
2914  old_dof_index))
2915  line->set_dof_index(
2916  d,
2917  new_numbers[indices_we_care_about
2919  old_dof_index)],
2920  fe_index);
2921  else
2922  line->set_dof_index(
2923  d, numbers::invalid_dof_index, fe_index);
2924  }
2925  }
2926  }
2927  }
2928 
2929  // at the end, restore the user
2930  // flags for the lines
2931  const_cast<::Triangulation<dim, spacedim> &>(
2932  dof_handler.get_triangulation())
2933  .load_user_flags_line(saved_line_user_flags);
2934  }
2935 
2936  // then deal with dofs on quads
2937  {
2938  std::vector<bool> saved_quad_user_flags;
2939  const_cast<::Triangulation<dim, spacedim> &>(
2940  dof_handler.get_triangulation())
2941  .save_user_flags_quad(saved_quad_user_flags);
2942  const_cast<::Triangulation<dim, spacedim> &>(
2943  dof_handler.get_triangulation())
2944  .clear_user_flags_quad();
2945 
2947  cell = dof_handler.begin_active();
2948  cell != dof_handler.end();
2949  ++cell)
2950  if (!cell->is_artificial())
2951  for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell;
2952  ++q)
2953  if (cell->quad(q)->user_flag_set() == false)
2954  {
2955  const typename hp::DoFHandler<dim,
2956  spacedim>::quad_iterator
2957  quad = cell->quad(q);
2958  quad->set_user_flag();
2959 
2960  const unsigned int n_active_fe_indices =
2961  quad->n_active_fe_indices();
2962 
2963  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2964  {
2965  const unsigned int fe_index =
2966  quad->nth_active_fe_index(f);
2967 
2968  for (unsigned int d = 0;
2969  d < dof_handler.get_fe(fe_index).dofs_per_quad;
2970  ++d)
2971  {
2972  const types::global_dof_index old_dof_index =
2973  quad->dof_index(d, fe_index);
2974  if (old_dof_index != numbers::invalid_dof_index)
2975  {
2976  // In the following blocks, we first check
2977  // whether we were given an IndexSet of DoFs
2978  // to touch. If not (the first 'if' case
2979  // here), then we are in the sequential case
2980  // and are allowed to touch all DoFs.
2981  //
2982  // If yes (the 'else' case), then we need to
2983  // distinguish whether the DoF whose number we
2984  // want to touch is in fact locally owned
2985  // (i.e., is in the index set) and then we can
2986  // actually assign it a new number; otherwise,
2987  // we have encountered a non-locally owned DoF
2988  // for which we don't know the new number yet
2989  // and so set it to an invalid index. This
2990  // will later be fixed up after the first
2991  // ghost exchange phase when we unify hp DoFs
2992  // on neighboring cells.
2993  if (indices_we_care_about.size() == 0)
2994  quad->set_dof_index(
2995  d, new_numbers[old_dof_index], fe_index);
2996  else
2997  {
2998  if (indices_we_care_about.is_element(
2999  old_dof_index))
3000  quad->set_dof_index(
3001  d,
3002  new_numbers[indices_we_care_about
3004  old_dof_index)],
3005  fe_index);
3006  else
3007  quad->set_dof_index(
3008  d,
3010  fe_index);
3011  }
3012  }
3013  }
3014  }
3015  }
3016 
3017  // at the end, restore the user flags for the quads
3018  const_cast<::Triangulation<dim, spacedim> &>(
3019  dof_handler.get_triangulation())
3020  .load_user_flags_quad(saved_quad_user_flags);
3021  }
3022  }
3023 
3024 
3025 
3037  template <class DoFHandlerType>
3038  static void
3039  renumber_dofs(const std::vector<types::global_dof_index> &new_numbers,
3040  const IndexSet &indices_we_care_about,
3041  DoFHandlerType &dof_handler,
3042  const bool check_validity)
3043  {
3044  if (DoFHandlerType::dimension == 1)
3045  Assert(indices_we_care_about == IndexSet(0), ExcNotImplemented());
3046 
3047  // renumber DoF indices on vertices, cells, and faces. this
3048  // can be done in parallel because the respective functions
3049  // work on separate data structures
3050  Threads::TaskGroup<> tasks;
3051  tasks += Threads::new_task([&]() {
3052  renumber_vertex_dofs(new_numbers,
3053  indices_we_care_about,
3054  dof_handler,
3055  check_validity);
3056  });
3057  tasks += Threads::new_task([&]() {
3058  renumber_face_dofs(new_numbers, indices_we_care_about, dof_handler);
3059  });
3060  tasks += Threads::new_task([&]() {
3061  renumber_cell_dofs(new_numbers, indices_we_care_about, dof_handler);
3062  });
3063  tasks.join_all();
3064 
3065  // update the cache used for cell dof indices
3066  update_all_active_cell_dof_indices_caches(dof_handler);
3067  }
3068 
3069 
3070 
3071  /* --------------------- renumber_mg_dofs functionality ----------------
3072  */
3073 
3081  template <int dim, int spacedim>
3082  static void
3083  renumber_vertex_mg_dofs(
3084  const std::vector<::types::global_dof_index> &new_numbers,
3085  const IndexSet & indices_we_care_about,
3086  DoFHandler<dim, spacedim> &dof_handler,
3087  const unsigned int level,
3088  const bool check_validity)
3089  {
3090  (void)check_validity;
3091  Assert(level < dof_handler.get_triangulation().n_levels(),
3092  ExcInternalError());
3093 
3094  for (typename std::vector<
3095  typename DoFHandler<dim, spacedim>::MGVertexDoFs>::iterator i =
3096  dof_handler.mg_vertex_dofs.begin();
3097  i != dof_handler.mg_vertex_dofs.end();
3098  ++i)
3099  // if the present vertex lives on the current level
3100  if ((i->get_coarsest_level() <= level) &&
3101  (i->get_finest_level() >= level))
3102  for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_vertex;
3103  ++d)
3104  {
3105  const ::types::global_dof_index idx =
3106  i->get_index(level,
3107  d,
3108  dof_handler.get_fe().dofs_per_vertex);
3109 
3110  if (idx != numbers::invalid_dof_index)
3111  {
3112  Assert(check_validity == false ||
3113  (indices_we_care_about.size() > 0 ?
3114  indices_we_care_about.is_element(idx) :
3115  (idx < new_numbers.size())),
3116  ExcInternalError());
3117  i->set_index(level,
3118  d,
3119  dof_handler.get_fe().dofs_per_vertex,
3120  (indices_we_care_about.size() == 0) ?
3121  (new_numbers[idx]) :
3122  (new_numbers[indices_we_care_about
3123  .index_within_set(idx)]));
3124  }
3125  }
3126  }
3127 
3128 
3129 
3137  template <int dim, int spacedim>
3138  static void
3139  renumber_cell_mg_dofs(
3140  const std::vector<::types::global_dof_index> &new_numbers,
3141  const IndexSet & indices_we_care_about,
3142  DoFHandler<dim, spacedim> &dof_handler,
3143  const unsigned int level)
3144  {
3145  for (std::vector<types::global_dof_index>::iterator i =
3146  dof_handler.mg_levels[level]->dof_object.dofs.begin();
3147  i != dof_handler.mg_levels[level]->dof_object.dofs.end();
3148  ++i)
3149  {
3150  if (*i != numbers::invalid_dof_index)
3151  {
3152  Assert((indices_we_care_about.size() > 0 ?
3153  indices_we_care_about.is_element(*i) :
3154  (*i < new_numbers.size())),
3155  ExcInternalError());
3156  *i =
3157  (indices_we_care_about.size() == 0) ?
3158  (new_numbers[*i]) :
3159  (new_numbers[indices_we_care_about.index_within_set(*i)]);
3160  }
3161  }
3162  }
3163 
3164 
3165 
3173  template <int spacedim>
3174  static void
3175  renumber_face_mg_dofs(
3176  const std::vector<types::global_dof_index> & /*new_numbers*/,
3177  const IndexSet & /*indices_we_care_about*/,
3178  DoFHandler<1, spacedim> & /*dof_handler*/,
3179  const unsigned int /*level*/,
3180  const bool /*check_validity*/)
3181  {
3182  // nothing to do in 1d because there are no separate faces
3183  }
3184 
3185 
3186 
3187  template <int spacedim>
3188  static void
3189  renumber_face_mg_dofs(
3190  const std::vector<::types::global_dof_index> &new_numbers,
3191  const IndexSet & indices_we_care_about,
3192  DoFHandler<2, spacedim> &dof_handler,
3193  const unsigned int level,
3194  const bool check_validity)
3195  {
3196  if (dof_handler.get_fe().dofs_per_line > 0)
3197  {
3198  // save user flags as they will be modified
3199  std::vector<bool> user_flags;
3200  dof_handler.get_triangulation().save_user_flags(user_flags);
3201  const_cast<::Triangulation<2, spacedim> &>(
3202  dof_handler.get_triangulation())
3203  .clear_user_flags();
3204 
3205  // flag all lines adjacent to cells of the current
3206  // level, as those lines logically belong to the same
3207  // level as the cell, at least for for isotropic
3208  // refinement
3209  typename DoFHandler<2, spacedim>::level_cell_iterator cell,
3210  endc = dof_handler.end(level);
3211  for (cell = dof_handler.begin(level); cell != endc; ++cell)
3212  if (cell->level_subdomain_id() !=
3214  for (unsigned int line = 0;
3215  line < GeometryInfo<2>::faces_per_cell;
3216  ++line)
3217  cell->face(line)->set_user_flag();
3218 
3219  for (typename DoFHandler<2, spacedim>::cell_iterator cell =
3220  dof_handler.begin();
3221  cell != dof_handler.end();
3222  ++cell)
3223  for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell;
3224  ++l)
3225  if (cell->line(l)->user_flag_set())
3226  {
3227  for (unsigned int d = 0;
3228  d < dof_handler.get_fe().dofs_per_line;
3229  ++d)
3230  {
3231  const ::types::global_dof_index idx =
3232  cell->line(l)->mg_dof_index(level, d);
3233  if (check_validity)
3235  ExcInternalError());
3236 
3237  if (idx != numbers::invalid_dof_index)
3238  cell->line(l)->set_mg_dof_index(
3239  level,
3240  d,
3241  ((indices_we_care_about.size() == 0) ?
3242  new_numbers[idx] :
3243  new_numbers[indices_we_care_about
3244  .index_within_set(idx)]));
3245  }
3246  cell->line(l)->clear_user_flag();
3247  }
3248  // finally, restore user flags
3249  const_cast<::Triangulation<2, spacedim> &>(
3250  dof_handler.get_triangulation())
3251  .load_user_flags(user_flags);
3252  }
3253  }
3254 
3255 
3256 
3257  template <int spacedim>
3258  static void
3259  renumber_face_mg_dofs(
3260  const std::vector<::types::global_dof_index> &new_numbers,
3261  const IndexSet & indices_we_care_about,
3262  DoFHandler<3, spacedim> &dof_handler,
3263  const unsigned int level,
3264  const bool check_validity)
3265  {
3266  if (dof_handler.get_fe().dofs_per_line > 0 ||
3267  dof_handler.get_fe().dofs_per_quad > 0)
3268  {
3269  // save user flags as they will be modified
3270  std::vector<bool> user_flags;
3271  dof_handler.get_triangulation().save_user_flags(user_flags);
3272  const_cast<::Triangulation<3, spacedim> &>(
3273  dof_handler.get_triangulation())
3274  .clear_user_flags();
3275 
3276  // flag all lines adjacent to cells of the current
3277  // level, as those lines logically belong to the same
3278  // level as the cell, at least for isotropic refinement
3279  typename DoFHandler<3, spacedim>::level_cell_iterator cell,
3280  endc = dof_handler.end(level);
3281  for (cell = dof_handler.begin(level); cell != endc; ++cell)
3282  if (cell->level_subdomain_id() !=
3284  for (unsigned int line = 0;
3285  line < GeometryInfo<3>::lines_per_cell;
3286  ++line)
3287  cell->line(line)->set_user_flag();
3288 
3289  for (typename DoFHandler<3, spacedim>::cell_iterator cell =
3290  dof_handler.begin();
3291  cell != dof_handler.end();
3292  ++cell)
3293  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell;
3294  ++l)
3295  if (cell->line(l)->user_flag_set())
3296  {
3297  for (unsigned int d = 0;
3298  d < dof_handler.get_fe().dofs_per_line;
3299  ++d)
3300  {
3301  const ::types::global_dof_index idx =
3302  cell->line(l)->mg_dof_index(level, d);
3303  if (check_validity)
3305  ExcInternalError());
3306 
3307  if (idx != numbers::invalid_dof_index)
3308  cell->line(l)->set_mg_dof_index(
3309  level,
3310  d,
3311  ((indices_we_care_about.size() == 0) ?
3312  new_numbers[idx] :
3313  new_numbers[indices_we_care_about
3314  .index_within_set(idx)]));
3315  }
3316  cell->line(l)->clear_user_flag();
3317  }
3318 
3319  // flag all quads adjacent to cells of the current level, as
3320  // those quads logically belong to the same level as the cell,
3321  // at least for isotropic refinement
3322  for (cell = dof_handler.begin(level); cell != endc; ++cell)
3323  if (cell->level_subdomain_id() !=
3325  for (unsigned int quad = 0;
3326  quad < GeometryInfo<3>::quads_per_cell;
3327  ++quad)
3328  cell->quad(quad)->set_user_flag();
3329 
3330  for (typename DoFHandler<3, spacedim>::cell_iterator cell =
3331  dof_handler.begin();
3332  cell != dof_handler.end();
3333  ++cell)
3334  for (unsigned int l = 0; l < GeometryInfo<3>::quads_per_cell;
3335  ++l)
3336  if (cell->quad(l)->user_flag_set())
3337  {
3338  for (unsigned int d = 0;
3339  d < dof_handler.get_fe().dofs_per_quad;
3340  ++d)
3341  {
3342  const ::types::global_dof_index idx =
3343  cell->quad(l)->mg_dof_index(level, d);
3344  if (check_validity)
3346  ExcInternalError());
3347 
3348  if (idx != numbers::invalid_dof_index)
3349  cell->quad(l)->set_mg_dof_index(
3350  level,
3351  d,
3352  ((indices_we_care_about.size() == 0) ?
3353  new_numbers[idx] :
3354  new_numbers[indices_we_care_about
3355  .index_within_set(idx)]));
3356  }
3357  cell->quad(l)->clear_user_flag();
3358  }
3359 
3360  // finally, restore user flags
3361  const_cast<::Triangulation<3, spacedim> &>(
3362  dof_handler.get_triangulation())
3363  .load_user_flags(user_flags);
3364  }
3365  }
3366 
3367 
3368 
3369  template <int dim, int spacedim>
3370  static void
3371  renumber_mg_dofs(
3372  const std::vector<::types::global_dof_index> &new_numbers,
3373  const IndexSet & indices_we_care_about,
3374  DoFHandler<dim, spacedim> &dof_handler,
3375  const unsigned int level,
3376  const bool check_validity)
3377  {
3378  Assert(level < dof_handler.get_triangulation().n_global_levels(),
3379  ExcInternalError());
3380 
3381  // renumber DoF indices on vertices, cells, and faces. this
3382  // can be done in parallel because the respective functions
3383  // work on separate data structures
3384  Threads::TaskGroup<> tasks;
3385  tasks += Threads::new_task([&]() {
3386  renumber_vertex_mg_dofs(new_numbers,
3387  indices_we_care_about,
3388  dof_handler,
3389  level,
3390  check_validity);
3391  });
3392  tasks += Threads::new_task([&]() {
3393  renumber_face_mg_dofs(new_numbers,
3394  indices_we_care_about,
3395  dof_handler,
3396  level,
3397  check_validity);
3398  });
3399  tasks += Threads::new_task([&]() {
3400  renumber_cell_mg_dofs(new_numbers,
3401  indices_we_care_about,
3402  dof_handler,
3403  level);
3404  });
3405  tasks.join_all();
3406  }
3407 
3408 
3409 
3410  template <int dim, int spacedim>
3411  static void
3412  renumber_mg_dofs(
3413  const std::vector<::types::global_dof_index> & /*new_numbers*/,
3414  const IndexSet & /*indices_we_care_about*/,
3415  hp::DoFHandler<dim, spacedim> & /*dof_handler*/,
3416  const unsigned int /*level*/,
3417  const bool /*check_validity*/)
3418  {
3419  Assert(false, ExcNotImplemented());
3420  }
3421  };
3422 
3423 
3424 
3425  /* --------------------- class Sequential ---------------- */
3426 
3427 
3428 
3429  template <class DoFHandlerType>
3430  Sequential<DoFHandlerType>::Sequential(DoFHandlerType &dof_handler)
3431  : dof_handler(&dof_handler)
3432  {}
3433 
3434 
3435 
3436  template <class DoFHandlerType>
3437  NumberCache
3439  {
3440  const types::global_dof_index n_initial_dofs =
3441  Implementation::distribute_dofs(numbers::invalid_subdomain_id,
3442  *dof_handler);
3443 
3444  const types::global_dof_index n_dofs =
3445  Implementation::unify_dof_indices(*dof_handler,
3446  n_initial_dofs,
3447  /*check_validity=*/true);
3448 
3449  // return a sequential, complete index set
3450  return NumberCache(n_dofs);
3451  }
3452 
3453 
3454 
3455  template <class DoFHandlerType>
3456  std::vector<NumberCache>
3458  {
3459  std::vector<bool> user_flags;
3460  dof_handler->get_triangulation().save_user_flags(user_flags);
3461 
3462  const_cast<::Triangulation<DoFHandlerType::dimension,
3463  DoFHandlerType::space_dimension> &>(
3464  dof_handler->get_triangulation())
3465  .clear_user_flags();
3466 
3467  std::vector<NumberCache> number_caches;
3468  number_caches.reserve(dof_handler->get_triangulation().n_levels());
3469  for (unsigned int level = 0;
3470  level < dof_handler->get_triangulation().n_levels();
3471  ++level)
3472  {
3473  // first distribute dofs on this level
3474  const types::global_dof_index n_level_dofs =
3475  Implementation::distribute_dofs_on_level(
3476  numbers::invalid_subdomain_id, *dof_handler, level);
3477 
3478  // then add a complete, sequential index set
3479  number_caches.emplace_back(n_level_dofs);
3480  }
3481 
3482  const_cast<::Triangulation<DoFHandlerType::dimension,
3483  DoFHandlerType::space_dimension> &>(
3484  dof_handler->get_triangulation())
3485  .load_user_flags(user_flags);
3486 
3487  return number_caches;
3488  }
3489 
3490 
3491 
3492  template <class DoFHandlerType>
3493  NumberCache
3495  const std::vector<types::global_dof_index> &new_numbers) const
3496  {
3497  Implementation::renumber_dofs(new_numbers,
3498  IndexSet(0),
3499  *dof_handler,
3500  /*check_validity=*/true);
3501 
3502  // return a sequential, complete index set. take into account that the
3503  // number of DoF indices may in fact be smaller than there were before
3504  // if some previously separately numbered dofs have been identified.
3505  // this is, for example, what the hp::DoFHandler does: it first
3506  // enumerates all DoFs on cells independently, and then unifies
3507  // some located at vertices or faces; this leaves us with fewer
3508  // DoFs than there were before, so use the largest index as
3509  // the one to determine the size of the index space
3510  return NumberCache(
3511  *std::max_element(new_numbers.begin(), new_numbers.end()) + 1);
3512  }
3513 
3514 
3515 
3516  template <class DoFHandlerType>
3517  NumberCache
3519  const unsigned int level,
3520  const std::vector<types::global_dof_index> &new_numbers) const
3521  {
3522  Implementation::renumber_mg_dofs(
3523  new_numbers, IndexSet(0), *dof_handler, level, true);
3524 
3525  // return a sequential, complete index set
3526  return NumberCache(new_numbers.size());
3527  }
3528 
3529 
3530  /* --------------------- class ParallelShared ---------------- */
3531 
3532 
3533  template <class DoFHandlerType>
3535  DoFHandlerType &dof_handler)
3536  : dof_handler(&dof_handler)
3537  {}
3538 
3539 
3540 
3541  namespace
3542  {
3551  template <class DoFHandlerType>
3552  std::vector<types::subdomain_id>
3553  get_dof_subdomain_association(const DoFHandlerType & dof_handler,
3554  const types::global_dof_index n_dofs,
3555  const unsigned int n_procs)
3556  {
3557  (void)n_procs;
3558  std::vector<types::subdomain_id> subdomain_association(
3560  std::vector<types::global_dof_index> local_dof_indices;
3561  local_dof_indices.reserve(DoFTools::max_dofs_per_cell(dof_handler));
3562 
3563  // loop over all cells and record which subdomain a DoF belongs to.
3564  // give to the smaller subdomain_id in case it is on an interface
3565  typename DoFHandlerType::active_cell_iterator
3566  cell = dof_handler.begin_active(),
3567  endc = dof_handler.end();
3568  for (; cell != endc; ++cell)
3569  {
3570  // get the owner of the cell; note that we have made sure above
3571  // that all cells are either locally owned or ghosts (not
3572  // artificial), so this call will always yield the true owner
3573  const types::subdomain_id subdomain_id = cell->subdomain_id();
3574  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
3575  local_dof_indices.resize(dofs_per_cell);
3576  cell->get_dof_indices(local_dof_indices);
3577 
3578  // set subdomain ids. if dofs already have their values set then
3579  // they must be on partition interfaces. In that case assign them
3580  // to the processor with the smaller subdomain id.
3581  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3582  if (subdomain_association[local_dof_indices[i]] ==
3584  subdomain_association[local_dof_indices[i]] = subdomain_id;
3585  else if (subdomain_association[local_dof_indices[i]] >
3586  subdomain_id)
3587  {
3588  subdomain_association[local_dof_indices[i]] = subdomain_id;
3589  }
3590  }
3591 
3592  Assert(std::find(subdomain_association.begin(),
3593  subdomain_association.end(),
3595  subdomain_association.end(),
3596  ExcInternalError());
3597 
3598  Assert(*std::max_element(subdomain_association.begin(),
3599  subdomain_association.end()) < n_procs,
3600  ExcInternalError());
3601 
3602  return subdomain_association;
3603  }
3604 
3605 
3612  template <class DoFHandlerType>
3613  std::vector<types::subdomain_id>
3614  get_dof_level_subdomain_association(
3615  const DoFHandlerType & dof_handler,
3616  const types::global_dof_index n_dofs_on_level,
3617  const unsigned int n_procs,
3618  const unsigned int level)
3619  {
3620  (void)n_procs;
3621  std::vector<types::subdomain_id> level_subdomain_association(
3622  n_dofs_on_level, numbers::invalid_subdomain_id);
3623  std::vector<types::global_dof_index> local_dof_indices;
3624  local_dof_indices.reserve(DoFTools::max_dofs_per_cell(dof_handler));
3625 
3626  // loop over all cells and record which subdomain a DoF belongs to.
3627  // interface goes to proccessor with smaller subdomain id
3628  typename DoFHandlerType::cell_iterator cell =
3629  dof_handler.begin(level),
3630  endc = dof_handler.end(level);
3631  for (; cell != endc; ++cell)
3632  {
3633  // get the owner of the cell; note that we have made sure above
3634  // that all cells are either locally owned or ghosts (not
3635  // artificial), so this call will always yield the true owner
3636  const types::subdomain_id level_subdomain_id =
3637  cell->level_subdomain_id();
3638  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
3639  local_dof_indices.resize(dofs_per_cell);
3640  cell->get_mg_dof_indices(local_dof_indices);
3641 
3642  // set level subdomain ids. if dofs already have their values set
3643  // then they must be on partition interfaces. In that case assign
3644  // them to the processor with the smaller subdomain id.
3645  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3646  if (level_subdomain_association[local_dof_indices[i]] ==
3648  level_subdomain_association[local_dof_indices[i]] =
3649  level_subdomain_id;
3650  else if (level_subdomain_association[local_dof_indices[i]] >
3651  level_subdomain_id)
3652  {
3653  level_subdomain_association[local_dof_indices[i]] =
3654  level_subdomain_id;
3655  }
3656  }
3657 
3658  Assert(std::find(level_subdomain_association.begin(),
3659  level_subdomain_association.end(),
3661  level_subdomain_association.end(),
3662  ExcInternalError());
3663 
3664  Assert(*std::max_element(level_subdomain_association.begin(),
3665  level_subdomain_association.end()) < n_procs,
3666  ExcInternalError());
3667 
3668  return level_subdomain_association;
3669  }
3670  } // namespace
3671 
3672 
3673 
3674  template <class DoFHandlerType>
3675  NumberCache
3677  {
3678  const unsigned int dim = DoFHandlerType::dimension;
3679  const unsigned int spacedim = DoFHandlerType::space_dimension;
3680 
3682  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
3683  &this->dof_handler->get_triangulation()));
3684  Assert(tr != nullptr, ExcInternalError());
3685 
3686  const unsigned int n_procs =
3688 
3689  // If the underlying shared::Tria allows artificial cells,
3690  // then save the current set of subdomain ids, and set
3691  // subdomain ids to the "true" owner of each cell. we later
3692  // restore these flags
3693  std::vector<types::subdomain_id> saved_subdomain_ids;
3694  if (tr->with_artificial_cells())
3695  {
3696  saved_subdomain_ids.resize(tr->n_active_cells());
3697 
3698  const std::vector<types::subdomain_id> &true_subdomain_ids =
3700 
3701  for (const auto &cell : tr->active_cell_iterators())
3702  {
3703  const unsigned int index = cell->active_cell_index();
3704  saved_subdomain_ids[index] = cell->subdomain_id();
3705  cell->set_subdomain_id(true_subdomain_ids[index]);
3706  }
3707  }
3708 
3709  // first let the sequential algorithm do its magic. it is going to
3710  // enumerate DoFs on all cells, regardless of owner
3711  const types::global_dof_index n_initial_dofs =
3712  Implementation::distribute_dofs(numbers::invalid_subdomain_id,
3713  *this->dof_handler);
3714 
3715  const types::global_dof_index n_dofs =
3716  Implementation::unify_dof_indices(*this->dof_handler,
3717  n_initial_dofs,
3718  /*check_validity=*/true);
3719 
3720  // then re-enumerate them based on their subdomain association.
3721  // for this, we first have to identify for each current DoF
3722  // index which subdomain they belong to. ideally, we would
3723  // like to call DoFRenumbering::subdomain_wise(), but
3724  // because the NumberCache of the current DoFHandler is not
3725  // fully set up yet, we can't quite do that. also, that
3726  // function has to deal with other kinds of triangulations as
3727  // well, whereas we here know what kind of triangulation
3728  // we have and can simplify the code accordingly
3729  std::vector<types::global_dof_index> new_dof_indices(
3730  n_dofs, enumeration_dof_index);
3731  {
3732  // first get the association of each dof with a subdomain and
3733  // determine the total number of subdomain ids used
3734  const std::vector<types::subdomain_id> subdomain_association =
3735  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3736 
3737  // then renumber the subdomains by first looking at those belonging
3738  // to subdomain 0, then those of subdomain 1, etc. note that the
3739  // algorithm is stable, i.e. if two dofs i,j have i<j and belong to
3740  // the same subdomain, then they will be in this order also after
3741  // reordering
3742  types::global_dof_index next_free_index = 0;
3743  for (types::subdomain_id subdomain = 0; subdomain < n_procs;
3744  ++subdomain)
3745  for (types::global_dof_index i = 0; i < n_dofs; ++i)
3746  if (subdomain_association[i] == subdomain)
3747  {
3748  Assert(new_dof_indices[i] == enumeration_dof_index,
3749  ExcInternalError());
3750  new_dof_indices[i] = next_free_index;
3751  ++next_free_index;
3752  }
3753 
3754  // we should have numbered all dofs
3755  Assert(next_free_index == n_dofs, ExcInternalError());
3756  Assert(std::find(new_dof_indices.begin(),
3757  new_dof_indices.end(),
3758  enumeration_dof_index) == new_dof_indices.end(),
3759  ExcInternalError());
3760  }
3761  // finally do the renumbering. we can use the sequential
3762  // version of the function because we do things on all
3763  // cells and all cells have their subdomain ids and DoFs
3764  // correctly set
3765  Implementation::renumber_dofs(new_dof_indices,
3766  IndexSet(0),
3767  *this->dof_handler,
3768  /*check_validity=*/true);
3769 
3770  // update the number cache. for this, we first have to find the
3771  // subdomain association for each DoF again following renumbering, from
3772  // which we can then compute the IndexSets of locally owned DoFs for all
3773  // processors. all other fields then follow from this
3774  //
3775  // given the way we enumerate degrees of freedom, the locally owned
3776  // ranges must all be contiguous and consecutive. this makes filling
3777  // the IndexSets cheap. an assertion at the top verifies that this
3778  // assumption is true
3779  const std::vector<types::subdomain_id> subdomain_association =
3780  get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3781 
3782  for (unsigned int i = 1; i < n_dofs; ++i)
3783  Assert(subdomain_association[i] >= subdomain_association[i - 1],
3784  ExcInternalError());
3785 
3786  std::vector<IndexSet> locally_owned_dofs_per_processor(
3787  n_procs, IndexSet(n_dofs));
3788  {
3789  // we know that the set of subdomain indices is contiguous from
3790  // the assertion above; find the start and end index for each
3791  // processor, taking into account that sometimes a processor
3792  // may not in fact have any DoFs at all. we do the latter
3793  // by just identifying contiguous ranges of subdomain_ids
3794  // and filling IndexSets for those subdomains; subdomains
3795  // that don't appear will lead to IndexSets that are simply
3796  // never touched and remain empty as initialized above.
3797  unsigned int start_index = 0;
3798  unsigned int end_index = 0;
3799  while (start_index < n_dofs)
3800  {
3801  while ((end_index) < n_dofs &&
3802  (subdomain_association[end_index] ==
3803  subdomain_association[start_index]))
3804  ++end_index;
3805 
3806  // we've now identified a range of same indices. set that
3807  // range in the corresponding IndexSet
3808  if (end_index > start_index)
3809  {
3810  const unsigned int subdomain_owner =
3811  subdomain_association[start_index];
3812  locally_owned_dofs_per_processor[subdomain_owner].add_range(
3813  start_index, end_index);
3814  }
3815 
3816  // then move on to thinking about the next range
3817  start_index = end_index;
3818  }
3819  }
3820 
3821  // finally, restore current subdomain ids
3822  if (tr->with_artificial_cells())
3823  for (const auto &cell : tr->active_cell_iterators())
3824  cell->set_subdomain_id(
3825  saved_subdomain_ids[cell->active_cell_index()]);
3826 
3827  // return a NumberCache object made up from the sets of locally
3828  // owned DoFs
3829  return NumberCache(
3830  locally_owned_dofs_per_processor,
3831  this->dof_handler->get_triangulation().locally_owned_subdomain());
3832  }
3833 
3834 
3835 
3836  template <class DoFHandlerType>
3837  std::vector<NumberCache>
3839  {
3840  const unsigned int dim = DoFHandlerType::dimension;
3841  const unsigned int spacedim = DoFHandlerType::space_dimension;
3842 
3844  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
3845  &this->dof_handler->get_triangulation()));
3846  Assert(tr != nullptr, ExcInternalError());
3847 
3848  const unsigned int n_procs =
3850  const unsigned int n_levels = tr->n_global_levels();
3851 
3852  std::vector<NumberCache> number_caches;
3853  number_caches.reserve(n_levels);
3854 
3855  // We create an index set for each level
3856  for (unsigned int lvl = 0; lvl < n_levels; ++lvl)
3857  {
3858  // If the underlying shared::Tria allows artificial cells,
3859  // then save the current set of level subdomain ids, and set
3860  // subdomain ids to the "true" owner of each cell. we later
3861  // restore these flags
3862  // Note: "allows_artificial_cells" is currently enforced for
3863  // MG computations.
3864  std::vector<types::subdomain_id> saved_level_subdomain_ids;
3865  saved_level_subdomain_ids.resize(tr->n_cells(lvl));
3866  {
3867  typename parallel::shared::Triangulation<dim,
3868  spacedim>::cell_iterator
3869  cell = this->dof_handler->get_triangulation().begin(lvl),
3870  endc = this->dof_handler->get_triangulation().end(lvl);
3871 
3872  const std::vector<types::subdomain_id> &true_level_subdomain_ids =
3874 
3875  for (unsigned int index = 0; cell != endc; ++cell, ++index)
3876  {
3877  saved_level_subdomain_ids[index] = cell->level_subdomain_id();
3878  cell->set_level_subdomain_id(true_level_subdomain_ids[index]);
3879  }
3880  }
3881 
3882  // Next let the sequential algorithm do its magic. it is going to
3883  // enumerate DoFs on all cells on the given level, regardless of
3884  // owner
3885  const types::global_dof_index n_dofs_on_level =
3886  Implementation::distribute_dofs_on_level(
3887  numbers::invalid_subdomain_id, *this->dof_handler, lvl);
3888 
3889  // then re-enumerate them based on their level subdomain
3890  // association. for this, we first have to identify for each current
3891  // DoF index which subdomain they belong to. ideally, we would like
3892  // to call DoFRenumbering::subdomain_wise(), but because the
3893  // NumberCache of the current DoFHandler is not fully set up yet, we
3894  // can't quite do that. also, that function has to deal with other
3895  // kinds of triangulations as well, whereas we here know what kind
3896  // of triangulation we have and can simplify the code accordingly
3897  std::vector<types::global_dof_index> new_dof_indices(
3898  n_dofs_on_level, numbers::invalid_dof_index);
3899  {
3900  // first get the association of each dof with a subdomain and
3901  // determine the total number of subdomain ids used
3902  const std::vector<types::subdomain_id>
3903  level_subdomain_association =
3904  get_dof_level_subdomain_association(*this->dof_handler,
3905  n_dofs_on_level,
3906  n_procs,
3907  lvl);
3908 
3909  // then renumber the subdomains by first looking at those
3910  // belonging to subdomain 0, then those of subdomain 1, etc. note
3911  // that the algorithm is stable, i.e. if two dofs i,j have i<j and
3912  // belong to the same subdomain, then they will be in this order
3913  // also after reordering
3914  types::global_dof_index next_free_index = 0;
3915  for (types::subdomain_id level_subdomain = 0;
3916  level_subdomain < n_procs;
3917  ++level_subdomain)
3918  for (types::global_dof_index i = 0; i < n_dofs_on_level; ++i)
3919  if (level_subdomain_association[i] == level_subdomain)
3920  {
3921  Assert(new_dof_indices[i] == numbers::invalid_dof_index,
3922  ExcInternalError());
3923  new_dof_indices[i] = next_free_index;
3924  ++next_free_index;
3925  }
3926 
3927  // we should have numbered all dofs
3928  Assert(next_free_index == n_dofs_on_level, ExcInternalError());
3929  Assert(std::find(new_dof_indices.begin(),
3930  new_dof_indices.end(),
3932  new_dof_indices.end(),
3933  ExcInternalError());
3934  }
3935 
3936  // finally do the renumbering. we can use the sequential
3937  // version of the function because we do things on all
3938  // cells and all cells have their subdomain ids and DoFs
3939  // correctly set
3940  Implementation::renumber_mg_dofs(
3941  new_dof_indices, IndexSet(0), *this->dof_handler, lvl, true);
3942 
3943  // update the number cache. for this, we first have to find the
3944  // level subdomain association for each DoF again following
3945  // renumbering, from which we can then compute the IndexSets of
3946  // locally owned DoFs for all processors. all other fields then
3947  // follow from this
3948  //
3949  // given the way we enumerate degrees of freedom, the locally owned
3950  // ranges must all be contiguous and consecutive. this makes filling
3951  // the IndexSets cheap. an assertion at the top verifies that this
3952  // assumption is true
3953  const std::vector<types::subdomain_id> level_subdomain_association =
3954  get_dof_level_subdomain_association(*this->dof_handler,
3955  n_dofs_on_level,
3956  n_procs,
3957  lvl);
3958 
3959  for (unsigned int i = 1; i < n_dofs_on_level; ++i)
3960  Assert(level_subdomain_association[i] >=
3961  level_subdomain_association[i - 1],
3962  ExcInternalError());
3963 
3964  std::vector<IndexSet> locally_owned_dofs_per_processor(
3965  n_procs, IndexSet(n_dofs_on_level));
3966  {
3967  // we know that the set of subdomain indices is contiguous from
3968  // the assertion above; find the start and end index for each
3969  // processor, taking into account that sometimes a processor
3970  // may not in fact have any DoFs at all. we do the latter
3971  // by just identifying contiguous ranges of level_subdomain_ids
3972  // and filling IndexSets for those subdomains; subdomains
3973  // that don't appear will lead to IndexSets that are simply
3974  // never touched and remain empty as initialized above.
3975  unsigned int start_index = 0;
3976  unsigned int end_index = 0;
3977  while (start_index < n_dofs_on_level)
3978  {
3979  while ((end_index) < n_dofs_on_level &&
3980  (level_subdomain_association[end_index] ==
3981  level_subdomain_association[start_index]))
3982  ++end_index;
3983 
3984  // we've now identified a range of same indices. set that
3985  // range in the corresponding IndexSet
3986  if (end_index > start_index)
3987  {
3988  const unsigned int level_subdomain_owner =
3989  level_subdomain_association[start_index];
3990  locally_owned_dofs_per_processor[level_subdomain_owner]
3991  .add_range(start_index, end_index);
3992  }
3993 
3994  // then move on to thinking about the next range
3995  start_index = end_index;
3996  }
3997  }
3998 
3999  // finally, restore current level subdomain ids
4000  {
4001  typename parallel::shared::Triangulation<dim,
4002  spacedim>::cell_iterator
4003  cell = this->dof_handler->get_triangulation().begin(lvl),
4004  endc = this->dof_handler->get_triangulation().end(lvl);
4005 
4006  for (unsigned int index = 0; cell != endc; ++cell, ++index)
4007  cell->set_level_subdomain_id(saved_level_subdomain_ids[index]);
4008 
4009  // add NumberCache for current level
4010  number_caches.emplace_back(
4011  NumberCache(locally_owned_dofs_per_processor,
4012  this->dof_handler->get_triangulation()
4013  .locally_owned_subdomain()));
4014  }
4015  }
4016 
4017  return number_caches;
4018  }
4019 
4020 
4021 
4022  template <class DoFHandlerType>
4023  NumberCache
4025  const std::vector<types::global_dof_index> &new_numbers) const
4026  {
4027 #ifndef DEAL_II_WITH_MPI
4028  (void)new_numbers;
4029  Assert(false, ExcNotImplemented());
4030  return NumberCache();
4031 #else
4032  const unsigned int dim = DoFHandlerType::dimension;
4033  const unsigned int spacedim = DoFHandlerType::space_dimension;
4034 
4035  // Similar to distribute_dofs() we need to have a special treatment in
4036  // case artificial cells are present.
4038  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
4039  &this->dof_handler->get_triangulation()));
4040  Assert(tr != nullptr, ExcInternalError());
4041 
4042  typename parallel::shared::Triangulation<dim,
4043  spacedim>::active_cell_iterator
4044  cell = this->dof_handler->get_triangulation().begin_active(),
4045  endc = this->dof_handler->get_triangulation().end();
4046  std::vector<types::subdomain_id> current_subdomain_ids(
4047  tr->n_active_cells());
4048  const std::vector<types::subdomain_id> &true_subdomain_ids =
4050  if (tr->with_artificial_cells())
4051  for (unsigned int index = 0; cell != endc; cell++, index++)
4052  {
4053  current_subdomain_ids[index] = cell->subdomain_id();
4054  cell->set_subdomain_id(true_subdomain_ids[index]);
4055  }
4056 
4057  std::vector<types::global_dof_index> global_gathered_numbers(
4058  this->dof_handler->n_dofs(), 0);
4059  // as we call DoFRenumbering::subdomain_wise (*dof_handler) from
4060  // distribute_dofs(), we need to support sequential-like input.
4061  // Distributed-like input from, for example, component_wise renumbering
4062  // is also supported.
4063  if (new_numbers.size() == this->dof_handler->n_dofs())
4064  {
4065  global_gathered_numbers = new_numbers;
4066  }
4067  else
4068  {
4069  Assert(new_numbers.size() ==
4070  this->dof_handler->locally_owned_dofs().n_elements(),
4071  ExcInternalError());
4072  const unsigned int n_cpu =
4074  std::vector<types::global_dof_index> gathered_new_numbers(
4075  this->dof_handler->n_dofs(), 0);
4077  this->dof_handler->get_triangulation()
4078  .locally_owned_subdomain(),
4079  ExcInternalError())
4080 
4081  // gather new numbers among processors into one vector
4082  {
4083  std::vector<types::global_dof_index> new_numbers_copy(
4084  new_numbers);
4085 
4086  // store the number of elements that are to be received from each
4087  // process
4088  std::vector<int> rcounts(n_cpu);
4089 
4090  types::global_dof_index shift = 0;
4091  // set rcounts based on new_numbers:
4092  int cur_count = new_numbers_copy.size();
4093  int ierr = MPI_Allgather(&cur_count,
4094  1,
4095  MPI_INT,
4096  rcounts.data(),
4097  1,
4098  MPI_INT,
4099  tr->get_communicator());
4100  AssertThrowMPI(ierr);
4101 
4102  // compute the displacements (relative to recvbuf)
4103  // at which to place the incoming data from process i
4104  std::vector<int> displacements(n_cpu);
4105  for (unsigned int i = 0; i < n_cpu; i++)
4106  {
4107  displacements[i] = shift;
4108  shift += rcounts[i];
4109  }
4110  Assert(new_numbers_copy.size() ==
4111  static_cast<unsigned int>(
4113  tr->get_communicator())]),
4114  ExcInternalError());
4115  ierr = MPI_Allgatherv(new_numbers_copy.data(),
4116  new_numbers_copy.size(),
4117  DEAL_II_DOF_INDEX_MPI_TYPE,
4118  gathered_new_numbers.data(),
4119  rcounts.data(),
4120  displacements.data(),
4121  DEAL_II_DOF_INDEX_MPI_TYPE,
4122  tr->get_communicator());
4123  AssertThrowMPI(ierr);
4124  }
4125 
4126  // put new numbers according to the current
4127  // locally_owned_dofs_per_processor IndexSets
4128  types::global_dof_index shift = 0;
4129  // flag_1 and flag_2 are
4130  // used to control that there is a
4131  // one-to-one relation between old and new DoFs.
4132  std::vector<unsigned int> flag_1(this->dof_handler->n_dofs(), 0);
4133  std::vector<unsigned int> flag_2(this->dof_handler->n_dofs(), 0);
4134  for (unsigned int i = 0; i < n_cpu; i++)
4135  {
4136  const IndexSet iset =
4137  this->dof_handler->locally_owned_dofs_per_processor()[i];
4138  for (types::global_dof_index ind = 0; ind < iset.n_elements();
4139  ind++)
4140  {
4141  const types::global_dof_index target =
4142  iset.nth_index_in_set(ind);
4143  const types::global_dof_index value =
4144  gathered_new_numbers[shift + ind];
4145  Assert(target < this->dof_handler->n_dofs(),
4146  ExcInternalError());
4147  Assert(value < this->dof_handler->n_dofs(),
4148  ExcInternalError());
4149  global_gathered_numbers[target] = value;
4150  flag_1[target]++;
4151  flag_2[value]++;
4152  }
4153  shift += iset.n_elements();
4154  }
4155 
4156  Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1,
4157  ExcInternalError());
4158  Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1,
4159  ExcInternalError());
4160  Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1,
4161  ExcInternalError());
4162  Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1,
4163  ExcInternalError());
4164  }
4165 
4166  // let the sequential algorithm do its magic; ignore the
4167  // return type, but reconstruct the number cache based on
4168  // which DoFs each process owns
4169  Implementation::renumber_dofs(global_gathered_numbers,
4170  IndexSet(0),
4171  *this->dof_handler,
4172  /*check_validity=*/true);
4173 
4174  const NumberCache number_cache(
4175  DoFTools::locally_owned_dofs_per_subdomain(*this->dof_handler),
4176  this->dof_handler->get_triangulation().locally_owned_subdomain());
4177 
4178  // restore artificial cells
4179  cell = tr->begin_active();
4180  if (tr->with_artificial_cells())
4181  for (unsigned int index = 0; cell != endc; cell++, index++)
4182  cell->set_subdomain_id(current_subdomain_ids[index]);
4183 
4184  return number_cache;
4185 #endif
4186  }
4187 
4188 
4189 
4190  template <class DoFHandlerType>
4191  NumberCache
4193  const unsigned int /*level*/,
4194  const std::vector<types::global_dof_index> & /*new_numbers*/) const
4195  {
4196  // multigrid is not currently implemented for shared triangulations
4197  Assert(false, ExcNotImplemented());
4198 
4199  return {};
4200  }
4201 
4202 
4203 
4204  /* --------------------- class ParallelDistributed ---------------- */
4205 
4206 #ifdef DEAL_II_WITH_P4EST
4207 
4208  namespace
4209  {
4210  template <int dim, int spacedim>
4211  void
4212  get_mg_dofindices_recursively(
4214  const typename ::internal::p4est::types<dim>::quadrant
4215  &p4est_cell,
4216  const typename DoFHandler<dim, spacedim>::level_cell_iterator
4217  &dealii_cell,
4218  const typename ::internal::p4est::types<dim>::quadrant
4219  & quadrant,
4220  std::vector<::types::global_dof_index> &dof_numbers_and_indices)
4221  {
4222  if (internal::p4est::quadrant_is_equal<dim>(p4est_cell, quadrant))
4223  {
4224  // why would somebody request a cell that is not ours?
4225  Assert(dealii_cell->level_subdomain_id() ==
4226  tria.locally_owned_subdomain(),
4227  ExcInternalError());
4228 
4229  std::vector<::types::global_dof_index> local_dof_indices(
4230  dealii_cell->get_fe().dofs_per_cell);
4231  dealii_cell->get_mg_dof_indices(local_dof_indices);
4232 
4233  dof_numbers_and_indices.push_back(
4234  dealii_cell->get_fe().dofs_per_cell);
4235  dof_numbers_and_indices.insert(dof_numbers_and_indices.end(),
4236  local_dof_indices.begin(),
4237  local_dof_indices.end());
4238  return; // we are done
4239  }
4240 
4241  if (!dealii_cell->has_children())
4242  return;
4243 
4244  if (!internal::p4est::quadrant_is_ancestor<dim>(p4est_cell, quadrant))
4245  return;
4246 
4247  typename ::internal::p4est::types<dim>::quadrant
4249  internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
4250 
4251  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
4252  ++c)
4253  get_mg_dofindices_recursively<dim, spacedim>(
4254  tria,
4255  p4est_child[c],
4256  dealii_cell->child(c),
4257  quadrant,
4258  dof_numbers_and_indices);
4259  }
4260 
4261 
4262 
4263  template <int dim, int spacedim>
4264  void
4265  find_marked_mg_ghost_cells_recursively(
4267  & tria,
4268  const unsigned int tree_index,
4269  const typename DoFHandler<dim, spacedim>::level_cell_iterator
4270  &dealii_cell,
4271  const typename ::internal::p4est::types<dim>::quadrant
4272  &p4est_cell,
4273  std::map<::types::subdomain_id,
4274  std::vector<std::pair<
4275  unsigned int,
4276  typename ::internal::p4est::types<dim>::quadrant>>>
4277  &neighbor_cell_list)
4278  {
4279  // recurse...
4280  if (dealii_cell->has_children())
4281  {
4282  typename ::internal::p4est::types<dim>::quadrant
4284  internal::p4est::init_quadrant_children<dim>(p4est_cell,
4285  p4est_child);
4286 
4287 
4288  for (unsigned int c = 0;
4289  c < GeometryInfo<dim>::max_children_per_cell;
4290  ++c)
4291  find_marked_mg_ghost_cells_recursively<dim, spacedim>(
4292  tria,
4293  tree_index,
4294  dealii_cell->child(c),
4295  p4est_child[c],
4296  neighbor_cell_list);
4297  }
4298 
4299  if (dealii_cell->user_flag_set() &&
4300  dealii_cell->level_subdomain_id() !=
4301  tria.locally_owned_subdomain())
4302  {
4303  neighbor_cell_list[dealii_cell->level_subdomain_id()]
4304  .emplace_back(tree_index, p4est_cell);
4305  }
4306  }
4307 
4308 
4309 
4310  template <int dim, int spacedim>
4311  void
4312  set_mg_dofindices_recursively(
4314  const typename ::internal::p4est::types<dim>::quadrant
4315  &p4est_cell,
4316  const typename DoFHandler<dim, spacedim>::level_cell_iterator
4317  &dealii_cell,
4318  const typename ::internal::p4est::types<dim>::quadrant
4319  & quadrant,
4320  ::types::global_dof_index *dofs)
4321  {
4322  if (internal::p4est::quadrant_is_equal<dim>(p4est_cell, quadrant))
4323  {
4324  Assert(dealii_cell->level_subdomain_id() !=
4325  ::numbers::artificial_subdomain_id,
4326  ExcInternalError());
4327 
4328  // update dof indices of cell
4329  std::vector<::types::global_dof_index> dof_indices(
4330  dealii_cell->get_fe().dofs_per_cell);
4331  dealii_cell->get_mg_dof_indices(dof_indices);
4332 
4333  bool complete = true;
4334  for (unsigned int i = 0; i < dof_indices.size(); ++i)
4335  if (dofs[i] != numbers::invalid_dof_index)
4336  {
4337  Assert((dof_indices[i] == (numbers::invalid_dof_index)) ||
4338  (dof_indices[i] == dofs[i]),
4339  ExcInternalError());
4340  dof_indices[i] = dofs[i];
4341  }
4342  else
4343  complete = false;
4344 
4345  if (!complete)
4346  const_cast<
4347  typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
4348  dealii_cell)
4349  ->set_user_flag();
4350  else
4351  const_cast<
4352  typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
4353  dealii_cell)
4354  ->clear_user_flag();
4355 
4356  const_cast<
4357  typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
4358  dealii_cell)
4359  ->set_mg_dof_indices(dof_indices);
4360  return;
4361  }
4362 
4363  if (!dealii_cell->has_children())
4364  return;
4365 
4366  if (!internal::p4est::quadrant_is_ancestor<dim>(p4est_cell, quadrant))
4367  return;
4368 
4369  typename ::internal::p4est::types<dim>::quadrant
4371  internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);
4372 
4373  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
4374  ++c)
4375  set_mg_dofindices_recursively<dim, spacedim>(
4376  tria, p4est_child[c], dealii_cell->child(c), quadrant, dofs);
4377  }
4378 
4379 
4380 
4381  template <int dim, int spacedim, class DoFHandlerType>
4382  void
4383  communicate_mg_ghost_cells(
4385  & tria,
4386  DoFHandlerType &dof_handler)
4387  {
4388  using QuadrantBufferType = std::vector<
4389  std::pair<unsigned int,
4390  typename ::internal::p4est::types<dim>::quadrant>>;
4391  // build list of cells to request for each neighbor
4392  std::set<::types::subdomain_id> level_ghost_owners =
4393  tria.level_ghost_owners();
4394  std::map<::types::subdomain_id, QuadrantBufferType>
4395  neighbor_cell_list;
4396  for (const auto level_ghost_owner : level_ghost_owners)
4397  neighbor_cell_list[level_ghost_owner] = {};
4398 
4399  for (typename DoFHandlerType::level_cell_iterator cell =
4400  dof_handler.begin(0);
4401  cell != dof_handler.end(0);
4402  ++cell)
4403  {
4404  typename ::internal::p4est::types<dim>::quadrant
4405  p4est_coarse_cell;
4406  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4407 
4408  types::coarse_cell_id coarse_cell_id = 0;
4409  try
4410  {
4411  coarse_cell_id = cell->id().get_coarse_cell_id();
4412  }
4413  catch (...)
4414  {
4415  // In the case of parallel::fullydistributed::Triangulation,
4416  // a dummy cell throws an exception which is caught here.
4417  // We ignore this cell here.
4418  continue;
4419  };
4420 
4421  find_marked_mg_ghost_cells_recursively<dim, spacedim>(
4422  tria,
4423  coarse_cell_id,
4424  cell,
4425  p4est_coarse_cell,
4426  neighbor_cell_list);
4427  }
4428  Assert(level_ghost_owners.size() == neighbor_cell_list.size(),
4429  ExcInternalError());
4430 
4431  //* send our requests:
4432  std::vector<MPI_Request> requests(level_ghost_owners.size());
4433  {
4434  unsigned int idx = 0;
4435  for (const auto &it : neighbor_cell_list)
4436  {
4437  // send the data about the relevant cells
4438  const int ierr =
4439  MPI_Isend(it.second.data(),
4440  it.second.size() * sizeof(it.second[0]),
4441  MPI_BYTE,
4442  it.first,
4443  10101,
4444  tria.get_communicator(),
4445  &requests[idx]);
4446  AssertThrowMPI(ierr);
4447  ++idx;
4448  }
4449  }
4450 
4451  //* receive requests and reply with the ghost indices
4452  std::vector<QuadrantBufferType> quadrant_data_to_send(
4453  level_ghost_owners.size());
4454  std::vector<std::vector<types::global_dof_index>>
4455  send_dof_numbers_and_indices(level_ghost_owners.size());
4456  std::vector<MPI_Request> reply_requests(level_ghost_owners.size());
4457 
4458  for (unsigned int idx = 0; idx < level_ghost_owners.size(); ++idx)
4459  {
4460  MPI_Status status;
4461  int len;
4462  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4463  10101,
4464  tria.get_communicator(),
4465  &status);
4466  AssertThrowMPI(ierr);
4467  ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4468  AssertThrowMPI(ierr);
4469  Assert(len % sizeof(quadrant_data_to_send[idx][0]) == 0,
4470  ExcInternalError());
4471  const unsigned int n_cells =
4472  len / sizeof(quadrant_data_to_send[idx][0]);
4473  quadrant_data_to_send[idx].resize(n_cells);
4474 
4475  ierr = MPI_Recv(quadrant_data_to_send[idx].data(),
4476  len,
4477  MPI_BYTE,
4478  status.MPI_SOURCE,
4479  status.MPI_TAG,
4480  tria.get_communicator(),
4481  &status);
4482  AssertThrowMPI(ierr);
4483 
4484  // store the dof indices for each cell
4485  for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells);
4486  ++c)
4487  {
4488  const auto temp =
4489  CellId(quadrant_data_to_send[idx][c].first, 0, NULL)
4490  .to_cell(tria);
4491 
4492  typename DoFHandlerType::level_cell_iterator cell(
4493  &dof_handler.get_triangulation(),
4494  0,
4495  temp->index(),
4496  &dof_handler);
4497 
4498  typename ::internal::p4est::types<dim>::quadrant
4499  p4est_coarse_cell;
4500  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4501 
4502  get_mg_dofindices_recursively<dim, spacedim>(
4503  tria,
4504  p4est_coarse_cell,
4505  cell,
4506  quadrant_data_to_send[idx][c].second,
4507  send_dof_numbers_and_indices[idx]);
4508  }
4509 
4510  // send reply
4511  ierr = MPI_Isend(send_dof_numbers_and_indices[idx].data(),
4512  send_dof_numbers_and_indices[idx].size(),
4513  DEAL_II_DOF_INDEX_MPI_TYPE,
4514  status.MPI_SOURCE,
4515  10102,
4516  tria.get_communicator(),
4517  &reply_requests[idx]);
4518  AssertThrowMPI(ierr);
4519  }
4520 
4521  //* finally receive the replies
4522  for (unsigned int idx = 0; idx < level_ghost_owners.size(); ++idx)
4523  {
4524  MPI_Status status;
4525  int len;
4526  int ierr = MPI_Probe(MPI_ANY_SOURCE,
4527  10102,
4528  tria.get_communicator(),
4529  &status);
4530  AssertThrowMPI(ierr);
4531  ierr = MPI_Get_count(&status, DEAL_II_DOF_INDEX_MPI_TYPE, &len);
4532  const QuadrantBufferType &quadrants =
4533  neighbor_cell_list[status.MPI_SOURCE];
4534  AssertThrowMPI(ierr);
4535  Assert((len > 0 && !quadrants.empty()) ||
4536  (len == 0 && quadrants.empty()),
4537  ExcInternalError());
4538  std::vector<types::global_dof_index>
4539  receive_dof_numbers_and_indices(len);
4540 
4541  ierr = MPI_Recv(receive_dof_numbers_and_indices.data(),
4542  len,
4543  DEAL_II_DOF_INDEX_MPI_TYPE,
4544  status.MPI_SOURCE,
4545  status.MPI_TAG,
4546  tria.get_communicator(),
4547  &status);
4548  AssertThrowMPI(ierr);
4549 
4550  // set the dof indices for each cell
4552  receive_dof_numbers_and_indices.data();
4553  for (const auto &it : quadrants)
4554  {
4555  const auto temp = CellId(it.first, 0, NULL).to_cell(tria);
4556 
4557  typename DoFHandlerType::level_cell_iterator cell(
4558  &tria, 0, temp->index(), &dof_handler);
4559 
4560  typename ::internal::p4est::types<dim>::quadrant
4561  p4est_coarse_cell;
4562  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);
4563 
4564  Assert(cell->get_fe().dofs_per_cell == dofs[0],
4565  ExcInternalError());
4566 
4567  set_mg_dofindices_recursively<dim, spacedim>(
4568  tria, p4est_coarse_cell, cell, it.second, dofs + 1);
4569  dofs += 1 + dofs[0];
4570  }
4571  Assert(dofs == receive_dof_numbers_and_indices.data() +
4572  receive_dof_numbers_and_indices.size(),
4573  ExcInternalError());
4574  }
4575 
4576  // complete all sends, so that we can safely destroy the
4577  // buffers.
4578  if (requests.size() > 0)
4579  {
4580  const int ierr = MPI_Waitall(requests.size(),
4581  requests.data(),
4582  MPI_STATUSES_IGNORE);
4583  AssertThrowMPI(ierr);
4584  }
4585  if (reply_requests.size() > 0)
4586  {
4587  const int ierr = MPI_Waitall(reply_requests.size(),
4588  reply_requests.data(),
4589  MPI_STATUSES_IGNORE);
4590  AssertThrowMPI(ierr);
4591  }
4592  }
4593 
4594 
4595 
4596  template <int spacedim>
4597  void
4598  communicate_mg_ghost_cells(
4601  {
4602  Assert(false, ExcNotImplemented());
4603  }
4604 
4605 
4606 
4607  template <int spacedim>
4608  void
4609  communicate_mg_ghost_cells(
4612  {
4613  Assert(false, ExcNotImplemented());
4614  }
4615 
4616 
4617 
4636  template <class DoFHandlerType>
4637  void
4638  communicate_dof_indices_on_marked_cells(
4639  const DoFHandlerType &dof_handler,
4640  const std::map<unsigned int, std::set<::types::subdomain_id>> &)
4641  {
4642 # ifndef DEAL_II_WITH_MPI
4643  (void)vertices_with_ghost_neighbors;
4644  Assert(false, ExcNotImplemented());
4645 # else
4646  const unsigned int dim = DoFHandlerType::dimension;
4647  const unsigned int spacedim = DoFHandlerType::space_dimension;
4648 
4649  // define functions that pack data on cells that are ghost cells
4650  // somewhere else, and unpack data on cells where we get information
4651  // from elsewhere
4652  auto pack =
4653  [](const typename DoFHandlerType::active_cell_iterator &cell)
4654  -> std_cxx17::optional<std::vector<types::global_dof_index>> {
4655  Assert(cell->is_locally_owned(), ExcInternalError());
4656 
4657  // first see whether we need to do anything at all on this cell.
4658  // this is determined by whether the user_flag is set on the
4659  // cell that indicates that the *complete* set of DoF indices
4660  // has not been sent
4661  if (cell->user_flag_set())
4662  {
4663  // get dof indices for the current cell
4664  std::vector<types::global_dof_index> local_dof_indices(
4665  cell->get_fe().dofs_per_cell);
4666  cell->get_dof_indices(local_dof_indices);
4667 
4668  // now see if there are dof indices that were previously
4669  // unknown. this can only happen in phase 1, and in
4670  // that case we know that the user flag must have been set
4671  //
4672  // in any case, if the cell *is* complete, we do not
4673  // need to send the data any more in the next phase. indicate
4674  // this by removing the user flag
4675  if (std::find(local_dof_indices.begin(),
4676  local_dof_indices.end(),
4678  local_dof_indices.end())
4679  {
4680  Assert(cell->user_flag_set(), ExcInternalError());
4681  }
4682  else
4683  cell->clear_user_flag();
4684 
4685  return local_dof_indices;
4686  }
4687  else
4688  {
4689  // the fact that the user flag wasn't set means that there is
4690  // nothing we need to send that hasn't been sent so far.
4691  // so return an empty array, but also verify that indeed
4692  // the cell is complete
4693 # ifdef DEBUG
4694  std::vector<types::global_dof_index> local_dof_indices(
4695  cell->get_fe().dofs_per_cell);
4696  cell->get_dof_indices(local_dof_indices);
4697 
4698  const bool is_complete =
4699  (std::find(local_dof_indices.begin(),
4700  local_dof_indices.end(),
4702  local_dof_indices.end());
4703  Assert(is_complete, ExcInternalError());
4704 # endif
4705  return std_cxx17::optional<
4706  std::vector<types::global_dof_index>>();
4707  }
4708  };
4709 
4710  auto unpack =
4711  [](const typename DoFHandlerType::active_cell_iterator &cell,
4712  const std::vector<types::global_dof_index> &received_dof_indices)
4713  -> void {
4714  // this function should only be called on ghost cells, and
4715  // on top of that, only on cells that have not been
4716  // completed -- which we indicate via the user flag.
4717  // check both
4718  Assert(cell->is_ghost(), ExcInternalError());
4719  Assert(cell->user_flag_set(), ExcInternalError());
4720 
4721  // if we just got an incomplete array of DoF indices, then we must
4722  // be in the first ghost exchange and the user flag must have been
4723  // set. we tested that already above.
4724  //
4725  // if we did get a complete array, then we may be in the first
4726  // or second ghost exchange, but in any case we need not exchange
4727  // another time. so delete the user flag
4728  const bool is_complete = (std::find(received_dof_indices.begin(),
4729  received_dof_indices.end(),
4731  received_dof_indices.end());
4732  if (is_complete)
4733  cell->clear_user_flag();
4734 
4735  // in any case, set the DoF indices on this cell. some
4736  // of the ones we received may still be invalid because
4737  // the sending processor did not know them yet, so we
4738  // need to merge the ones we get with those that are
4739  // already set here and may have already been known. for
4740  // those that we already know *and* get, they must obviously
4741  // agree
4742  //
4743  // before getting the local dof indices, we need to update the
4744  // cell dof indices cache because we may have set dof indices
4745  // on a neighboring ghost cell before this one, which may have
4746  // affected the dof indices we know about the current cell
4747  std::vector<types::global_dof_index> local_dof_indices(
4748  cell->get_fe().dofs_per_cell);
4749  cell->update_cell_dof_indices_cache();
4750  cell->get_dof_indices(local_dof_indices);
4751 
4752  for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
4753  if (local_dof_indices[i] == numbers::invalid_dof_index)
4754  local_dof_indices[i] = received_dof_indices[i];
4755  else
4756  // we already know the dof index. check that there
4757  // is no conflict
4758  Assert((received_dof_indices[i] ==
4760  (received_dof_indices[i] == local_dof_indices[i]),
4761  ExcInternalError());
4762 
4763  const_cast<typename DoFHandlerType::active_cell_iterator &>(cell)
4764  ->set_dof_indices(local_dof_indices);
4765  };
4766 
4768  std::vector<types::global_dof_index>,
4769  DoFHandlerType>(dof_handler, pack, unpack);
4770 
4771  // finally update the cell DoF indices caches to make sure
4772  // our internal data structures are consistent
4773  update_all_active_cell_dof_indices_caches(dof_handler);
4774 
4775 
4776  // have a barrier so that sends between two calls to this
4777  // function are not mixed up.
4778  //
4779  // this is necessary because above we just see if there are
4780  // messages and then receive them, without discriminating
4781  // where they come from and whether they were sent in phase
4782  // 1 or 2 (the function is called twice in a row). the need
4783  // for a global communication step like this barrier could
4784  // be avoided by receiving messages specifically from those
4785  // processors from which we expect messages, and by using
4786  // different tags for phase 1 and 2, but the cost of a
4787  // barrier is negligible compared to everything else we do
4788  // here
4789  if (const auto *triangulation = dynamic_cast<
4791  &dof_handler.get_triangulation()))
4792  {
4793  const int ierr = MPI_Barrier(triangulation->get_communicator());
4794  AssertThrowMPI(ierr);
4795  }
4796  else
4797  {
4798  Assert(false,
4799  ExcMessage(
4800  "The function communicate_dof_indices_on_marked_cells() "
4801  "only works with parallel distributed triangulations."));
4802  }
4803 # endif
4804  }
4805 
4806 
4807 
4808  } // namespace
4809 
4810 #endif // DEAL_II_WITH_P4EST
4811 
4812 
4813 
4814  template <class DoFHandlerType>
4816  DoFHandlerType &dof_handler)
4817  : dof_handler(&dof_handler)
4818  {}
4819 
4820 
4821 
4822  template <class DoFHandlerType>
4823  NumberCache
4825  {
4826 #ifndef DEAL_II_WITH_P4EST
4827  Assert(false, ExcNotImplemented());
4828  return NumberCache();
4829 #else
4830  const unsigned int dim = DoFHandlerType::dimension;
4831  const unsigned int spacedim = DoFHandlerType::space_dimension;
4832 
4835  *>(const_cast<::Triangulation<dim, spacedim> *>(
4836  &dof_handler->get_triangulation())));
4837  Assert(triangulation != nullptr, ExcInternalError());
4838 
4839  const types::subdomain_id subdomain_id =
4840  triangulation->locally_owned_subdomain();
4841 
4842 
4843  /*
4844  The following algorithm has a number of stages that are all
4845  documented in the paper that describes the parallel::distributed
4846  functionality:
4847 
4848  1/ locally enumerate dofs on locally owned cells
4849  2/ eliminate dof duplicates on all cells.
4850  un-numerate those that are on interfaces with ghost
4851  cells and that we don't own based on the tie-breaking
4852  criterion. unify dofs afterwards.
4853  3/ unify dofs and re-enumerate the remaining valid ones.
4854  the end result is that we only enumerate locally owned
4855  DoFs
4856  4/ shift indices so that each processor has a unique
4857  range of indices
4858  5/ for all locally owned cells that are ghost
4859  cells somewhere else, send our own DoF indices
4860  to the appropriate set of other processors.
4861  overwrite invalid DoF indices on ghost interfaces
4862  with the corresponding valid ones that we now know.
4863  6/ send DoF indices again to get the correct indices
4864  on ghost cells that we may not have known earlier
4865  */
4866 
4867  // --------- Phase 1: enumerate dofs on locally owned cells
4868  const types::global_dof_index n_initial_local_dofs =
4869  Implementation::distribute_dofs(subdomain_id, *dof_handler);
4870 
4871  // --------- Phase 2: eliminate dof duplicates on all cells:
4872  // - un-numerate dofs on interfaces to ghost cells
4873  // that we don't own
4874  // - in case of hp::DoFHandler, unify dofs
4875  std::vector<::types::global_dof_index> renumbering(
4876  n_initial_local_dofs, enumeration_dof_index);
4877 
4878  // first, we invalidate degrees of freedom that belong to processors
4879  // of a lower rank, from which we will receive the final (and lower)
4880  // degrees of freedom later.
4881  Implementation::
4882  invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering(
4883  renumbering, subdomain_id, *dof_handler);
4884 
4885  // then, we identify DoF duplicates if a hp::DoFHandler is used
4886  std::vector<std::map<types::global_dof_index, types::global_dof_index>>
4887  all_constrained_indices(dim);
4888  Implementation::compute_dof_identities(all_constrained_indices,
4889  *dof_handler);
4890 
4891  // --------- Phase 3: re-enumerate the valid degrees of freedom
4892  // consecutively. thus, we finally receive the
4893  // correct number of locally owned DoFs after
4894  // this step.
4895  //
4896  // the order in which we handle Phases 2 and 3 is important,
4897  // since we want to clarify ownership of degrees of freedom before
4898  // we actually unify and enumerate their indices. otherwise, we could
4899  // end up having a degee of freedom to which only invalid indices will
4900  // be assigned.
4901  const types::global_dof_index n_locally_owned_dofs =
4902  Implementation::enumerate_dof_indices_for_renumbering(
4903  renumbering, all_constrained_indices, *dof_handler);
4904 
4905  // --------- Phase 4: shift indices so that each processor has a unique
4906  // range of indices
4907  ::types::global_dof_index my_shift = 0;
4908  const int ierr =
4909  MPI_Exscan(DEAL_II_MPI_CONST_CAST(&n_locally_owned_dofs),
4910  &my_shift,
4911  1,
4912  DEAL_II_DOF_INDEX_MPI_TYPE,
4913  MPI_SUM,
4914  triangulation->get_communicator());
4915  AssertThrowMPI(ierr);
4916 
4917  // make dof indices globally consecutive
4918  for (auto &new_index : renumbering)
4919  if (new_index != numbers::invalid_dof_index)
4920  new_index += my_shift;
4921 
4922  // now re-enumerate all dofs to this shifted and condensed
4923  // numbering form. we renumber some dofs as invalid, so
4924  // choose the nocheck-version.
4925  Implementation::renumber_dofs(renumbering,
4926  IndexSet(0),
4927  *dof_handler,
4928  /*check_validity=*/false);
4929 
4930  // now a little bit of housekeeping
4931  const ::types::global_dof_index n_global_dofs =
4932  Utilities::MPI::sum(n_locally_owned_dofs,
4933  triangulation->get_communicator());
4934 
4935  NumberCache number_cache;
4936  number_cache.n_global_dofs = n_global_dofs;
4937  number_cache.n_locally_owned_dofs = n_locally_owned_dofs;
4938  number_cache.locally_owned_dofs = IndexSet(n_global_dofs);
4939  number_cache.locally_owned_dofs.add_range(my_shift,
4940  my_shift +
4941  n_locally_owned_dofs);
4942  number_cache.locally_owned_dofs.compress();
4943 
4944  // this ends the phase where we enumerate degrees of freedom on
4945  // each processor. what is missing is communicating DoF indices
4946  // on ghost cells
4947 
4948  // --------- Phase 5: for all locally owned cells that are ghost
4949  // cells somewhere else, send our own DoF indices
4950  // to the appropriate set of other processors
4951  {
4952  std::vector<bool> user_flags;
4953  triangulation->save_user_flags(user_flags);
4954  triangulation->clear_user_flags();
4955 
4956  // figure out which cells are ghost cells on which we have
4957  // to exchange DoF indices
4958  const std::map<unsigned int, std::set<::types::subdomain_id>>
4959  vertices_with_ghost_neighbors =
4960  triangulation->compute_vertices_with_ghost_neighbors();
4961 
4962  // mark all cells that either have to send data (locally
4963  // owned cells that are adjacent to ghost neighbors in some
4964  // way) or receive data (all ghost cells) via the user flags
4965  for (const auto &cell : dof_handler->active_cell_iterators())
4966  if (cell->is_locally_owned())
4967  {
4968  for (unsigned int v = 0;
4969  v < GeometryInfo<dim>::vertices_per_cell;
4970  ++v)
4971  if (vertices_with_ghost_neighbors.find(cell->vertex_index(
4972  v)) != vertices_with_ghost_neighbors.end())
4973  {
4974  cell->set_user_flag();
4975  break;
4976  }
4977  }
4978  else if (cell->is_ghost())
4979  cell->set_user_flag();
4980 
4981 
4982 
4983  // Send and receive cells. After this, only the local cells
4984  // are marked, that received new data. This has to be
4985  // communicated in a second communication step.
4986  //
4987  // as explained in the 'distributed' paper, this has to be
4988  // done twice
4989  communicate_dof_indices_on_marked_cells(
4990  *dof_handler, vertices_with_ghost_neighbors);
4991 
4992  // in case of hp::DoFHandlers, we may have received valid
4993  // indices of degrees of freedom that are dominated by a fe
4994  // object adjacent to a ghost interface.
4995  // thus, we overwrite the remaining invalid indices with
4996  // the valid ones in this step.
4997  Implementation::merge_invalid_dof_indices_on_ghost_interfaces(
4998  *dof_handler);
4999 
5000  // --------- Phase 6: all locally owned cells have their correct
5001  // DoF indices set. however, some ghost cells
5002  // may still have invalid ones. thus, exchange
5003  // one more time.
5004  communicate_dof_indices_on_marked_cells(
5005  *dof_handler, vertices_with_ghost_neighbors);
5006 
5007  // at this point, we must have taken care of the data transfer
5008  // on all cells we had previously marked. verify this
5009 # ifdef DEBUG
5010  for (const auto &cell : dof_handler->active_cell_iterators())
5011  Assert(cell->user_flag_set() == false, ExcInternalError());
5012 # endif
5013 
5014  triangulation->load_user_flags(user_flags);
5015  }
5016 
5017 # ifdef DEBUG
5018  // check that we are really done
5019  {
5020  std::vector<::types::global_dof_index> local_dof_indices;
5021 
5022  for (const auto &cell : dof_handler->active_cell_iterators())
5023  if (!cell->is_artificial())
5024  {
5025  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
5026  cell->get_dof_indices(local_dof_indices);
5027  if (local_dof_indices.end() !=
5028  std::find(local_dof_indices.begin(),
5029  local_dof_indices.end(),
5031  {
5032  if (cell->is_ghost())
5033  {
5034  Assert(false,
5035  ExcMessage(
5036  "A ghost cell ended up with incomplete "
5037  "DoF index information. This should not "
5038  "have happened!"));
5039  }
5040  else
5041  {
5042  Assert(
5043  false,
5044  ExcMessage(
5045  "A locally owned cell ended up with incomplete "
5046  "DoF index information. This should not "
5047  "have happened!"));
5048  }
5049  }
5050  }
5051  }
5052 # endif // DEBUG
5053  return number_cache;
5054 #endif // DEAL_II_WITH_P4EST
5055  }
5056 
5057 
5058 
5059  template <class DoFHandlerType>
5060  std::vector<NumberCache>
5062  {
5063 #ifndef DEAL_II_WITH_P4EST
5064  Assert(false, ExcNotImplemented());
5065  return std::vector<NumberCache>();
5066 #else
5067  const unsigned int dim = DoFHandlerType::dimension;
5068  const unsigned int spacedim = DoFHandlerType::space_dimension;
5069 
5072  *>(const_cast<::Triangulation<dim, spacedim> *>(
5073  &dof_handler->get_triangulation())));
5074  Assert(triangulation != nullptr, ExcInternalError());
5075 
5077  ExcMessage(
5078  "Multigrid DoFs can only be distributed on a parallel "
5079  "Triangulation if the flag construct_multigrid_hierarchy "
5080  "is set in the constructor."));
5081 
5082  // loop over all levels that exist globally (across all
5083  // processors), even if the current processor does not in fact
5084  // have any cells on that level or if the local part of the
5085  // Triangulation has fewer levels. we need to do this because
5086  // we need to communicate across all processors on all levels
5087  const unsigned int n_levels = triangulation->n_global_levels();
5088  std::vector<NumberCache> number_caches;
5089  number_caches.reserve(n_levels);
5090  for (unsigned int level = 0; level < n_levels; ++level)
5091  {
5092  NumberCache level_number_cache;
5093 
5094  //* 1. distribute on own subdomain
5095  const unsigned int n_initial_local_dofs =
5096  Implementation::distribute_dofs_on_level(
5097  triangulation->locally_owned_subdomain(), *dof_handler, level);
5098 
5099  //* 2. iterate over ghostcells and kill dofs that are not
5100  // owned by us
5101  std::vector<::types::global_dof_index> renumbering(
5102  n_initial_local_dofs);
5103  for (::types::global_dof_index i = 0; i < renumbering.size();
5104  ++i)
5105  renumbering[i] = i;
5106 
5107  if (level < triangulation->n_levels())
5108  {
5109  std::vector<::types::global_dof_index> local_dof_indices;
5110 
5111  typename DoFHandlerType::level_cell_iterator
5112  cell = dof_handler->begin(level),
5113  endc = dof_handler->end(level);
5114 
5115  for (; cell != endc; ++cell)
5116  if (cell->level_subdomain_id() !=
5118  (cell->level_subdomain_id() <
5119  triangulation->locally_owned_subdomain()))
5120  {
5121  // we found a neighboring ghost cell whose
5122  // subdomain is "stronger" than our own
5123  // subdomain
5124 
5125  // delete all dofs that live there and that we
5126  // have previously assigned a number to
5127  // (i.e. the ones on the interface)
5128  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
5129  cell->get_mg_dof_indices(local_dof_indices);
5130  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
5131  ++i)
5132  if (local_dof_indices[i] != numbers::invalid_dof_index)
5133  renumbering[local_dof_indices[i]] =
5135  }
5136  }
5137 
5138  level_number_cache.n_locally_owned_dofs = 0;
5139  for (types::global_dof_index &index : renumbering)
5140  if (index != numbers::invalid_dof_index)
5141  index = level_number_cache.n_locally_owned_dofs++;
5142 
5143  //* 3. communicate local dofcount and shift ids to make
5144  // them unique
5145  ::types::global_dof_index my_shift = 0;
5146  int ierr = MPI_Exscan(DEAL_II_MPI_CONST_CAST(
5147  &level_number_cache.n_locally_owned_dofs),
5148  &my_shift,
5149  1,
5150  DEAL_II_DOF_INDEX_MPI_TYPE,
5151  MPI_SUM,
5152  triangulation->get_communicator());
5153  AssertThrowMPI(ierr);
5154 
5155  // The last processor knows about the total number of dofs, so we
5156  // can use a cheaper broadcast rather than an MPI_Allreduce via
5157  // MPI::sum().
5158  level_number_cache.n_global_dofs =
5159  my_shift + level_number_cache.n_locally_owned_dofs;
5160  ierr = MPI_Bcast(&level_number_cache.n_global_dofs,
5161  1,
5162  DEAL_II_DOF_INDEX_MPI_TYPE,
5164  triangulation->get_communicator()) -
5165  1,
5166  triangulation->get_communicator());
5167 
5168  // shift indices
5169  for (types::global_dof_index &index : renumbering)
5170  if (index != numbers::invalid_dof_index)
5171  index += my_shift;
5172 
5173  // now re-enumerate all dofs to this shifted and condensed
5174  // numbering form. we renumber some dofs as invalid, so
5175  // choose the nocheck-version of the function
5176  //
5177  // of course there is nothing for us to renumber if the
5178  // level we are currently dealing with doesn't even exist
5179  // within the current triangulation, so skip renumbering
5180  // in that case
5181  if (level < triangulation->n_levels())
5182  Implementation::renumber_mg_dofs(
5183  renumbering, IndexSet(0), *dof_handler, level, false);
5184 
5185  // now a little bit of housekeeping
5186  level_number_cache.locally_owned_dofs =
5187  IndexSet(level_number_cache.n_global_dofs);
5188  level_number_cache.locally_owned_dofs.add_range(
5189  my_shift, my_shift + level_number_cache.n_locally_owned_dofs);
5190  level_number_cache.locally_owned_dofs.compress();
5191 
5192  number_caches.emplace_back(level_number_cache);
5193  }
5194 
5195 
5196  //* communicate ghost DoFs
5197  // We mark all ghost cells by setting the user_flag and then request
5198  // these cells from the corresponding owners. As this information
5199  // can be incomplete,
5200  {
5201  std::vector<bool> user_flags;
5202  triangulation->save_user_flags(user_flags);
5203  triangulation->clear_user_flags();
5204 
5205  // mark all ghost cells for transfer
5206  {
5207  typename DoFHandlerType::level_cell_iterator cell,
5208  endc = dof_handler->end();
5209  for (cell = dof_handler->begin(); cell != endc; ++cell)
5210  if (cell->level_subdomain_id() !=
5211  ::numbers::artificial_subdomain_id &&
5212  !cell->is_locally_owned_on_level())
5213  cell->set_user_flag();
5214  }
5215 
5216  // Phase 1. Request all marked cells from corresponding owners. If we
5217  // managed to get every DoF, remove the user_flag, otherwise we
5218  // will request them again in the step below.
5219  communicate_mg_ghost_cells(*triangulation, *dof_handler);
5220 
5221  // have a barrier so that sends from above and below this
5222  // place are not mixed up.
5223  //
5224  // this is necessary because above we just see if there are
5225  // messages and then receive them, without discriminating
5226  // where they come from and whether they were sent in phase
5227  // 1 or 2 in communicate_mg_ghost_cells() on another
5228  // processor. the need for a global communication step like
5229  // this barrier could be avoided by receiving messages
5230  // specifically from those processors from which we expect
5231  // messages, and by using different tags for phase 1 and 2,
5232  // but the cost of a barrier is negligible compared to
5233  // everything else we do here
5234  const int ierr = MPI_Barrier(triangulation->get_communicator());
5235  AssertThrowMPI(ierr);
5236 
5237  // Phase 2, only request the cells that were not completed
5238  // in Phase 1.
5239  communicate_mg_ghost_cells(*triangulation, *dof_handler);
5240 
5241 # ifdef DEBUG
5242  // make sure we have removed all flags:
5243  {
5244  typename DoFHandlerType::level_cell_iterator cell,
5245  endc = dof_handler->end();
5246  for (cell = dof_handler->begin(); cell != endc; ++cell)
5247  if (cell->level_subdomain_id() !=
5248  ::numbers::artificial_subdomain_id &&
5249  !cell->is_locally_owned_on_level())
5250  Assert(cell->user_flag_set() == false, ExcInternalError());
5251  }
5252 # endif
5253 
5254  triangulation->load_user_flags(user_flags);
5255  }
5256 
5257 
5258 
5259 # ifdef DEBUG
5260  // check that we are really done
5261  {
5262  std::vector<::types::global_dof_index> local_dof_indices;
5263  typename DoFHandlerType::level_cell_iterator cell,
5264  endc = dof_handler->end();
5265 
5266  for (cell = dof_handler->begin(); cell != endc; ++cell)
5267  if (cell->level_subdomain_id() !=
5268  ::numbers::artificial_subdomain_id)
5269  {
5270  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
5271  cell->get_mg_dof_indices(local_dof_indices);
5272  if (local_dof_indices.end() !=
5273  std::find(local_dof_indices.begin(),
5274  local_dof_indices.end(),
5276  {
5277  Assert(false, ExcMessage("not all DoFs got distributed!"));
5278  }
5279  }
5280  }
5281 # endif // DEBUG
5282 
5283  return number_caches;
5284 
5285 #endif // DEAL_II_WITH_P4EST
5286  }
5287 
5288 
5289  template <class DoFHandlerType>
5290  NumberCache
5292  const std::vector<::types::global_dof_index> &new_numbers) const
5293  {
5294  (void)new_numbers;
5295 
5296  Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
5297  ExcInternalError());
5298 
5299 #ifndef DEAL_II_WITH_P4EST
5300  Assert(false, ExcNotImplemented());
5301  return NumberCache();
5302 #else
5303  const unsigned int dim = DoFHandlerType::dimension;
5304  const unsigned int spacedim = DoFHandlerType::space_dimension;
5305 
5308  *>(const_cast<::Triangulation<dim, spacedim> *>(
5309  &dof_handler->get_triangulation())));
5310  Assert(triangulation != nullptr, ExcInternalError());
5311 
5312 
5313  // We start by checking whether only the numbering within the MPI
5314  // ranks changed. In that case, we can apply the renumbering with some
5315  // local renumbering only (this is similar to the renumber_mg_dofs()
5316  // function below)
5317  bool locally_owned_set_changes = false;
5318  for (types::global_dof_index i : new_numbers)
5319  if (dof_handler->locally_owned_dofs().is_element(i) == false)
5320  {
5321  locally_owned_set_changes = true;
5322  break;
5323  }
5324 
5325  if (Utilities::MPI::sum(static_cast<unsigned int>(
5326  locally_owned_set_changes),
5327  triangulation->get_communicator()) == 0)
5328  {
5329  // Since only the order within the local subdomains has changed,
5330  // all we need to do is to propagate the knowledge about the
5331  // numbers from the locally owned dofs (given by the new_numbers
5332  // array) to all ghosted dofs on neighboring processors. We can do
5333  // this by ghost layer exchange routines as in parallel vectors:
5334  // We create an IndexSet for the relevant dofs and then export
5335  // into an array of those values via Utilities::MPI::Partitioner.
5336  IndexSet relevant_dofs;
5338  relevant_dofs);
5339  std::vector<types::global_dof_index> ghosted_new_numbers(
5340  relevant_dofs.n_elements());
5341  {
5342  Utilities::MPI::Partitioner partitioner(
5343  dof_handler->locally_owned_dofs(),
5344  relevant_dofs,
5345  triangulation->get_communicator());
5346 
5347  // choose some number that makes it unlikely to get conflicts
5348  // with other ongoing non-blocking communication (there
5349  // shouldn't be any at this place in most programs).
5350  const unsigned int communication_channel = 19;
5351  std::vector<types::global_dof_index> temp_array(
5352  partitioner.n_import_indices());
5353  std::vector<MPI_Request> requests;
5354  partitioner.export_to_ghosted_array_start(
5355  communication_channel,
5356  make_array_view(new_numbers),
5357  make_array_view(temp_array),
5359  ghosted_new_numbers.data() + new_numbers.size(),
5360  partitioner.n_ghost_indices()),
5361  requests);
5362  partitioner.export_to_ghosted_array_finish(
5364  ghosted_new_numbers.data() + new_numbers.size(),
5365  partitioner.n_ghost_indices()),
5366  requests);
5367 
5368  // we need to fill the indices of the locally owned part into
5369  // the new numbers array, which is not provided by the parallel
5370  // partitioner. their right position is somewhere in the middle
5371  // of the array, so we first copy the ghosted part from smaller
5372  // ranks to the front, then insert the data in the middle.
5373  unsigned int n_ghosts_on_smaller_ranks = 0;
5374  for (std::pair<unsigned int, unsigned int> t :
5375  partitioner.ghost_targets())
5376  {
5377  if (t.first > partitioner.this_mpi_process())
5378  break;
5379  n_ghosts_on_smaller_ranks += t.second;
5380  }
5381  if (n_ghosts_on_smaller_ranks > 0)
5382  {
5383  Assert(ghosted_new_numbers.data() != nullptr,
5384  ExcInternalError());
5385  std::memmove(ghosted_new_numbers.data(),
5386  ghosted_new_numbers.data() + new_numbers.size(),
5387  sizeof(types::global_dof_index) *
5388  n_ghosts_on_smaller_ranks);
5389  }
5390  if (new_numbers.size() > 0)
5391  {
5392  Assert(new_numbers.data() != nullptr, ExcInternalError());
5393  std::memcpy(ghosted_new_numbers.data() +
5394  n_ghosts_on_smaller_ranks,
5395  new_numbers.data(),
5396  sizeof(types::global_dof_index) *
5397  new_numbers.size());
5398  }
5399  }
5400 
5401  // In case we do not carry any relevant dof (but only some remote
5402  // processor), we do not need to call the renumbering. We call the
5403  // version without validity check because vertex dofs will be
5404  // set already in the artificial region.
5405  if (relevant_dofs.n_elements() > 0)
5406  Implementation::renumber_dofs(ghosted_new_numbers,
5407  relevant_dofs,
5408  *dof_handler,
5409  /*check_validity=*/false);
5410 
5411  NumberCache number_cache;
5412  number_cache.locally_owned_dofs = dof_handler->locally_owned_dofs();
5413  number_cache.n_global_dofs = dof_handler->n_dofs();
5414  number_cache.n_locally_owned_dofs =
5415  number_cache.locally_owned_dofs.n_elements();
5416  return number_cache;
5417  }
5418  else
5419  {
5420  // Now back to the more complicated case
5421  //
5422  // First figure out the new set of locally owned DoF indices.
5423  // If we own no DoFs, we still need to go through this function,
5424  // but we can skip this calculation.
5425  //
5426  // The IndexSet::add_indices() function is substantially more
5427  // efficient if the set of indices is already sorted because
5428  // it can then insert ranges instead of individual elements.
5429  // consequently, pre-sort the array of new indices
5430  IndexSet my_locally_owned_new_dof_indices(dof_handler->n_dofs());
5431  if (dof_handler->n_locally_owned_dofs() > 0)
5432  {
5433  std::vector<::types::global_dof_index>
5434  new_numbers_sorted = new_numbers;
5435  std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
5436 
5437  my_locally_owned_new_dof_indices.add_indices(
5438  new_numbers_sorted.begin(), new_numbers_sorted.end());
5439  my_locally_owned_new_dof_indices.compress();
5440 
5441  Assert(my_locally_owned_new_dof_indices.n_elements() ==
5442  new_numbers.size(),
5443  ExcInternalError());
5444  }
5445 
5446  // delete all knowledge of DoF indices that are not locally
5447  // owned. we do so by getting DoF indices on cells, checking
5448  // whether they are locally owned, if not, setting them to
5449  // an invalid value, and then setting them again on the current
5450  // cell
5451  //
5452  // DoFs we (i) know about, and (ii) don't own locally must be
5453  // located either on ghost cells, or on the interface between a
5454  // locally owned cell and a ghost cell. In any case, it is
5455  // sufficient to kill them only from the ghost side cell, so loop
5456  // only over ghost cells
5457  {
5458  std::vector<::types::global_dof_index> local_dof_indices;
5459 
5460  for (auto cell : dof_handler->active_cell_iterators())
5461  if (cell->is_ghost())
5462  {
5463  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
5464  cell->get_dof_indices(local_dof_indices);
5465 
5466  for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
5467  ++i)
5468  // delete a DoF index if it has not already been deleted
5469  // (e.g., by visiting a neighboring cell, if it is on the
5470  // boundary), and if we don't own it
5471  if ((local_dof_indices[i] !=
5473  (!dof_handler->locally_owned_dofs().is_element(
5474  local_dof_indices[i])))
5475  local_dof_indices[i] = numbers::invalid_dof_index;
5476 
5477  cell->set_dof_indices(local_dof_indices);
5478  }
5479  }
5480 
5481 
5482  // renumber. Skip when there is nothing to do because we own no DoF.
5483  if (dof_handler->locally_owned_dofs().n_elements() > 0)
5484  Implementation::renumber_dofs(new_numbers,
5485  dof_handler->locally_owned_dofs(),
5486  *dof_handler,
5487  /*check_validity=*/false);
5488 
5489  // Communicate newly assigned DoF indices to other processors
5490  // and get the same information for our own ghost cells.
5491  //
5492  // This is the same as phase 5+6 in the distribute_dofs() algorithm,
5493  // taking into account that we have to unify a few DoFs in between
5494  // then communication phases if we do hp numbering
5495  {
5496  std::vector<bool> user_flags;
5497  triangulation->save_user_flags(user_flags);
5498  triangulation->clear_user_flags();
5499 
5500  // mark all own cells for transfer
5501  for (const auto &cell : dof_handler->active_cell_iterators())
5502  if (!cell->is_artificial())
5503  cell->set_user_flag();
5504 
5505  // figure out which cells are ghost cells on which we have
5506  // to exchange DoF indices
5507  const std::map<unsigned int,
5508  std::set<::types::subdomain_id>>
5509  vertices_with_ghost_neighbors =
5510  triangulation->compute_vertices_with_ghost_neighbors();
5511 
5512 
5513  // Send and receive cells. After this, only the local cells
5514  // are marked, that received new data. This has to be
5515  // communicated in a second communication step.
5516  //
5517  // as explained in the 'distributed' paper, this has to be
5518  // done twice
5519  communicate_dof_indices_on_marked_cells(
5520  *dof_handler, vertices_with_ghost_neighbors);
5521 
5522  // in case of hp::DoFHandlers, we may have received valid
5523  // indices of degrees of freedom that are dominated by a fe
5524  // object adjacent to a ghost interface.
5525  // thus, we overwrite the remaining invalid indices with
5526  // the valid ones in this step.
5527  Implementation::merge_invalid_dof_indices_on_ghost_interfaces(
5528  *dof_handler);
5529 
5530  communicate_dof_indices_on_marked_cells(
5531  *dof_handler, vertices_with_ghost_neighbors);
5532 
5533  triangulation->load_user_flags(user_flags);
5534  }
5535 
5536  NumberCache number_cache;
5537  number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
5538  number_cache.n_global_dofs = dof_handler->n_dofs();
5539  number_cache.n_locally_owned_dofs =
5540  number_cache.locally_owned_dofs.n_elements();
5541  return number_cache;
5542  }
5543 #endif
5544  }
5545 
5546 
5547 
5548  template <class DoFHandlerType>
5549  NumberCache
5551  const unsigned int level,
5552  const std::vector<types::global_dof_index> &new_numbers) const
5553  {
5554  // we only implement the case where the multigrid numbers are
5555  // renumbered within the processor's partition, rather than the most
5556  // general case
5557  const IndexSet index_set = dof_handler->locally_owned_mg_dofs(level);
5558 
5559 #ifdef DEAL_II_WITH_MPI
5560 
5561  constexpr int dim = DoFHandlerType::dimension;
5562  constexpr int spacedim = DoFHandlerType::space_dimension;
5564  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
5565  &this->dof_handler->get_triangulation()));
5566  Assert(tr != nullptr, ExcInternalError());
5567 
5568  const unsigned int my_rank =
5570 
5571 # ifdef DEBUG
5572  for (types::global_dof_index i : new_numbers)
5573  {
5574  Assert(index_set.is_element(i),
5576  "Renumberings that change the locally owned mg dofs "
5577  "partitioning are currently not implemented for "
5578  "the multigrid levels"));
5579  }
5580 # endif
5581 
5582  // we need to access all locally relevant degrees of freedom. we
5583  // use Utilities::MPI::Partitioner for handling the data exchange
5584  // of the new numbers, which is simply the extraction of ghost data
5585  IndexSet relevant_dofs;
5587  level,
5588  relevant_dofs);
5589  std::vector<types::global_dof_index> ghosted_new_numbers(
5590  relevant_dofs.n_elements());
5591  {
5592  Utilities::MPI::Partitioner partitioner(index_set,
5593  relevant_dofs,
5594  tr->get_communicator());
5595  std::vector<types::global_dof_index> temp_array(
5596  partitioner.n_import_indices());
5597  const unsigned int communication_channel = 17;
5598  std::vector<MPI_Request> requests;
5599  partitioner.export_to_ghosted_array_start(
5600  communication_channel,
5601  make_array_view(new_numbers),
5602  make_array_view(temp_array),
5603  ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
5604  new_numbers.size(),
5605  partitioner.n_ghost_indices()),
5606  requests);
5607  partitioner.export_to_ghosted_array_finish(
5608  ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
5609  new_numbers.size(),
5610  partitioner.n_ghost_indices()),
5611  requests);
5612 
5613  // we need to fill the indices of the locally owned part into the
5614  // new numbers array. their right position is somewhere in the
5615  // middle of the array, so we first copy the ghosted part from
5616  // smaller ranks to the front, then insert the data in the middle.
5617  unsigned int n_ghosts_on_smaller_ranks = 0;
5618  for (std::pair<unsigned int, unsigned int> t :
5619  partitioner.ghost_targets())
5620  {
5621  if (t.first > my_rank)
5622  break;
5623  n_ghosts_on_smaller_ranks += t.second;
5624  }
5625  if (n_ghosts_on_smaller_ranks > 0)
5626  {
5627  Assert(ghosted_new_numbers.data() != nullptr, ExcInternalError());
5628  std::memmove(ghosted_new_numbers.data(),
5629  ghosted_new_numbers.data() + new_numbers.size(),
5630  sizeof(types::global_dof_index) *
5631  n_ghosts_on_smaller_ranks);
5632  }
5633  if (new_numbers.size() > 0)
5634  {
5635  Assert(new_numbers.data() != nullptr, ExcInternalError());
5636  std::memcpy(ghosted_new_numbers.data() +
5637  n_ghosts_on_smaller_ranks,
5638  new_numbers.data(),
5639  sizeof(types::global_dof_index) * new_numbers.size());
5640  }
5641  }
5642 
5643  // in case we do not own any of the given level (but only some remote
5644  // processor), we do not need to call the renumbering
5645  if (level < this->dof_handler->get_triangulation().n_levels() &&
5646  relevant_dofs.n_elements() > 0)
5647  Implementation::renumber_mg_dofs(
5648  ghosted_new_numbers, relevant_dofs, *dof_handler, level, true);
5649 #else
5650  (void)new_numbers;
5651  Assert(false, ExcNotImplemented());
5652 #endif
5653 
5654  NumberCache number_cache;
5655  number_cache.locally_owned_dofs = index_set;
5656  number_cache.n_global_dofs = dof_handler->n_dofs();
5657  number_cache.n_locally_owned_dofs =
5658  number_cache.locally_owned_dofs.n_elements();
5659  return number_cache;
5660  }
5661  } // namespace Policy
5662  } // namespace DoFHandlerImplementation
5663 } // namespace internal
5664 
5665 
5666 
5667 /*-------------- Explicit Instantiations -------------------------------*/
5668 #include "dof_handler_policy.inst"
5669 
5670 
5671 DEAL_II_NAMESPACE_CLOSE
unsigned int n_active_cells() const
Definition: tria.cc:12601
void get_active_fe_indices(std::vector< unsigned int > &active_fe_indices) const
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1402
const Triangulation< dim, spacedim > & get_triangulation() const
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
const std::set< types::subdomain_id > & level_ghost_owners() const
Definition: tria_base.cc:350
static const unsigned int invalid_unsigned_int
Definition: types.h:187
virtual bool is_multilevel_hierarchy_constructed() const =0
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandlerType &dof_handler)
Definition: dof_tools.cc:1272
const types::subdomain_id invalid_subdomain_id
Definition: types.h:279
void clear_user_flags()
Definition: tria.cc:11100
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:931
typename ActiveSelector::quad_iterator quad_iterator
Definition: dof_handler.h:271
Triangulation< dim, spacedim >::cell_iterator to_cell(const Triangulation< dim, spacedim > &tria) const
Definition: cell_id.cc:159
unsigned int n_cells() const
Definition: tria.cc:12593
size_type nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:1847
Task< RT > new_task(const std::function< RT()> &function)
const std::vector< types::subdomain_id > & get_true_subdomain_ids_of_cells() const
Definition: shared_tria.cc:324
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12111
void load_user_flags(std::istream &in)
Definition: tria.cc:11160
const hp::FECollection< dim, spacedim > & get_fe() const
cell_iterator end() const
Definition: dof_handler.cc:960
const unsigned int dofs_per_quad
Definition: fe_base.h:237
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:950
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1670
active_cell_iterator begin_active(const unsigned int level=0) const
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11939
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
unsigned int coarse_cell_id
Definition: types.h:109
const unsigned int dofs_per_line
Definition: fe_base.h:231
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
const std::vector< types::subdomain_id > & get_true_level_subdomain_ids_of_cells(const unsigned int level) const
Definition: shared_tria.cc:333
void extract_locally_relevant_level_dofs(const DoFHandlerType &dof_handler, const unsigned int level, IndexSet &dof_set)
Definition: dof_tools.cc:1120
virtual std::vector< NumberCache > distribute_mg_dofs() const override
size_type size() const
Definition: index_set.h:1600
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:312
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:944
unsigned int n_locally_owned_dofs() const
virtual std::vector< NumberCache > distribute_mg_dofs() const override
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::quad_iterator quad_iterator
Definition: dof_handler.h:270
virtual std::vector< NumberCache > distribute_mg_dofs() const override
unsigned int subdomain_id
Definition: types.h:43
T sum(const T &t, const MPI_Comm &mpi_communicator)
const unsigned int dofs_per_hex
Definition: fe_base.h:243
#define Assert(cond, exc)
Definition: exceptions.h:1407
IteratorRange< active_cell_iterator > active_cell_iterators() const
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1888
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
unsigned int max_dofs_per_cell(const DoFHandler< dim, spacedim > &dh)
types::global_dof_index n_dofs() const
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1078
types::subdomain_id locally_owned_subdomain() const override
Definition: tria_base.cc:332
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:928
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > faces
Definition: dof_handler.h:1422
virtual MPI_Comm get_communicator() const
Definition: tria_base.cc:173
const hp::FECollection< dim, spacedim > & get_fe_collection() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Definition: hp.h:117
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > levels
Definition: dof_handler.h:1410
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:74
Definition: cell_id.h:68
void add_range(const size_type begin, const size_type end)
Definition: index_set.h:1640
virtual std::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors() const
Definition: tria_base.cc:359
unsigned int global_dof_index
Definition: types.h:89
const types::subdomain_id artificial_subdomain_id
Definition: types.h:296
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
void compress() const
Definition: index_set.h:1608
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const
Definition: fe.cc:939
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1695
const std::vector< IndexSet > & locally_owned_dofs_per_processor() const
typename ActiveSelector::line_iterator line_iterator
Definition: dof_handler.h:247
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1122
virtual unsigned int n_global_levels() const override
Definition: tria_base.cc:130
const Triangulation< dim, spacedim > & get_triangulation() const
void export_to_ghosted_array_start(const unsigned int communication_channel, const ArrayView< const Number, MemorySpaceType > &locally_owned_array, const ArrayView< Number, MemorySpaceType > &temporary_storage, const ArrayView< Number, MemorySpaceType > &ghost_array, std::vector< MPI_Request > &requests) const
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:85
void save_user_flags(std::ostream &out) const
Definition: tria.cc:11111
static ::ExceptionBase & ExcNotImplemented()
Iterator points to a valid object.
const unsigned int dofs_per_vertex
Definition: fe_base.h:225
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:1167
static unsigned int n_threads()
Definition: table.h:699
bool is_element(const size_type index) const
Definition: index_set.h:1732
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:356
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:328
const types::global_dof_index invalid_dof_index
Definition: types.h:202
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 IndexSet & locally_owned_dofs() const
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
size_type n_elements() const
Definition: index_set.h:1799
cell_iterator end() const
std::vector< types::global_dof_index > vertex_dofs
Definition: dof_handler.h:1396
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
typename ActiveSelector::line_iterator line_iterator
Definition: dof_handler.h:246
static ::ExceptionBase & ExcInternalError()