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