Reference documentation for deal.II version GIT 276cfaa7ee 2022-08-16 14:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_handler.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2021 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 #include <deal.II/base/config.h>
17 
20 #include <deal.II/base/mpi.templates.h>
22 
23 #include <deal.II/distributed/cell_data_transfer.templates.h>
27 
30 
33 #include <deal.II/grid/tria.h>
37 
38 #include <algorithm>
39 #include <memory>
40 #include <set>
41 #include <unordered_set>
42 
44 
45 template <int dim, int spacedim>
47 
48 
49 template <int dim, int spacedim>
52 
53 
54 namespace internal
55 {
56  template <int dim, int spacedim>
57  std::string
58  policy_to_string(const ::internal::DoFHandlerImplementation::Policy::
59  PolicyBase<dim, spacedim> &policy)
60  {
61  std::string policy_name;
62  if (dynamic_cast<const typename ::internal::DoFHandlerImplementation::
63  Policy::Sequential<dim, spacedim> *>(&policy) ||
64  dynamic_cast<const typename ::internal::DoFHandlerImplementation::
65  Policy::Sequential<dim, spacedim> *>(&policy))
66  policy_name = "Policy::Sequential<";
67  else if (dynamic_cast<
68  const typename ::internal::DoFHandlerImplementation::
69  Policy::ParallelDistributed<dim, spacedim> *>(&policy) ||
70  dynamic_cast<
71  const typename ::internal::DoFHandlerImplementation::
72  Policy::ParallelDistributed<dim, spacedim> *>(&policy))
73  policy_name = "Policy::ParallelDistributed<";
74  else if (dynamic_cast<
75  const typename ::internal::DoFHandlerImplementation::
76  Policy::ParallelShared<dim, spacedim> *>(&policy) ||
77  dynamic_cast<
78  const typename ::internal::DoFHandlerImplementation::
79  Policy::ParallelShared<dim, spacedim> *>(&policy))
80  policy_name = "Policy::ParallelShared<";
81  else
83  policy_name += Utilities::int_to_string(dim) + "," +
84  Utilities::int_to_string(spacedim) + ">";
85  return policy_name;
86  }
87 
88 
89  namespace DoFHandlerImplementation
90  {
96  {
101  template <int spacedim>
102  static unsigned int
104  {
105  return std::min(static_cast<types::global_dof_index>(
106  3 * dof_handler.fe_collection.max_dofs_per_vertex() +
107  2 * dof_handler.fe_collection.max_dofs_per_line()),
108  dof_handler.n_dofs());
109  }
110 
111  template <int spacedim>
112  static unsigned int
114  {
115  // get these numbers by drawing pictures
116  // and counting...
117  // example:
118  // | | |
119  // --x-----x--x--X--
120  // | | | |
121  // | x--x--x
122  // | | | |
123  // --x--x--*--x--x--
124  // | | | |
125  // x--x--x |
126  // | | | |
127  // --X--x--x-----x--
128  // | | |
129  // x = vertices connected with center vertex *;
130  // = total of 19
131  // (the X vertices are connected with * if
132  // the vertices adjacent to X are hanging
133  // nodes)
134  // count lines -> 28 (don't forget to count
135  // mother and children separately!)
136  types::global_dof_index max_couplings;
137  switch (dof_handler.tria->max_adjacent_cells())
138  {
139  case 4:
140  max_couplings =
141  19 * dof_handler.fe_collection.max_dofs_per_vertex() +
142  28 * dof_handler.fe_collection.max_dofs_per_line() +
143  8 * dof_handler.fe_collection.max_dofs_per_quad();
144  break;
145  case 5:
146  max_couplings =
147  21 * dof_handler.fe_collection.max_dofs_per_vertex() +
148  31 * dof_handler.fe_collection.max_dofs_per_line() +
149  9 * dof_handler.fe_collection.max_dofs_per_quad();
150  break;
151  case 6:
152  max_couplings =
153  28 * dof_handler.fe_collection.max_dofs_per_vertex() +
154  42 * dof_handler.fe_collection.max_dofs_per_line() +
155  12 * dof_handler.fe_collection.max_dofs_per_quad();
156  break;
157  case 7:
158  max_couplings =
159  30 * dof_handler.fe_collection.max_dofs_per_vertex() +
160  45 * dof_handler.fe_collection.max_dofs_per_line() +
161  13 * dof_handler.fe_collection.max_dofs_per_quad();
162  break;
163  case 8:
164  max_couplings =
165  37 * dof_handler.fe_collection.max_dofs_per_vertex() +
166  56 * dof_handler.fe_collection.max_dofs_per_line() +
167  16 * dof_handler.fe_collection.max_dofs_per_quad();
168  break;
169 
170  // the following numbers are not based on actual counting but by
171  // extrapolating the number sequences from the previous ones (for
172  // example, for n_dofs_per_vertex(), the sequence above is 19, 21,
173  // 28, 30, 37, and is continued as follows):
174  case 9:
175  max_couplings =
176  39 * dof_handler.fe_collection.max_dofs_per_vertex() +
177  59 * dof_handler.fe_collection.max_dofs_per_line() +
178  17 * dof_handler.fe_collection.max_dofs_per_quad();
179  break;
180  case 10:
181  max_couplings =
182  46 * dof_handler.fe_collection.max_dofs_per_vertex() +
183  70 * dof_handler.fe_collection.max_dofs_per_line() +
184  20 * dof_handler.fe_collection.max_dofs_per_quad();
185  break;
186  case 11:
187  max_couplings =
188  48 * dof_handler.fe_collection.max_dofs_per_vertex() +
189  73 * dof_handler.fe_collection.max_dofs_per_line() +
190  21 * dof_handler.fe_collection.max_dofs_per_quad();
191  break;
192  case 12:
193  max_couplings =
194  55 * dof_handler.fe_collection.max_dofs_per_vertex() +
195  84 * dof_handler.fe_collection.max_dofs_per_line() +
196  24 * dof_handler.fe_collection.max_dofs_per_quad();
197  break;
198  case 13:
199  max_couplings =
200  57 * dof_handler.fe_collection.max_dofs_per_vertex() +
201  87 * dof_handler.fe_collection.max_dofs_per_line() +
202  25 * dof_handler.fe_collection.max_dofs_per_quad();
203  break;
204  case 14:
205  max_couplings =
206  63 * dof_handler.fe_collection.max_dofs_per_vertex() +
207  98 * dof_handler.fe_collection.max_dofs_per_line() +
208  28 * dof_handler.fe_collection.max_dofs_per_quad();
209  break;
210  case 15:
211  max_couplings =
212  65 * dof_handler.fe_collection.max_dofs_per_vertex() +
213  103 * dof_handler.fe_collection.max_dofs_per_line() +
214  29 * dof_handler.fe_collection.max_dofs_per_quad();
215  break;
216  case 16:
217  max_couplings =
218  72 * dof_handler.fe_collection.max_dofs_per_vertex() +
219  114 * dof_handler.fe_collection.max_dofs_per_line() +
220  32 * dof_handler.fe_collection.max_dofs_per_quad();
221  break;
222 
223  default:
224  Assert(false, ExcNotImplemented());
225  max_couplings = 0;
226  }
227  return std::min(max_couplings, dof_handler.n_dofs());
228  }
229 
230  template <int spacedim>
231  static unsigned int
233  {
234  // TODO:[?] Invent significantly better estimates than the ones in this
235  // function
236 
237  // doing the same thing here is a rather complicated thing, compared
238  // to the 2d case, since it is hard to draw pictures with several
239  // refined hexahedra :-) so I presently only give a coarse
240  // estimate for the case that at most 8 hexes meet at each vertex
241  //
242  // can anyone give better estimate here?
243  const unsigned int max_adjacent_cells =
244  dof_handler.tria->max_adjacent_cells();
245 
246  types::global_dof_index max_couplings;
247  if (max_adjacent_cells <= 8)
248  max_couplings =
249  7 * 7 * 7 * dof_handler.fe_collection.max_dofs_per_vertex() +
250  7 * 6 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_line() +
251  9 * 4 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_quad() +
252  27 * dof_handler.fe_collection.max_dofs_per_hex();
253  else
254  {
255  Assert(false, ExcNotImplemented());
256  max_couplings = 0;
257  }
258 
259  return std::min(max_couplings, dof_handler.n_dofs());
260  }
261 
266  template <int dim, int spacedim>
267  static void
269  {
270  dof_handler.object_dof_indices.clear();
271  dof_handler.object_dof_indices.resize(dof_handler.tria->n_levels());
272  dof_handler.object_dof_indices.shrink_to_fit();
273 
274  dof_handler.object_dof_ptr.clear();
275  dof_handler.object_dof_ptr.resize(dof_handler.tria->n_levels());
276  dof_handler.object_dof_ptr.shrink_to_fit();
277 
278  dof_handler.cell_dof_cache_indices.clear();
279  dof_handler.cell_dof_cache_indices.resize(dof_handler.tria->n_levels());
280  dof_handler.cell_dof_cache_indices.shrink_to_fit();
281 
282  dof_handler.cell_dof_cache_ptr.clear();
283  dof_handler.cell_dof_cache_ptr.resize(dof_handler.tria->n_levels());
284  dof_handler.cell_dof_cache_ptr.shrink_to_fit();
285  }
286 
290  template <int dim, int spacedim>
291  static void
293  const unsigned int n_inner_dofs_per_cell)
294  {
295  for (unsigned int i = 0; i < dof_handler.tria->n_levels(); ++i)
296  {
297  // 1) object_dof_indices
298  dof_handler.object_dof_ptr[i][dim].assign(
299  dof_handler.tria->n_raw_cells(i) + 1, 0);
300 
301  for (const auto &cell :
302  dof_handler.tria->cell_iterators_on_level(i))
303  if (cell->is_active() && !cell->is_artificial())
304  dof_handler.object_dof_ptr[i][dim][cell->index() + 1] =
305  n_inner_dofs_per_cell;
306 
307  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i); ++j)
308  dof_handler.object_dof_ptr[i][dim][j + 1] +=
309  dof_handler.object_dof_ptr[i][dim][j];
310 
311  dof_handler.object_dof_indices[i][dim].resize(
312  dof_handler.object_dof_ptr[i][dim].back(),
314 
315  // 2) cell_dof_cache_indices
316  dof_handler.cell_dof_cache_ptr[i].assign(
317  dof_handler.tria->n_raw_cells(i) + 1, 0);
318 
319  for (const auto &cell :
320  dof_handler.tria->cell_iterators_on_level(i))
321  if (cell->is_active() && !cell->is_artificial())
322  dof_handler.cell_dof_cache_ptr[i][cell->index() + 1] =
323  dof_handler.get_fe().n_dofs_per_cell();
324 
325  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i); ++j)
326  dof_handler.cell_dof_cache_ptr[i][j + 1] +=
327  dof_handler.cell_dof_cache_ptr[i][j];
328 
329  dof_handler.cell_dof_cache_indices[i].resize(
330  dof_handler.cell_dof_cache_ptr[i].back(),
332  }
333  }
334 
339  template <int dim, int spacedim, typename T>
340  static void
342  const unsigned int structdim,
343  const unsigned int n_raw_entities,
344  const T & cell_process)
345  {
346  if (dof_handler.tria->n_cells() == 0)
347  return;
348 
349  dof_handler.object_dof_ptr[0][structdim].assign(n_raw_entities + 1, -1);
350  // determine for each entity the number of dofs
351  for (const auto &cell : dof_handler.tria->cell_iterators())
352  if (cell->is_active() && !cell->is_artificial())
353  cell_process(
354  cell,
355  [&](const unsigned int n_dofs_per_entity,
356  const unsigned int index) {
357  auto &n_dofs_per_entity_target =
358  dof_handler.object_dof_ptr[0][structdim][index + 1];
359 
360  // make sure that either the entity has not been visited or
361  // the entity has the same number of dofs assigned
362  Assert((n_dofs_per_entity_target ==
363  static_cast<
365  -1) ||
366  n_dofs_per_entity_target == n_dofs_per_entity),
368 
369  n_dofs_per_entity_target = n_dofs_per_entity;
370  });
371 
372  // convert the absolute numbers to CRS
373  dof_handler.object_dof_ptr[0][structdim][0] = 0;
374  for (unsigned int i = 1; i < n_raw_entities + 1; ++i)
375  {
376  if (dof_handler.object_dof_ptr[0][structdim][i] ==
377  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
378  -1))
379  dof_handler.object_dof_ptr[0][structdim][i] =
380  dof_handler.object_dof_ptr[0][structdim][i - 1];
381  else
382  dof_handler.object_dof_ptr[0][structdim][i] +=
383  dof_handler.object_dof_ptr[0][structdim][i - 1];
384  }
385 
386  // allocate memory for indices
387  dof_handler.object_dof_indices[0][structdim].resize(
388  dof_handler.object_dof_ptr[0][structdim].back(),
390  }
391 
398  template <int dim, int spacedim>
399  static void
401  {
402  reset_to_empty_objects(dof_handler);
403 
404  const auto &fe = dof_handler.get_fe();
405 
406  // cell
407  reserve_cells(dof_handler,
408  dim == 1 ? fe.n_dofs_per_line() :
409  (dim == 2 ? fe.n_dofs_per_quad(0) :
410  fe.n_dofs_per_hex()));
411 
412  // vertices
413  reserve_subentities(dof_handler,
414  0,
415  dof_handler.tria->n_vertices(),
416  [&](const auto &cell, const auto &process) {
417  for (const auto vertex_index :
418  cell->vertex_indices())
419  process(fe.n_dofs_per_vertex(),
420  cell->vertex_index(vertex_index));
421  });
422 
423  // lines
424  if (dim == 2 || dim == 3)
425  reserve_subentities(dof_handler,
426  1,
427  dof_handler.tria->n_raw_lines(),
428  [&](const auto &cell, const auto &process) {
429  for (const auto line_index :
430  cell->line_indices())
431  process(fe.n_dofs_per_line(),
432  cell->line(line_index)->index());
433  });
434 
435  // quads
436  if (dim == 3)
437  reserve_subentities(dof_handler,
438  2,
439  dof_handler.tria->n_raw_quads(),
440  [&](const auto &cell, const auto &process) {
441  for (const auto face_index :
442  cell->face_indices())
443  process(fe.n_dofs_per_quad(face_index),
444  cell->face(face_index)->index());
445  });
446  }
447 
448  template <int spacedim>
449  static void
451  {
452  Assert(dof_handler.get_triangulation().n_levels() > 0,
453  ExcMessage("Invalid triangulation"));
454  dof_handler.clear_mg_space();
455 
456  const ::Triangulation<1, spacedim> &tria =
457  dof_handler.get_triangulation();
458  const unsigned int dofs_per_line =
459  dof_handler.get_fe().n_dofs_per_line();
460  const unsigned int n_levels = tria.n_levels();
461 
462  for (unsigned int i = 0; i < n_levels; ++i)
463  {
464  dof_handler.mg_levels.emplace_back(
466  dof_handler.mg_levels.back()->dof_object.dofs =
467  std::vector<types::global_dof_index>(tria.n_raw_lines(i) *
468  dofs_per_line,
470  }
471 
472  const unsigned int n_vertices = tria.n_vertices();
473 
474  dof_handler.mg_vertex_dofs.resize(n_vertices);
475 
476  std::vector<unsigned int> max_level(n_vertices, 0);
477  std::vector<unsigned int> min_level(n_vertices, n_levels);
478 
479  for (typename ::Triangulation<1, spacedim>::cell_iterator cell =
480  tria.begin();
481  cell != tria.end();
482  ++cell)
483  {
484  const unsigned int level = cell->level();
485 
486  for (const auto vertex : cell->vertex_indices())
487  {
488  const unsigned int vertex_index = cell->vertex_index(vertex);
489 
490  if (min_level[vertex_index] > level)
491  min_level[vertex_index] = level;
492 
493  if (max_level[vertex_index] < level)
494  max_level[vertex_index] = level;
495  }
496  }
497 
498  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
499  if (tria.vertex_used(vertex))
500  {
501  Assert(min_level[vertex] < n_levels, ExcInternalError());
502  Assert(max_level[vertex] >= min_level[vertex],
503  ExcInternalError());
504  dof_handler.mg_vertex_dofs[vertex].init(
505  min_level[vertex],
506  max_level[vertex],
507  dof_handler.get_fe().n_dofs_per_vertex());
508  }
509 
510  else
511  {
512  Assert(min_level[vertex] == n_levels, ExcInternalError());
513  Assert(max_level[vertex] == 0, ExcInternalError());
514  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
515  }
516  }
517 
518  template <int spacedim>
519  static void
521  {
522  Assert(dof_handler.get_triangulation().n_levels() > 0,
523  ExcMessage("Invalid triangulation"));
524  dof_handler.clear_mg_space();
525 
526  const ::FiniteElement<2, spacedim> &fe = dof_handler.get_fe();
527  const ::Triangulation<2, spacedim> &tria =
528  dof_handler.get_triangulation();
529  const unsigned int n_levels = tria.n_levels();
530 
531  for (unsigned int i = 0; i < n_levels; ++i)
532  {
533  dof_handler.mg_levels.emplace_back(
534  std::make_unique<
536  dof_handler.mg_levels.back()->dof_object.dofs =
537  std::vector<types::global_dof_index>(
538  tria.n_raw_quads(i) *
539  fe.n_dofs_per_quad(0 /*note: in 2D there is only one quad*/),
541  }
542 
543  dof_handler.mg_faces =
544  std::make_unique<internal::DoFHandlerImplementation::DoFFaces<2>>();
545  dof_handler.mg_faces->lines.dofs =
546  std::vector<types::global_dof_index>(tria.n_raw_lines() *
547  fe.n_dofs_per_line(),
549 
550  const unsigned int n_vertices = tria.n_vertices();
551 
552  dof_handler.mg_vertex_dofs.resize(n_vertices);
553 
554  std::vector<unsigned int> max_level(n_vertices, 0);
555  std::vector<unsigned int> min_level(n_vertices, n_levels);
556 
557  for (typename ::Triangulation<2, spacedim>::cell_iterator cell =
558  tria.begin();
559  cell != tria.end();
560  ++cell)
561  {
562  const unsigned int level = cell->level();
563 
564  for (const auto vertex : cell->vertex_indices())
565  {
566  const unsigned int vertex_index = cell->vertex_index(vertex);
567 
568  if (min_level[vertex_index] > level)
569  min_level[vertex_index] = level;
570 
571  if (max_level[vertex_index] < level)
572  max_level[vertex_index] = level;
573  }
574  }
575 
576  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
577  if (tria.vertex_used(vertex))
578  {
579  Assert(min_level[vertex] < n_levels, ExcInternalError());
580  Assert(max_level[vertex] >= min_level[vertex],
581  ExcInternalError());
582  dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
583  max_level[vertex],
584  fe.n_dofs_per_vertex());
585  }
586 
587  else
588  {
589  Assert(min_level[vertex] == n_levels, ExcInternalError());
590  Assert(max_level[vertex] == 0, ExcInternalError());
591  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
592  }
593  }
594 
595  template <int spacedim>
596  static void
598  {
599  Assert(dof_handler.get_triangulation().n_levels() > 0,
600  ExcMessage("Invalid triangulation"));
601  dof_handler.clear_mg_space();
602 
603  const ::FiniteElement<3, spacedim> &fe = dof_handler.get_fe();
604  const ::Triangulation<3, spacedim> &tria =
605  dof_handler.get_triangulation();
606  const unsigned int n_levels = tria.n_levels();
607 
608  for (unsigned int i = 0; i < n_levels; ++i)
609  {
610  dof_handler.mg_levels.emplace_back(
611  std::make_unique<
613  dof_handler.mg_levels.back()->dof_object.dofs =
614  std::vector<types::global_dof_index>(tria.n_raw_hexs(i) *
615  fe.n_dofs_per_hex(),
617  }
618 
619  dof_handler.mg_faces =
620  std::make_unique<internal::DoFHandlerImplementation::DoFFaces<3>>();
621  dof_handler.mg_faces->lines.dofs =
622  std::vector<types::global_dof_index>(tria.n_raw_lines() *
623  fe.n_dofs_per_line(),
625 
626  // TODO: the implementation makes the assumption that all faces have the
627  // same number of dofs
628  AssertDimension(fe.n_unique_faces(), 1);
629  dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index>(
630  tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/),
632 
633  const unsigned int n_vertices = tria.n_vertices();
634 
635  dof_handler.mg_vertex_dofs.resize(n_vertices);
636 
637  std::vector<unsigned int> max_level(n_vertices, 0);
638  std::vector<unsigned int> min_level(n_vertices, n_levels);
639 
640  for (typename ::Triangulation<3, spacedim>::cell_iterator cell =
641  tria.begin();
642  cell != tria.end();
643  ++cell)
644  {
645  const unsigned int level = cell->level();
646 
647  for (const auto vertex : cell->vertex_indices())
648  {
649  const unsigned int vertex_index = cell->vertex_index(vertex);
650 
651  if (min_level[vertex_index] > level)
652  min_level[vertex_index] = level;
653 
654  if (max_level[vertex_index] < level)
655  max_level[vertex_index] = level;
656  }
657  }
658 
659  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
660  if (tria.vertex_used(vertex))
661  {
662  Assert(min_level[vertex] < n_levels, ExcInternalError());
663  Assert(max_level[vertex] >= min_level[vertex],
664  ExcInternalError());
665  dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
666  max_level[vertex],
667  fe.n_dofs_per_vertex());
668  }
669 
670  else
671  {
672  Assert(min_level[vertex] == n_levels, ExcInternalError());
673  Assert(max_level[vertex] == 0, ExcInternalError());
674  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
675  }
676  }
677  };
678  } // namespace DoFHandlerImplementation
679 
680 
681 
682  namespace hp
683  {
684  namespace DoFHandlerImplementation
685  {
691  {
697  template <int dim, int spacedim>
698  static void
700  DoFHandler<dim, spacedim> &dof_handler)
701  {
702  (void)dof_handler;
703  for (const auto &cell : dof_handler.active_cell_iterators())
704  if (cell->is_locally_owned())
705  Assert(
706  !cell->future_fe_index_set(),
707  ExcMessage(
708  "There shouldn't be any cells flagged for p-adaptation when partitioning."));
709  }
710 
711 
712 
717  template <int dim, int spacedim>
718  static void
720  {
721  // The final step in all of the reserve_space() functions is to set
722  // up vertex dof information. since vertices are sequentially
723  // numbered, what we do first is to set up an array in which
724  // we record whether a vertex is associated with any of the
725  // given fe's, by setting a bit. in a later step, we then
726  // actually allocate memory for the required dofs
727  //
728  // in the following, we only need to consider vertices that are
729  // adjacent to either a locally owned or a ghost cell; we never
730  // store anything on vertices that are only surrounded by
731  // artificial cells. so figure out that subset of vertices
732  // first
733  std::vector<bool> locally_used_vertices(
734  dof_handler.tria->n_vertices(), false);
735  for (const auto &cell : dof_handler.active_cell_iterators())
736  if (!cell->is_artificial())
737  for (const auto v : cell->vertex_indices())
738  locally_used_vertices[cell->vertex_index(v)] = true;
739 
740  std::vector<std::vector<bool>> vertex_fe_association(
741  dof_handler.fe_collection.size(),
742  std::vector<bool>(dof_handler.tria->n_vertices(), false));
743 
744  for (const auto &cell : dof_handler.active_cell_iterators())
745  if (!cell->is_artificial())
746  for (const auto v : cell->vertex_indices())
747  vertex_fe_association[cell->active_fe_index()]
748  [cell->vertex_index(v)] = true;
749 
750  // in debug mode, make sure that each vertex is associated
751  // with at least one FE (note that except for unused
752  // vertices, all vertices are actually active). this is of
753  // course only true for vertices that are part of either
754  // ghost or locally owned cells
755 #ifdef DEBUG
756  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
757  if (locally_used_vertices[v] == true)
758  if (dof_handler.tria->vertex_used(v) == true)
759  {
760  unsigned int fe = 0;
761  for (; fe < dof_handler.fe_collection.size(); ++fe)
762  if (vertex_fe_association[fe][v] == true)
763  break;
764  Assert(fe != dof_handler.fe_collection.size(),
765  ExcInternalError());
766  }
767 #endif
768 
769  const unsigned int d = 0;
770  const unsigned int l = 0;
771 
772  dof_handler.hp_object_fe_ptr[d].clear();
773  dof_handler.hp_object_fe_indices[d].clear();
774  dof_handler.object_dof_ptr[l][d].clear();
775  dof_handler.object_dof_indices[l][d].clear();
776 
777  dof_handler.hp_object_fe_ptr[d].reserve(
778  dof_handler.tria->n_vertices() + 1);
779 
780  unsigned int vertex_slots_needed = 0;
781  unsigned int fe_slots_needed = 0;
782 
783  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
784  {
785  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
786 
787  if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
788  {
789  for (unsigned int fe = 0;
790  fe < dof_handler.fe_collection.size();
791  ++fe)
792  if (vertex_fe_association[fe][v] == true)
793  {
794  fe_slots_needed++;
795  vertex_slots_needed +=
796  dof_handler.get_fe(fe).n_dofs_per_vertex();
797  }
798  }
799  }
800 
801  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
802 
803  dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
804  dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
805 
806  dof_handler.object_dof_indices[l][d].reserve(vertex_slots_needed);
807 
808  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
809  if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
810  {
811  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
812  ++fe)
813  if (vertex_fe_association[fe][v] == true)
814  {
815  dof_handler.hp_object_fe_indices[d].push_back(fe);
816  dof_handler.object_dof_ptr[l][d].push_back(
817  dof_handler.object_dof_indices[l][d].size());
818 
819  for (unsigned int i = 0;
820  i < dof_handler.get_fe(fe).n_dofs_per_vertex();
821  i++)
822  dof_handler.object_dof_indices[l][d].push_back(
824  }
825  }
826 
827 
828  dof_handler.object_dof_ptr[l][d].push_back(
829  dof_handler.object_dof_indices[l][d].size());
830 
831  AssertDimension(vertex_slots_needed,
832  dof_handler.object_dof_indices[l][d].size());
833  AssertDimension(fe_slots_needed,
834  dof_handler.hp_object_fe_indices[d].size());
835  AssertDimension(fe_slots_needed + 1,
836  dof_handler.object_dof_ptr[l][d].size());
837  AssertDimension(dof_handler.tria->n_vertices() + 1,
838  dof_handler.hp_object_fe_ptr[d].size());
839 
840  dof_handler.object_dof_indices[l][d].assign(
841  vertex_slots_needed, numbers::invalid_dof_index);
842  }
843 
844 
845 
850  template <int dim, int spacedim>
851  static void
853  {
854  (void)dof_handler;
855  // count how much space we need on each level for the cell
856  // dofs and set the dof_*_offsets data. initially set the
857  // latter to an invalid index, and only later set it to
858  // something reasonable for active dof_handler.cells
859  //
860  // note that for dof_handler.cells, the situation is simpler
861  // than for other (lower dimensional) objects since exactly
862  // one finite element is used for it
863  for (unsigned int level = 0; level < dof_handler.tria->n_levels();
864  ++level)
865  {
866  dof_handler.object_dof_ptr[level][dim] =
867  std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
868  dof_handler.tria->n_raw_cells(level),
869  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
870  -1));
871  dof_handler.cell_dof_cache_ptr[level] =
872  std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
873  dof_handler.tria->n_raw_cells(level),
874  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
875  -1));
876 
877  types::global_dof_index next_free_dof = 0;
878  types::global_dof_index cache_size = 0;
879 
880  for (auto cell :
882  if (cell->is_active() && !cell->is_artificial())
883  {
884  dof_handler.object_dof_ptr[level][dim][cell->index()] =
885  next_free_dof;
886  next_free_dof +=
887  cell->get_fe().template n_dofs_per_object<dim>();
888 
889  dof_handler.cell_dof_cache_ptr[level][cell->index()] =
890  cache_size;
891  cache_size += cell->get_fe().n_dofs_per_cell();
892  }
893 
894  dof_handler.object_dof_indices[level][dim] =
895  std::vector<types::global_dof_index>(
896  next_free_dof, numbers::invalid_dof_index);
897  dof_handler.cell_dof_cache_indices[level] =
898  std::vector<types::global_dof_index>(
899  cache_size, numbers::invalid_dof_index);
900  }
901  }
902 
903 
904 
909  template <int dim, int spacedim>
910  static void
912  {
913  // FACE DOFS
914  //
915  // Count face dofs, then allocate as much space
916  // as we need and prime the linked list for faces (see the
917  // description in hp::DoFLevel) with the indices we will
918  // need. Note that our task is more complicated than for the
919  // cell case above since two adjacent cells may have different
920  // active FE indices, in which case we need to allocate
921  // *two* sets of face dofs for the same face. But they don't
922  // *have* to be different, and so we need to prepare for this
923  // as well.
924  //
925  // The way we do things is that we loop over all active cells (these
926  // are the only ones that have DoFs anyway) and all their faces. We
927  // note in the vector face_touched whether we have previously
928  // visited a face and if so skip it
929  {
930  std::vector<bool> face_touched(dim == 2 ?
931  dof_handler.tria->n_raw_lines() :
932  dof_handler.tria->n_raw_quads());
933 
934  const unsigned int d = dim - 1;
935  const unsigned int l = 0;
936 
937  dof_handler.hp_object_fe_ptr[d].clear();
938  dof_handler.hp_object_fe_indices[d].clear();
939  dof_handler.object_dof_ptr[l][d].clear();
940  dof_handler.object_dof_indices[l][d].clear();
941 
942  dof_handler.hp_object_fe_ptr[d].resize(
943  dof_handler.tria->n_raw_faces() + 1);
944 
945  // An array to hold how many slots (see the hp::DoFLevel
946  // class) we will have to store on each level
947  unsigned int n_face_slots = 0;
948 
949  for (const auto &cell : dof_handler.active_cell_iterators())
950  if (!cell->is_artificial())
951  for (const auto face : cell->face_indices())
952  if (!face_touched[cell->face(face)->index()])
953  {
954  unsigned int fe_slots_needed = 0;
955 
956  if (cell->at_boundary(face) ||
957  cell->face(face)->has_children() ||
958  cell->neighbor_is_coarser(face) ||
959  (!cell->at_boundary(face) &&
960  cell->neighbor(face)->is_artificial()) ||
961  (!cell->at_boundary(face) &&
962  !cell->neighbor(face)->is_artificial() &&
963  (cell->active_fe_index() ==
964  cell->neighbor(face)->active_fe_index())))
965  {
966  fe_slots_needed = 1;
967  n_face_slots +=
968  dof_handler.get_fe(cell->active_fe_index())
969  .template n_dofs_per_object<dim - 1>(face);
970  }
971  else
972  {
973  fe_slots_needed = 2;
974  n_face_slots +=
975  dof_handler.get_fe(cell->active_fe_index())
976  .template n_dofs_per_object<dim - 1>(face) +
977  dof_handler
978  .get_fe(cell->neighbor(face)->active_fe_index())
979  .template n_dofs_per_object<dim - 1>(
980  cell->neighbor_face_no(face));
981  }
982 
983  // mark this face as visited
984  face_touched[cell->face(face)->index()] = true;
985 
986  dof_handler
987  .hp_object_fe_ptr[d][cell->face(face)->index() + 1] =
988  fe_slots_needed;
989  }
990 
991  for (unsigned int i = 1; i < dof_handler.hp_object_fe_ptr[d].size();
992  i++)
993  dof_handler.hp_object_fe_ptr[d][i] +=
994  dof_handler.hp_object_fe_ptr[d][i - 1];
995 
996 
997  dof_handler.hp_object_fe_indices[d].resize(
998  dof_handler.hp_object_fe_ptr[d].back());
999  dof_handler.object_dof_ptr[l][d].resize(
1000  dof_handler.hp_object_fe_ptr[d].back() + 1);
1001 
1002  dof_handler.object_dof_indices[l][d].reserve(n_face_slots);
1003 
1004 
1005  // With the memory now allocated, loop over the
1006  // dof_handler cells again and prime the _offset values as
1007  // well as the fe_index fields
1008  face_touched = std::vector<bool>(face_touched.size());
1009 
1010  for (const auto &cell : dof_handler.active_cell_iterators())
1011  if (!cell->is_artificial())
1012  for (const auto face : cell->face_indices())
1013  if (!face_touched[cell->face(face)->index()])
1014  {
1015  // Same decision tree as before
1016  if (cell->at_boundary(face) ||
1017  cell->face(face)->has_children() ||
1018  cell->neighbor_is_coarser(face) ||
1019  (!cell->at_boundary(face) &&
1020  cell->neighbor(face)->is_artificial()) ||
1021  (!cell->at_boundary(face) &&
1022  !cell->neighbor(face)->is_artificial() &&
1023  (cell->active_fe_index() ==
1024  cell->neighbor(face)->active_fe_index())))
1025  {
1026  const unsigned int fe = cell->active_fe_index();
1027  const unsigned int n_dofs =
1028  dof_handler.get_fe(fe)
1029  .template n_dofs_per_object<dim - 1>(face);
1030  const unsigned int offset =
1031  dof_handler
1032  .hp_object_fe_ptr[d][cell->face(face)->index()];
1033 
1034  dof_handler.hp_object_fe_indices[d][offset] = fe;
1035  dof_handler.object_dof_ptr[l][d][offset + 1] = n_dofs;
1036 
1037  for (unsigned int i = 0; i < n_dofs; ++i)
1038  dof_handler.object_dof_indices[l][d].push_back(
1040  }
1041  else
1042  {
1043  unsigned int fe_1 = cell->active_fe_index();
1044  unsigned int face_no_1 = face;
1045  unsigned int fe_2 =
1046  cell->neighbor(face)->active_fe_index();
1047  unsigned int face_no_2 = cell->neighbor_face_no(face);
1048 
1049  if (fe_2 < fe_1)
1050  {
1051  std::swap(fe_1, fe_2);
1052  std::swap(face_no_1, face_no_2);
1053  }
1054 
1055  const unsigned int n_dofs_1 =
1056  dof_handler.get_fe(fe_1)
1057  .template n_dofs_per_object<dim - 1>(face_no_1);
1058 
1059  const unsigned int n_dofs_2 =
1060  dof_handler.get_fe(fe_2)
1061  .template n_dofs_per_object<dim - 1>(face_no_2);
1062 
1063  const unsigned int offset =
1064  dof_handler
1065  .hp_object_fe_ptr[d][cell->face(face)->index()];
1066 
1067  dof_handler.hp_object_fe_indices[d].push_back(
1068  cell->active_fe_index());
1069  dof_handler.object_dof_ptr[l][d].push_back(
1070  dof_handler.object_dof_indices[l][d].size());
1071 
1072  dof_handler.hp_object_fe_indices[d][offset + 0] =
1073  fe_1;
1074  dof_handler.hp_object_fe_indices[d][offset + 1] =
1075  fe_2;
1076  dof_handler.object_dof_ptr[l][d][offset + 1] =
1077  n_dofs_1;
1078  dof_handler.object_dof_ptr[l][d][offset + 2] =
1079  n_dofs_2;
1080 
1081 
1082  for (unsigned int i = 0; i < n_dofs_1 + n_dofs_2; ++i)
1083  dof_handler.object_dof_indices[l][d].push_back(
1085  }
1086 
1087  // mark this face as visited
1088  face_touched[cell->face(face)->index()] = true;
1089  }
1090 
1091  for (unsigned int i = 1;
1092  i < dof_handler.object_dof_ptr[l][d].size();
1093  i++)
1094  dof_handler.object_dof_ptr[l][d][i] +=
1095  dof_handler.object_dof_ptr[l][d][i - 1];
1096  }
1097  }
1098 
1099 
1100 
1107  template <int spacedim>
1108  static void
1110  {
1111  Assert(dof_handler.fe_collection.size() > 0,
1113  Assert(dof_handler.tria->n_levels() > 0,
1114  ExcMessage("The current Triangulation must not be empty."));
1115  Assert(dof_handler.tria->n_levels() ==
1116  dof_handler.hp_cell_future_fe_indices.size(),
1117  ExcInternalError());
1118 
1120  reset_to_empty_objects(dof_handler);
1121 
1122  Threads::TaskGroup<> tasks;
1123  tasks +=
1124  Threads::new_task(&reserve_space_cells<1, spacedim>, dof_handler);
1125  tasks += Threads::new_task(&reserve_space_vertices<1, spacedim>,
1126  dof_handler);
1127  tasks.join_all();
1128  }
1129 
1130 
1131 
1132  template <int spacedim>
1133  static void
1135  {
1136  Assert(dof_handler.fe_collection.size() > 0,
1138  Assert(dof_handler.tria->n_levels() > 0,
1139  ExcMessage("The current Triangulation must not be empty."));
1140  Assert(dof_handler.tria->n_levels() ==
1141  dof_handler.hp_cell_future_fe_indices.size(),
1142  ExcInternalError());
1143 
1145  reset_to_empty_objects(dof_handler);
1146 
1147  Threads::TaskGroup<> tasks;
1148  tasks +=
1149  Threads::new_task(&reserve_space_cells<2, spacedim>, dof_handler);
1150  tasks +=
1151  Threads::new_task(&reserve_space_faces<2, spacedim>, dof_handler);
1152  tasks += Threads::new_task(&reserve_space_vertices<2, spacedim>,
1153  dof_handler);
1154  tasks.join_all();
1155  }
1156 
1157 
1158 
1159  template <int spacedim>
1160  static void
1162  {
1163  Assert(dof_handler.fe_collection.size() > 0,
1165  Assert(dof_handler.tria->n_levels() > 0,
1166  ExcMessage("The current Triangulation must not be empty."));
1167  Assert(dof_handler.tria->n_levels() ==
1168  dof_handler.hp_cell_future_fe_indices.size(),
1169  ExcInternalError());
1170 
1172  reset_to_empty_objects(dof_handler);
1173 
1174  Threads::TaskGroup<> tasks;
1175  tasks +=
1176  Threads::new_task(&reserve_space_cells<3, spacedim>, dof_handler);
1177  tasks +=
1178  Threads::new_task(&reserve_space_faces<3, spacedim>, dof_handler);
1179  tasks += Threads::new_task(&reserve_space_vertices<3, spacedim>,
1180  dof_handler);
1181 
1182  // While the tasks above are running, we can turn to line dofs
1183 
1184  // the situation here is pretty much like with vertices:
1185  // there can be an arbitrary number of finite elements
1186  // associated with each line.
1187  //
1188  // the algorithm we use is somewhat similar to what we do in
1189  // reserve_space_vertices()
1190  {
1191  // what we do first is to set up an array in which we
1192  // record whether a line is associated with any of the
1193  // given fe's, by setting a bit. in a later step, we
1194  // then actually allocate memory for the required dofs
1195  std::vector<std::vector<bool>> line_fe_association(
1196  dof_handler.fe_collection.size(),
1197  std::vector<bool>(dof_handler.tria->n_raw_lines(), false));
1198 
1199  for (const auto &cell : dof_handler.active_cell_iterators())
1200  if (!cell->is_artificial())
1201  for (const auto l : cell->line_indices())
1202  line_fe_association[cell->active_fe_index()]
1203  [cell->line_index(l)] = true;
1204 
1205  // first check which of the lines is used at all,
1206  // i.e. is associated with a finite element. we do this
1207  // since not all lines may actually be used, in which
1208  // case we do not have to allocate any memory at all
1209  std::vector<bool> line_is_used(dof_handler.tria->n_raw_lines(),
1210  false);
1211  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1212  ++line)
1213  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
1214  ++fe)
1215  if (line_fe_association[fe][line] == true)
1216  {
1217  line_is_used[line] = true;
1218  break;
1219  }
1220 
1221 
1222 
1223  const unsigned int d = 1;
1224  const unsigned int l = 0;
1225 
1226  dof_handler.hp_object_fe_ptr[d].clear();
1227  dof_handler.hp_object_fe_indices[d].clear();
1228  dof_handler.object_dof_ptr[l][d].clear();
1229  dof_handler.object_dof_indices[l][d].clear();
1230 
1231  dof_handler.hp_object_fe_ptr[d].reserve(
1232  dof_handler.tria->n_raw_lines() + 1);
1233 
1234  unsigned int line_slots_needed = 0;
1235  unsigned int fe_slots_needed = 0;
1236 
1237  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1238  ++line)
1239  {
1240  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1241 
1242  if (line_is_used[line] == true)
1243  {
1244  for (unsigned int fe = 0;
1245  fe < dof_handler.fe_collection.size();
1246  ++fe)
1247  if (line_fe_association[fe][line] == true)
1248  {
1249  fe_slots_needed++;
1250  line_slots_needed +=
1251  dof_handler.get_fe(fe).n_dofs_per_line();
1252  }
1253  }
1254  }
1255 
1256  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1257 
1258  // make sure that all entries have been set
1259  AssertDimension(dof_handler.hp_object_fe_ptr[d].size(),
1260  dof_handler.tria->n_raw_lines() + 1);
1261 
1262  dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
1263  dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
1264 
1265  dof_handler.object_dof_indices[l][d].reserve(line_slots_needed);
1266 
1267  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1268  ++line)
1269  if (line_is_used[line] == true)
1270  {
1271  for (unsigned int fe = 0;
1272  fe < dof_handler.fe_collection.size();
1273  ++fe)
1274  if (line_fe_association[fe][line] == true)
1275  {
1276  dof_handler.hp_object_fe_indices[d].push_back(fe);
1277  dof_handler.object_dof_ptr[l][d].push_back(
1278  dof_handler.object_dof_indices[l][d].size());
1279 
1280  for (unsigned int i = 0;
1281  i < dof_handler.get_fe(fe).n_dofs_per_line();
1282  i++)
1283  dof_handler.object_dof_indices[l][d].push_back(
1285  }
1286  }
1287 
1288  dof_handler.object_dof_ptr[l][d].push_back(
1289  dof_handler.object_dof_indices[l][d].size());
1290 
1291  // make sure that all entries have been set
1292  AssertDimension(dof_handler.hp_object_fe_indices[d].size(),
1293  fe_slots_needed);
1294  AssertDimension(dof_handler.object_dof_ptr[l][d].size(),
1295  fe_slots_needed + 1);
1296  AssertDimension(dof_handler.object_dof_indices[l][d].size(),
1297  line_slots_needed);
1298  }
1299 
1300  // Ensure that everything is done at this point.
1301  tasks.join_all();
1302  }
1303 
1304 
1305 
1317  template <int dim, int spacedim>
1318  static void
1320  {
1321  Assert(
1322  dof_handler.hp_capability_enabled == true,
1324 
1325  using active_fe_index_type =
1326  typename ::DoFHandler<dim, spacedim>::active_fe_index_type;
1327 
1328  if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1329  dynamic_cast<
1330  const ::parallel::shared::Triangulation<dim, spacedim>
1331  *>(&dof_handler.get_triangulation()))
1332  {
1333  // we have a shared triangulation. in this case, every processor
1334  // knows about all cells, but every processor only has knowledge
1335  // about the active FE index on the cells it owns.
1336  //
1337  // we can create a complete set of active FE indices by letting
1338  // every processor create a vector of indices for all cells,
1339  // filling only those on the cells it owns and setting the indices
1340  // on the other cells to zero. then we add all of these vectors
1341  // up, and because every vector entry has exactly one processor
1342  // that owns it, the sum is correct
1343  std::vector<active_fe_index_type> active_fe_indices(
1344  tr->n_active_cells(), 0u);
1345  for (const auto &cell : dof_handler.active_cell_iterators())
1346  if (cell->is_locally_owned())
1347  active_fe_indices[cell->active_cell_index()] =
1348  cell->active_fe_index();
1349 
1350  Utilities::MPI::sum(active_fe_indices,
1351  tr->get_communicator(),
1352  active_fe_indices);
1353 
1354  // now go back and fill the active FE index on all other
1355  // cells. we would like to call cell->set_active_fe_index(),
1356  // but that function does not allow setting these indices on
1357  // non-locally_owned cells. so we have to work around the
1358  // issue a little bit by accessing the underlying data
1359  // structures directly
1360  for (const auto &cell : dof_handler.active_cell_iterators())
1361  if (!cell->is_locally_owned())
1362  dof_handler
1363  .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1364  active_fe_indices[cell->active_cell_index()];
1365  }
1366  else if (const ::parallel::
1367  DistributedTriangulationBase<dim, spacedim> *tr =
1368  dynamic_cast<
1369  const ::parallel::
1370  DistributedTriangulationBase<dim, spacedim> *>(
1371  &dof_handler.get_triangulation()))
1372  {
1373  // For completely distributed meshes, use the function that is
1374  // able to move data from locally owned cells on one processor to
1375  // the corresponding ghost cells on others. To this end, we need
1376  // to have functions that can pack and unpack the data we want to
1377  // transport -- namely, the single unsigned int active_fe_index
1378  // objects
1379  auto pack =
1380  [](const typename ::DoFHandler<dim, spacedim>::
1381  active_cell_iterator &cell) -> active_fe_index_type {
1382  return cell->active_fe_index();
1383  };
1384 
1385  auto unpack =
1386  [&dof_handler](
1387  const typename ::DoFHandler<dim, spacedim>::
1388  active_cell_iterator & cell,
1389  const active_fe_index_type active_fe_index) -> void {
1390  // we would like to say
1391  // cell->set_active_fe_index(active_fe_index);
1392  // but this is not allowed on cells that are not
1393  // locally owned, and we are on a ghost cell
1394  dof_handler
1395  .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1396  active_fe_index;
1397  };
1398 
1400  active_fe_index_type,
1401  ::DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1402  }
1403  else
1404  {
1405  // a sequential triangulation. there is nothing we need to do here
1406  Assert(
1407  (dynamic_cast<
1408  const ::parallel::TriangulationBase<dim, spacedim> *>(
1409  &dof_handler.get_triangulation()) == nullptr),
1410  ExcInternalError());
1411  }
1412  }
1413 
1414 
1415 
1429  template <int dim, int spacedim>
1430  static void
1432  {
1433  Assert(
1434  dof_handler.hp_capability_enabled == true,
1436 
1437  using active_fe_index_type =
1438  typename ::DoFHandler<dim, spacedim>::active_fe_index_type;
1439 
1440  if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1441  dynamic_cast<
1442  const ::parallel::shared::Triangulation<dim, spacedim>
1443  *>(&dof_handler.get_triangulation()))
1444  {
1445  std::vector<active_fe_index_type> future_fe_indices(
1446  tr->n_active_cells(), 0u);
1447  for (const auto &cell : dof_handler.active_cell_iterators() |
1449  future_fe_indices[cell->active_cell_index()] =
1450  dof_handler
1451  .hp_cell_future_fe_indices[cell->level()][cell->index()];
1452 
1453  Utilities::MPI::sum(future_fe_indices,
1454  tr->get_communicator(),
1455  future_fe_indices);
1456 
1457  for (const auto &cell : dof_handler.active_cell_iterators())
1458  if (!cell->is_locally_owned())
1459  dof_handler
1460  .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1461  future_fe_indices[cell->active_cell_index()];
1462  }
1463  else if (const ::parallel::
1464  DistributedTriangulationBase<dim, spacedim> *tr =
1465  dynamic_cast<
1466  const ::parallel::
1467  DistributedTriangulationBase<dim, spacedim> *>(
1468  &dof_handler.get_triangulation()))
1469  {
1470  auto pack =
1471  [&dof_handler](
1472  const typename ::DoFHandler<dim, spacedim>::
1473  active_cell_iterator &cell) -> active_fe_index_type {
1474  return dof_handler
1475  .hp_cell_future_fe_indices[cell->level()][cell->index()];
1476  };
1477 
1478  auto unpack =
1479  [&dof_handler](
1480  const typename ::DoFHandler<dim, spacedim>::
1481  active_cell_iterator & cell,
1482  const active_fe_index_type future_fe_index) -> void {
1483  dof_handler
1484  .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1485  future_fe_index;
1486  };
1487 
1489  active_fe_index_type,
1490  ::DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1491  }
1492  else
1493  {
1494  Assert(
1495  (dynamic_cast<
1496  const ::parallel::TriangulationBase<dim, spacedim> *>(
1497  &dof_handler.get_triangulation()) == nullptr),
1498  ExcInternalError());
1499  }
1500  }
1501 
1502 
1503 
1524  template <int dim, int spacedim>
1525  static void
1527  DoFHandler<dim, spacedim> &dof_handler)
1528  {
1529  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1530 
1531  for (const auto &cell : dof_handler.active_cell_iterators())
1532  if (cell->is_locally_owned())
1533  {
1534  if (cell->refine_flag_set())
1535  {
1536  // Store the active FE index of each cell that will be
1537  // refined to and distribute it later on its children.
1538  // Pick their future index if flagged for p-refinement.
1539  fe_transfer->refined_cells_fe_index.insert(
1540  {cell, cell->future_fe_index()});
1541  }
1542  else if (cell->coarsen_flag_set())
1543  {
1544  // From all cells that will be coarsened, determine their
1545  // parent and calculate its proper active FE index, so that
1546  // it can be set after refinement. But first, check if that
1547  // particular cell has a parent at all.
1548  Assert(cell->level() > 0, ExcInternalError());
1549  const auto &parent = cell->parent();
1550 
1551  // Check if the active FE index for the current cell has
1552  // been determined already.
1553  if (fe_transfer->coarsened_cells_fe_index.find(parent) ==
1554  fe_transfer->coarsened_cells_fe_index.end())
1555  {
1556  // Find a suitable active FE index for the parent cell
1557  // based on the 'least dominant finite element' of its
1558  // children. Consider the childrens' hypothetical future
1559  // index when they have been flagged for p-refinement.
1560 #ifdef DEBUG
1561  for (const auto &child : parent->child_iterators())
1562  Assert(child->is_active() &&
1563  child->coarsen_flag_set(),
1564  typename ::Triangulation<
1565  dim>::ExcInconsistentCoarseningFlags());
1566 #endif
1567 
1568  const unsigned int fe_index = ::internal::hp::
1569  DoFHandlerImplementation::Implementation::
1570  dominated_future_fe_on_children<dim, spacedim>(
1571  parent);
1572 
1573  fe_transfer->coarsened_cells_fe_index.insert(
1574  {parent, fe_index});
1575  }
1576  }
1577  else
1578  {
1579  // No h-refinement is scheduled for this cell.
1580  // However, it may have p-refinement indicators, so we
1581  // choose a new active FE index based on its flags.
1582  if (cell->future_fe_index_set() == true)
1583  fe_transfer->persisting_cells_fe_index.insert(
1584  {cell, cell->future_fe_index()});
1585  }
1586  }
1587  }
1588 
1589 
1590 
1595  template <int dim, int spacedim>
1596  static void
1598  DoFHandler<dim, spacedim> &dof_handler)
1599  {
1600  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1601 
1602  // Set active FE indices on persisting cells.
1603  for (const auto &persist : fe_transfer->persisting_cells_fe_index)
1604  {
1605  const auto &cell = persist.first;
1606 
1607  if (cell->is_locally_owned())
1608  {
1609  Assert(cell->is_active(), ExcInternalError());
1610  cell->set_active_fe_index(persist.second);
1611  }
1612  }
1613 
1614  // Distribute active FE indices from all refined cells on their
1615  // respective children.
1616  for (const auto &refine : fe_transfer->refined_cells_fe_index)
1617  {
1618  const auto &parent = refine.first;
1619 
1620  for (const auto &child : parent->child_iterators())
1621  if (child->is_locally_owned())
1622  {
1623  Assert(child->is_active(), ExcInternalError());
1624  child->set_active_fe_index(refine.second);
1625  }
1626  }
1627 
1628  // Set active FE indices on coarsened cells that have been determined
1629  // before the actual coarsening happened.
1630  for (const auto &coarsen : fe_transfer->coarsened_cells_fe_index)
1631  {
1632  const auto &cell = coarsen.first;
1633 
1634  if (cell->is_locally_owned())
1635  {
1636  Assert(cell->is_active(), ExcInternalError());
1637  cell->set_active_fe_index(coarsen.second);
1638  }
1639  }
1640  }
1641 
1642 
1653  template <int dim, int spacedim>
1654  static unsigned int
1657  const std::vector<unsigned int> & children_fe_indices,
1658  const ::hp::FECollection<dim, spacedim> &fe_collection)
1659  {
1660  Assert(!children_fe_indices.empty(), ExcInternalError());
1661 
1662  // convert vector to set
1663  const std::set<unsigned int> children_fe_indices_set(
1664  children_fe_indices.begin(), children_fe_indices.end());
1665 
1666  const unsigned int dominated_fe_index =
1667  fe_collection.find_dominated_fe_extended(children_fe_indices_set,
1668  /*codim=*/0);
1669 
1670  Assert(dominated_fe_index != numbers::invalid_unsigned_int,
1672 
1673  return dominated_fe_index;
1674  }
1675 
1676 
1684  template <int dim, int spacedim>
1685  static unsigned int
1687  const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
1688  {
1689  Assert(
1690  !parent->is_active(),
1691  ExcMessage(
1692  "You ask for information on children of this cell which is only "
1693  "available for active cells. This cell has no children."));
1694 
1695  const auto &dof_handler = parent->get_dof_handler();
1696  Assert(
1697  dof_handler.has_hp_capabilities(),
1699 
1700  std::set<unsigned int> future_fe_indices_children;
1701  for (const auto &child : parent->child_iterators())
1702  {
1703  Assert(
1704  child->is_active(),
1705  ExcMessage(
1706  "You ask for information on children of this cell which is only "
1707  "available for active cells. One of its children is not active."));
1708 
1709  // Ghost siblings might occur on parallel::shared::Triangulation
1710  // objects. The public interface does not allow to access future
1711  // FE indices on ghost cells. However, we need this information
1712  // here and thus call the internal function that does not check
1713  // for cell ownership. This requires that future FE indices have
1714  // been communicated prior to calling this function.
1715  const unsigned int future_fe_index_child =
1716  ::internal::DoFCellAccessorImplementation::
1717  Implementation::future_fe_index<dim, spacedim, false>(*child);
1718 
1719  future_fe_indices_children.insert(future_fe_index_child);
1720  }
1721  Assert(!future_fe_indices_children.empty(), ExcInternalError());
1722 
1723  const unsigned int future_fe_index =
1724  dof_handler.fe_collection.find_dominated_fe_extended(
1725  future_fe_indices_children,
1726  /*codim=*/0);
1727 
1728  Assert(future_fe_index != numbers::invalid_unsigned_int,
1730 
1731  return future_fe_index;
1732  }
1733  };
1734 
1735 
1736 
1740  template <int dim, int spacedim>
1741  void
1743  {
1744  Implementation::communicate_future_fe_indices<dim, spacedim>(
1745  dof_handler);
1746  }
1747 
1748 
1749 
1753  template <int dim, int spacedim>
1754  unsigned int
1756  const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
1757  {
1758  return Implementation::dominated_future_fe_on_children<dim, spacedim>(
1759  parent);
1760  }
1761  } // namespace DoFHandlerImplementation
1762  } // namespace hp
1763 } // namespace internal
1764 
1765 
1766 
1767 template <int dim, int spacedim>
1769  : hp_capability_enabled(true)
1770  , tria(nullptr, typeid(*this).name())
1771  , mg_faces(nullptr)
1772 {}
1773 
1774 
1775 
1776 template <int dim, int spacedim>
1778  : DoFHandler()
1779 {
1780  reinit(tria);
1781 }
1782 
1783 
1784 
1785 template <int dim, int spacedim>
1787 {
1788  // unsubscribe all attachments to signals of the underlying triangulation
1789  for (auto &connection : this->tria_listeners)
1790  connection.disconnect();
1791  this->tria_listeners.clear();
1792 
1793  for (auto &connection : this->tria_listeners_for_transfer)
1794  connection.disconnect();
1795  this->tria_listeners_for_transfer.clear();
1796 
1797  // release allocated memory
1798  // virtual functions called in constructors and destructors never use the
1799  // override in a derived class
1800  // for clarity be explicit on which function is called
1802 
1803  // also release the policy. this needs to happen before the
1804  // current object disappears because the policy objects
1805  // store references to the DoFhandler object they work on
1806  this->policy.reset();
1807 }
1808 
1809 
1810 
1811 template <int dim, int spacedim>
1812 void
1814 {
1815  //
1816  // call destructor
1817  //
1818  // remove association with old triangulation
1819  for (auto &connection : this->tria_listeners)
1820  connection.disconnect();
1821  this->tria_listeners.clear();
1822 
1823  for (auto &connection : this->tria_listeners_for_transfer)
1824  connection.disconnect();
1825  this->tria_listeners_for_transfer.clear();
1826 
1827  // release allocated memory and policy
1829  this->policy.reset();
1830 
1831  // reset the finite element collection
1832  this->fe_collection = hp::FECollection<dim, spacedim>();
1833 
1834  //
1835  // call constructor
1836  //
1837  // establish connection to new triangulation
1838  this->tria = &tria;
1839  this->setup_policy();
1840 
1841  // start in hp-mode and let distribute_dofs toggle it if necessary
1842  hp_capability_enabled = true;
1843  this->connect_to_triangulation_signals();
1844  this->create_active_fe_table();
1845 }
1846 
1847 
1848 
1849 /*------------------------ Cell iterator functions ------------------------*/
1850 
1851 template <int dim, int spacedim>
1853 DoFHandler<dim, spacedim>::begin(const unsigned int level) const
1854 {
1856  this->get_triangulation().begin(level);
1857  if (cell == this->get_triangulation().end(level))
1858  return end(level);
1859  return cell_iterator(*cell, this);
1860 }
1861 
1862 
1863 
1864 template <int dim, int spacedim>
1866 DoFHandler<dim, spacedim>::begin_active(const unsigned int level) const
1867 {
1868  // level is checked in begin
1869  cell_iterator i = begin(level);
1870  if (i.state() != IteratorState::valid)
1871  return i;
1872  while (i->has_children())
1873  if ((++i).state() != IteratorState::valid)
1874  return i;
1875  return i;
1876 }
1877 
1878 
1879 
1880 template <int dim, int spacedim>
1883 {
1884  return cell_iterator(&this->get_triangulation(), -1, -1, this);
1885 }
1886 
1887 
1888 
1889 template <int dim, int spacedim>
1891 DoFHandler<dim, spacedim>::end(const unsigned int level) const
1892 {
1894  this->get_triangulation().end(level);
1895  if (cell.state() != IteratorState::valid)
1896  return end();
1897  return cell_iterator(*cell, this);
1898 }
1899 
1900 
1901 
1902 template <int dim, int spacedim>
1904 DoFHandler<dim, spacedim>::end_active(const unsigned int level) const
1905 {
1907  this->get_triangulation().end_active(level);
1908  if (cell.state() != IteratorState::valid)
1909  return active_cell_iterator(end());
1910  return active_cell_iterator(*cell, this);
1911 }
1912 
1913 
1914 
1915 template <int dim, int spacedim>
1917 DoFHandler<dim, spacedim>::begin_mg(const unsigned int level) const
1918 {
1919  Assert(this->has_level_dofs(),
1920  ExcMessage("You can only iterate over mg "
1921  "levels if mg dofs got distributed."));
1923  this->get_triangulation().begin(level);
1924  if (cell == this->get_triangulation().end(level))
1925  return end_mg(level);
1926  return level_cell_iterator(*cell, this);
1927 }
1928 
1929 
1930 
1931 template <int dim, int spacedim>
1933 DoFHandler<dim, spacedim>::end_mg(const unsigned int level) const
1934 {
1935  Assert(this->has_level_dofs(),
1936  ExcMessage("You can only iterate over mg "
1937  "levels if mg dofs got distributed."));
1939  this->get_triangulation().end(level);
1940  if (cell.state() != IteratorState::valid)
1941  return end();
1942  return level_cell_iterator(*cell, this);
1943 }
1944 
1945 
1946 
1947 template <int dim, int spacedim>
1950 {
1951  return level_cell_iterator(&this->get_triangulation(), -1, -1, this);
1952 }
1953 
1954 
1955 
1956 template <int dim, int spacedim>
1959 {
1961  begin(), end());
1962 }
1963 
1964 
1965 
1966 template <int dim, int spacedim>
1969 {
1970  return IteratorRange<
1971  typename DoFHandler<dim, spacedim>::active_cell_iterator>(begin_active(),
1972  end());
1973 }
1974 
1975 
1976 
1977 template <int dim, int spacedim>
1980 {
1982  begin_mg(), end_mg());
1983 }
1984 
1985 
1986 
1987 template <int dim, int spacedim>
1990  const unsigned int level) const
1991 {
1993  begin(level), end(level));
1994 }
1995 
1996 
1997 
1998 template <int dim, int spacedim>
2001  const unsigned int level) const
2002 {
2003  return IteratorRange<
2005  begin_active(level), end_active(level));
2006 }
2007 
2008 
2009 
2010 template <int dim, int spacedim>
2013  const unsigned int level) const
2014 {
2016  begin_mg(level), end_mg(level));
2017 }
2018 
2019 
2020 
2021 //---------------------------------------------------------------------------
2022 
2023 
2024 
2025 template <int dim, int spacedim>
2028 {
2029  Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2030  ExcNotImplementedWithHP());
2031 
2032  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2033 
2034  std::unordered_set<types::global_dof_index> boundary_dofs;
2035  std::vector<types::global_dof_index> dofs_on_face;
2036  dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2037 
2038  const IndexSet &owned_dofs = locally_owned_dofs();
2039 
2040  // loop over all faces to check whether they are at a
2041  // boundary. note that we need not take special care of single
2042  // lines in 3d (using @p{cell->has_boundary_lines}), since we do
2043  // not support boundaries of dimension dim-2, and so every
2044  // boundary line is also part of a boundary face.
2045  for (const auto &cell : this->active_cell_iterators())
2046  if (cell->is_locally_owned() && cell->at_boundary())
2047  {
2048  for (const auto iface : cell->face_indices())
2049  {
2050  const auto face = cell->face(iface);
2051  if (face->at_boundary())
2052  {
2053  const unsigned int dofs_per_face =
2054  cell->get_fe().n_dofs_per_face(iface);
2055  dofs_on_face.resize(dofs_per_face);
2056 
2057  face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2058  for (unsigned int i = 0; i < dofs_per_face; ++i)
2059  {
2060  const unsigned int global_idof_index = dofs_on_face[i];
2061  if (owned_dofs.is_element(global_idof_index))
2062  {
2063  boundary_dofs.insert(global_idof_index);
2064  }
2065  }
2066  }
2067  }
2068  }
2069  return boundary_dofs.size();
2070 }
2071 
2072 
2073 
2074 template <int dim, int spacedim>
2077  const std::set<types::boundary_id> &boundary_ids) const
2078 {
2079  Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2080  ExcNotImplementedWithHP());
2081 
2082  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2083  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2084  boundary_ids.end(),
2086 
2087  // same as above, but with additional checks for set of boundary
2088  // indicators
2089  std::unordered_set<types::global_dof_index> boundary_dofs;
2090  std::vector<types::global_dof_index> dofs_on_face;
2091  dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2092 
2093  const IndexSet &owned_dofs = locally_owned_dofs();
2094 
2095  for (const auto &cell : this->active_cell_iterators())
2096  if (cell->is_locally_owned() && cell->at_boundary())
2097  {
2098  for (const auto iface : cell->face_indices())
2099  {
2100  const auto face = cell->face(iface);
2101  const unsigned int boundary_id = face->boundary_id();
2102  if (face->at_boundary() &&
2103  (boundary_ids.find(boundary_id) != boundary_ids.end()))
2104  {
2105  const unsigned int dofs_per_face =
2106  cell->get_fe().n_dofs_per_face(iface);
2107  dofs_on_face.resize(dofs_per_face);
2108 
2109  face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2110  for (unsigned int i = 0; i < dofs_per_face; ++i)
2111  {
2112  const unsigned int global_idof_index = dofs_on_face[i];
2113  if (owned_dofs.is_element(global_idof_index))
2114  {
2115  boundary_dofs.insert(global_idof_index);
2116  }
2117  }
2118  }
2119  }
2120  }
2121  return boundary_dofs.size();
2122 }
2123 
2124 
2125 
2126 template <int dim, int spacedim>
2127 std::size_t
2129 {
2130  std::size_t mem = MemoryConsumption::memory_consumption(this->tria) +
2131  MemoryConsumption::memory_consumption(this->fe_collection) +
2132  MemoryConsumption::memory_consumption(this->number_cache);
2133 
2134  mem += MemoryConsumption::memory_consumption(cell_dof_cache_indices) +
2135  MemoryConsumption::memory_consumption(cell_dof_cache_ptr) +
2136  MemoryConsumption::memory_consumption(object_dof_indices) +
2137  MemoryConsumption::memory_consumption(object_dof_ptr) +
2138  MemoryConsumption::memory_consumption(hp_object_fe_indices) +
2139  MemoryConsumption::memory_consumption(hp_object_fe_ptr) +
2140  MemoryConsumption::memory_consumption(hp_cell_active_fe_indices) +
2141  MemoryConsumption::memory_consumption(hp_cell_future_fe_indices);
2142 
2143 
2144  if (hp_capability_enabled)
2145  {
2146  // nothing to add
2147  }
2148  else
2149  {
2150  // collect size of multigrid data structures
2151 
2152  mem += MemoryConsumption::memory_consumption(this->block_info_object);
2153 
2154  for (unsigned int level = 0; level < this->mg_levels.size(); ++level)
2155  mem += this->mg_levels[level]->memory_consumption();
2156 
2157  if (this->mg_faces != nullptr)
2158  mem += MemoryConsumption::memory_consumption(*this->mg_faces);
2159 
2160  for (unsigned int i = 0; i < this->mg_vertex_dofs.size(); ++i)
2161  mem += sizeof(MGVertexDoFs) +
2162  (1 + this->mg_vertex_dofs[i].get_finest_level() -
2163  this->mg_vertex_dofs[i].get_coarsest_level()) *
2164  sizeof(types::global_dof_index);
2165  }
2166 
2167  return mem;
2168 }
2169 
2170 
2171 
2172 template <int dim, int spacedim>
2173 void
2175  const FiniteElement<dim, spacedim> &fe)
2176 {
2177  this->distribute_dofs(hp::FECollection<dim, spacedim>(fe));
2178 }
2179 
2180 
2181 
2182 template <int dim, int spacedim>
2183 void
2186 {
2187  Assert(
2188  this->tria != nullptr,
2189  ExcMessage(
2190  "You need to set the Triangulation in the DoFHandler using reinit() or "
2191  "in the constructor before you can distribute DoFs."));
2192  Assert(this->tria->n_levels() > 0,
2193  ExcMessage("The Triangulation you are using is empty!"));
2194  Assert(ff.size() > 0, ExcMessage("The hp::FECollection given is empty!"));
2195 
2196  //
2197  // register the new finite element collection
2198  //
2199  // don't create a new object if the one we have is identical
2200  if (this->fe_collection != ff)
2201  {
2202  this->fe_collection = hp::FECollection<dim, spacedim>(ff);
2203 
2204  const bool contains_multiple_fes = (this->fe_collection.size() > 1);
2205 
2206  // disable hp-mode if only a single finite element has been registered
2207  if (hp_capability_enabled && !contains_multiple_fes)
2208  {
2209  hp_capability_enabled = false;
2210 
2211  // unsubscribe connections to signals that are only relevant for
2212  // hp-mode, since we only have a single element here
2213  for (auto &connection : this->tria_listeners_for_transfer)
2214  connection.disconnect();
2215  this->tria_listeners_for_transfer.clear();
2216 
2217  // release active and future finite element tables
2218  this->hp_cell_active_fe_indices.clear();
2219  this->hp_cell_active_fe_indices.shrink_to_fit();
2220  this->hp_cell_future_fe_indices.clear();
2221  this->hp_cell_future_fe_indices.shrink_to_fit();
2222  }
2223 
2224  // re-enabling hp-mode is not permitted since the active and future FE
2225  // tables are no longer available
2226  AssertThrow(
2227  hp_capability_enabled || !contains_multiple_fes,
2228  ExcMessage(
2229  "You cannot re-enable hp-capabilities after you registered a single "
2230  "finite element. Please call reinit() or create a new DoFHandler "
2231  "object instead."));
2232  }
2233 
2234  //
2235  // enumerate all degrees of freedom
2236  //
2237  if (hp_capability_enabled)
2238  {
2239  // make sure every processor knows the active FE indices
2240  // on both its own cells and all ghost cells
2243 
2244 #ifdef DEBUG
2245  // make sure that the FE collection is large enough to
2246  // cover all FE indices presently in use on the mesh
2247  for (const auto &cell : this->active_cell_iterators())
2248  {
2249  if (!cell->is_artificial())
2250  Assert(cell->active_fe_index() < this->fe_collection.size(),
2251  ExcInvalidFEIndex(cell->active_fe_index(),
2252  this->fe_collection.size()));
2253  if (cell->is_locally_owned())
2254  Assert(cell->future_fe_index() < this->fe_collection.size(),
2255  ExcInvalidFEIndex(cell->future_fe_index(),
2256  this->fe_collection.size()));
2257  }
2258 #endif
2259  }
2260 
2261  {
2262  // We would like to enumerate all dofs for shared::Triangulations. If an
2263  // underlying shared::Tria allows artificial cells, we need to restore the
2264  // true cell owners temporarily.
2265  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
2266  // the current set of subdomain ids, set subdomain ids to the "true" owner
2267  // of each cell upon construction of the TemporarilyRestoreSubdomainIds
2268  // object, and later restore these flags when it is destroyed.
2270  spacedim>
2271  subdomain_modifier(this->get_triangulation());
2272 
2273  // Adjust size of levels to the triangulation. Note that we still have to
2274  // allocate space for all degrees of freedom on this mesh (including ghost
2275  // and cells that are entirely stored on different processors), though we
2276  // may not assign numbers to some of them (i.e. they will remain at
2277  // invalid_dof_index). We need to allocate the space because we will want
2278  // to be able to query the dof_indices on each cell, and simply be told
2279  // that we don't know them on some cell (i.e. get back invalid_dof_index)
2280  if (hp_capability_enabled)
2282  *this);
2283  else
2285  }
2286 
2287  // hand the actual work over to the policy
2288  this->number_cache = this->policy->distribute_dofs();
2289 
2290  // do some housekeeping: compress indices
2291  // if(hp_capability_enabled)
2292  // {
2293  // Threads::TaskGroup<> tg;
2294  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2295  // tg += Threads::new_task(
2296  // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2297  // *this->levels_hp[level],
2298  // this->fe_collection);
2299  // tg.join_all();
2300  // }
2301 
2302  // Initialize the block info object only if this is a sequential
2303  // triangulation. It doesn't work correctly yet if it is parallel and has not
2304  // yet been implemented for hp-mode.
2305  if (!hp_capability_enabled &&
2307  *>(&*this->tria) == nullptr)
2308  this->block_info_object.initialize(*this, false, true);
2309 }
2310 
2311 
2312 
2313 template <int dim, int spacedim>
2314 void
2316 {
2317  AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2318 
2319  Assert(
2320  this->object_dof_indices.size() > 0,
2321  ExcMessage(
2322  "Distribute active DoFs using distribute_dofs() before calling distribute_mg_dofs()."));
2323 
2324  Assert(
2325  ((this->tria->get_mesh_smoothing() &
2328  ExcMessage(
2329  "The mesh smoothing requirement 'limit_level_difference_at_vertices' has to be set for using multigrid!"));
2330 
2331  this->clear_mg_space();
2332 
2334  this->mg_number_cache = this->policy->distribute_mg_dofs();
2335 
2336  // initialize the block info object only if this is a sequential
2337  // triangulation. it doesn't work correctly yet if it is parallel
2338  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2339  &*this->tria) == nullptr)
2340  this->block_info_object.initialize(*this, true, false);
2341 }
2342 
2343 
2344 
2345 template <int dim, int spacedim>
2346 void
2348 {
2349  AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2350 
2351  this->block_info_object.initialize_local(*this);
2352 }
2353 
2354 
2355 
2356 template <int dim, int spacedim>
2357 void
2359 {
2360  // decide whether we need a sequential or a parallel distributed policy
2361  if (dynamic_cast<const ::parallel::shared::Triangulation<dim, spacedim>
2362  *>(&this->get_triangulation()) != nullptr)
2363  this->policy = std::make_unique<internal::DoFHandlerImplementation::Policy::
2364  ParallelShared<dim, spacedim>>(*this);
2365  else if (dynamic_cast<
2366  const ::parallel::DistributedTriangulationBase<dim, spacedim>
2367  *>(&this->get_triangulation()) == nullptr)
2368  this->policy = std::make_unique<
2370  *this);
2371  else
2372  this->policy =
2373  std::make_unique<internal::DoFHandlerImplementation::Policy::
2374  ParallelDistributed<dim, spacedim>>(*this);
2375 }
2376 
2377 
2378 
2379 template <int dim, int spacedim>
2380 void
2382 {
2383  // release memory
2384  this->clear_space();
2385  this->clear_mg_space();
2386 }
2387 
2388 
2389 
2390 template <int dim, int spacedim>
2391 void
2393 {
2394  cell_dof_cache_indices.clear();
2395 
2396  cell_dof_cache_ptr.clear();
2397 
2398  object_dof_indices.clear();
2399 
2400  object_dof_ptr.clear();
2401 
2402  this->number_cache.clear();
2403 
2404  this->hp_cell_active_fe_indices.clear();
2405  this->hp_cell_future_fe_indices.clear();
2406 }
2407 
2408 
2409 
2410 template <int dim, int spacedim>
2411 void
2413 {
2414  this->mg_levels.clear();
2415  this->mg_faces.reset();
2416 
2417  std::vector<MGVertexDoFs> tmp;
2418 
2419  std::swap(this->mg_vertex_dofs, tmp);
2420 
2421  this->mg_number_cache.clear();
2422 }
2423 
2424 
2425 
2426 template <int dim, int spacedim>
2427 void
2429  const std::vector<types::global_dof_index> &new_numbers)
2430 {
2431  if (hp_capability_enabled)
2432  {
2433  Assert(this->hp_cell_future_fe_indices.size() > 0,
2434  ExcMessage(
2435  "You need to distribute DoFs before you can renumber them."));
2436 
2437  AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2438 
2439 #ifdef DEBUG
2440  // assert that the new indices are consecutively numbered if we are
2441  // working on a single processor. this doesn't need to
2442  // hold in the case of a parallel mesh since we map the interval
2443  // [0...n_dofs()) into itself but only globally, not on each processor
2444  if (this->n_locally_owned_dofs() == this->n_dofs())
2445  {
2446  std::vector<types::global_dof_index> tmp(new_numbers);
2447  std::sort(tmp.begin(), tmp.end());
2448  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2450  for (; p != tmp.end(); ++p, ++i)
2451  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2452  }
2453  else
2454  for (const auto new_number : new_numbers)
2455  Assert(new_number < this->n_dofs(),
2456  ExcMessage(
2457  "New DoF index is not less than the total number of dofs."));
2458 #endif
2459 
2460  // uncompress the internal storage scheme of dofs on cells so that
2461  // we can access dofs in turns. uncompress in parallel, starting
2462  // with the most expensive levels (the highest ones)
2463  //{
2464  // Threads::TaskGroup<> tg;
2465  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2466  // tg += Threads::new_task(
2467  // &::internal::hp::DoFLevel::uncompress_data<dim, spacedim>,
2468  // *this->levels_hp[level],
2469  // this->fe_collection);
2470  // tg.join_all();
2471  //}
2472 
2473  // do the renumbering
2474  this->number_cache = this->policy->renumber_dofs(new_numbers);
2475 
2476  // now re-compress the dof indices
2477  //{
2478  // Threads::TaskGroup<> tg;
2479  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2480  // tg += Threads::new_task(
2481  // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2482  // *this->levels_hp[level],
2483  // this->fe_collection);
2484  // tg.join_all();
2485  //}
2486  }
2487  else
2488  {
2489  Assert(this->object_dof_indices.size() > 0,
2490  ExcMessage(
2491  "You need to distribute DoFs before you can renumber them."));
2492 
2493 #ifdef DEBUG
2494  if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
2495  &*this->tria) != nullptr)
2496  {
2497  Assert(new_numbers.size() == this->n_dofs() ||
2498  new_numbers.size() == this->n_locally_owned_dofs(),
2499  ExcMessage("Incorrect size of the input array."));
2500  }
2501  else if (dynamic_cast<
2503  &*this->tria) != nullptr)
2504  {
2505  AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2506  }
2507  else
2508  {
2509  AssertDimension(new_numbers.size(), this->n_dofs());
2510  }
2511 
2512  // assert that the new indices are consecutively numbered if we are
2513  // working on a single processor. this doesn't need to
2514  // hold in the case of a parallel mesh since we map the interval
2515  // [0...n_dofs()) into itself but only globally, not on each processor
2516  if (this->n_locally_owned_dofs() == this->n_dofs())
2517  {
2518  std::vector<types::global_dof_index> tmp(new_numbers);
2519  std::sort(tmp.begin(), tmp.end());
2520  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2522  for (; p != tmp.end(); ++p, ++i)
2523  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2524  }
2525  else
2526  for (const auto new_number : new_numbers)
2527  Assert(new_number < this->n_dofs(),
2528  ExcMessage(
2529  "New DoF index is not less than the total number of dofs."));
2530 #endif
2531 
2532  this->number_cache = this->policy->renumber_dofs(new_numbers);
2533  }
2534 }
2535 
2536 
2537 
2538 template <int dim, int spacedim>
2539 void
2541  const unsigned int level,
2542  const std::vector<types::global_dof_index> &new_numbers)
2543 {
2544  AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2545 
2546  Assert(
2547  this->mg_levels.size() > 0 && this->object_dof_indices.size() > 0,
2548  ExcMessage(
2549  "You need to distribute active and level DoFs before you can renumber level DoFs."));
2550  AssertIndexRange(level, this->get_triangulation().n_global_levels());
2551  AssertDimension(new_numbers.size(),
2552  this->locally_owned_mg_dofs(level).n_elements());
2553 
2554 #ifdef DEBUG
2555  // assert that the new indices are consecutively numbered if we are working
2556  // on a single processor. this doesn't need to hold in the case of a
2557  // parallel mesh since we map the interval [0...n_dofs(level)) into itself
2558  // but only globally, not on each processor
2559  if (this->n_locally_owned_dofs() == this->n_dofs())
2560  {
2561  std::vector<types::global_dof_index> tmp(new_numbers);
2562  std::sort(tmp.begin(), tmp.end());
2563  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2565  for (; p != tmp.end(); ++p, ++i)
2566  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2567  }
2568  else
2569  for (const auto new_number : new_numbers)
2570  Assert(new_number < this->n_dofs(level),
2571  ExcMessage(
2572  "New DoF index is not less than the total number of dofs."));
2573 #endif
2574 
2575  this->mg_number_cache[level] =
2576  this->policy->renumber_mg_dofs(level, new_numbers);
2577 }
2578 
2579 
2580 
2581 template <int dim, int spacedim>
2582 unsigned int
2584 {
2585  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2586 
2587  switch (dim)
2588  {
2589  case 1:
2590  return this->fe_collection.max_dofs_per_vertex();
2591  case 2:
2592  return (3 * this->fe_collection.max_dofs_per_vertex() +
2593  2 * this->fe_collection.max_dofs_per_line());
2594  case 3:
2595  // we need to take refinement of one boundary face into
2596  // consideration here; in fact, this function returns what
2597  // #max_coupling_between_dofs<2> returns
2598  //
2599  // we assume here, that only four faces meet at the boundary;
2600  // this assumption is not justified and needs to be fixed some
2601  // time. fortunately, omitting it for now does no harm since
2602  // the matrix will cry foul if its requirements are not
2603  // satisfied
2604  return (19 * this->fe_collection.max_dofs_per_vertex() +
2605  28 * this->fe_collection.max_dofs_per_line() +
2606  8 * this->fe_collection.max_dofs_per_quad());
2607  default:
2608  Assert(false, ExcNotImplemented());
2609  return 0;
2610  }
2611 }
2612 
2613 
2614 
2615 template <int dim, int spacedim>
2616 unsigned int
2618 {
2619  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2622 }
2623 
2624 
2625 
2626 template <int dim, int spacedim>
2627 void
2629  const std::vector<unsigned int> &active_fe_indices)
2630 {
2631  Assert(active_fe_indices.size() == this->get_triangulation().n_active_cells(),
2632  ExcDimensionMismatch(active_fe_indices.size(),
2633  this->get_triangulation().n_active_cells()));
2634 
2635  this->create_active_fe_table();
2636  // we could set the values directly, since they are stored as
2637  // protected data of this object, but for simplicity we use the
2638  // cell-wise access. this way we also have to pass some debug-mode
2639  // tests which we would have to duplicate ourselves otherwise
2640  for (const auto &cell : this->active_cell_iterators())
2641  if (cell->is_locally_owned())
2642  cell->set_active_fe_index(active_fe_indices[cell->active_cell_index()]);
2643 }
2644 
2645 
2646 
2647 template <int dim, int spacedim>
2648 void
2650  std::vector<unsigned int> &active_fe_indices) const
2651 {
2652  active_fe_indices.resize(this->get_triangulation().n_active_cells());
2653 
2654  // we could try to extract the values directly, since they are
2655  // stored as protected data of this object, but for simplicity we
2656  // use the cell-wise access.
2657  for (const auto &cell : this->active_cell_iterators())
2658  if (!cell->is_artificial())
2659  active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
2660 }
2661 
2662 
2663 
2664 template <int dim, int spacedim>
2665 void
2667 {
2668  // make sure this is called during initialization in hp-mode
2669  Assert(hp_capability_enabled, ExcOnlyAvailableWithHP());
2670 
2671  // connect functions to signals of the underlying triangulation
2672  this->tria_listeners.push_back(this->tria->signals.create.connect(
2673  [this]() { this->reinit(*(this->tria)); }));
2674  this->tria_listeners.push_back(
2675  this->tria->signals.clear.connect([this]() { this->clear(); }));
2676 
2677  // attach corresponding callback functions dealing with the transfer of
2678  // active FE indices depending on the type of triangulation
2679  if (dynamic_cast<
2680  const ::parallel::fullydistributed::Triangulation<dim, spacedim>
2681  *>(&this->get_triangulation()))
2682  {
2683  // no transfer of active FE indices for this class
2684  }
2685  else if (dynamic_cast<
2686  const ::parallel::distributed::Triangulation<dim, spacedim>
2687  *>(&this->get_triangulation()))
2688  {
2689  // repartitioning signals
2690  this->tria_listeners_for_transfer.push_back(
2691  this->tria->signals.pre_distributed_repartition.connect([this]() {
2692  internal::hp::DoFHandlerImplementation::Implementation::
2693  ensure_absence_of_future_fe_indices<dim, spacedim>(*this);
2694  }));
2695  this->tria_listeners_for_transfer.push_back(
2696  this->tria->signals.pre_distributed_repartition.connect(
2697  [this]() { this->pre_distributed_transfer_action(); }));
2698  this->tria_listeners_for_transfer.push_back(
2699  this->tria->signals.post_distributed_repartition.connect(
2700  [this] { this->post_distributed_transfer_action(); }));
2701 
2702  // refinement signals
2703  this->tria_listeners_for_transfer.push_back(
2704  this->tria->signals.post_p4est_refinement.connect(
2705  [this]() { this->pre_distributed_transfer_action(); }));
2706  this->tria_listeners_for_transfer.push_back(
2707  this->tria->signals.post_distributed_refinement.connect(
2708  [this]() { this->post_distributed_transfer_action(); }));
2709 
2710  // serialization signals
2711  this->tria_listeners_for_transfer.push_back(
2712  this->tria->signals.post_distributed_save.connect(
2713  [this]() { this->active_fe_index_transfer.reset(); }));
2714  this->tria_listeners_for_transfer.push_back(
2715  this->tria->signals.post_distributed_load.connect(
2716  [this]() { this->update_active_fe_table(); }));
2717  }
2718  else if (dynamic_cast<
2719  const ::parallel::shared::Triangulation<dim, spacedim> *>(
2720  &this->get_triangulation()) != nullptr)
2721  {
2722  // partitioning signals
2723  this->tria_listeners_for_transfer.push_back(
2724  this->tria->signals.pre_partition.connect([this]() {
2725  internal::hp::DoFHandlerImplementation::Implementation::
2726  ensure_absence_of_future_fe_indices(*this);
2727  }));
2728 
2729  // refinement signals
2730  this->tria_listeners_for_transfer.push_back(
2731  this->tria->signals.pre_refinement.connect([this]() {
2732  internal::hp::DoFHandlerImplementation::Implementation::
2733  communicate_future_fe_indices(*this);
2734  }));
2735  this->tria_listeners_for_transfer.push_back(
2736  this->tria->signals.pre_refinement.connect(
2737  [this] { this->pre_transfer_action(); }));
2738  this->tria_listeners_for_transfer.push_back(
2739  this->tria->signals.post_refinement.connect(
2740  [this] { this->post_transfer_action(); }));
2741  }
2742  else
2743  {
2744  // refinement signals
2745  this->tria_listeners_for_transfer.push_back(
2746  this->tria->signals.pre_refinement.connect(
2747  [this] { this->pre_transfer_action(); }));
2748  this->tria_listeners_for_transfer.push_back(
2749  this->tria->signals.post_refinement.connect(
2750  [this] { this->post_transfer_action(); }));
2751  }
2752 }
2753 
2754 
2755 
2756 template <int dim, int spacedim>
2757 void
2759 {
2760  AssertThrow(hp_capability_enabled == true, ExcOnlyAvailableWithHP());
2761 
2762 
2763  // Create sufficiently many hp::DoFLevels.
2764  // while (this->levels_hp.size() < this->tria->n_levels())
2765  // this->levels_hp.emplace_back(new ::internal::hp::DoFLevel);
2766 
2767  this->hp_cell_active_fe_indices.resize(this->tria->n_levels());
2768  this->hp_cell_future_fe_indices.resize(this->tria->n_levels());
2769 
2770  // then make sure that on each level we have the appropriate size
2771  // of active FE indices; preset them to zero, i.e. the default FE
2772  for (unsigned int level = 0; level < this->hp_cell_future_fe_indices.size();
2773  ++level)
2774  {
2775  if (this->hp_cell_active_fe_indices[level].size() == 0 &&
2776  this->hp_cell_future_fe_indices[level].size() == 0)
2777  {
2778  this->hp_cell_active_fe_indices[level].resize(
2779  this->tria->n_raw_cells(level), 0);
2780  this->hp_cell_future_fe_indices[level].resize(
2781  this->tria->n_raw_cells(level), invalid_active_fe_index);
2782  }
2783  else
2784  {
2785  // Either the active FE indices have size zero because
2786  // they were just created, or the correct size. Other
2787  // sizes indicate that something went wrong.
2788  Assert(this->hp_cell_active_fe_indices[level].size() ==
2789  this->tria->n_raw_cells(level) &&
2790  this->hp_cell_future_fe_indices[level].size() ==
2791  this->tria->n_raw_cells(level),
2792  ExcInternalError());
2793  }
2794 
2795  // it may be that the previous table was compressed; in that
2796  // case, restore the correct active FE index. the fact that
2797  // this no longer matches the indices in the table is of no
2798  // importance because the current function is called at a
2799  // point where we have to recreate the dof_indices tables in
2800  // the levels anyway
2801  // this->levels_hp[level]->normalize_active_fe_indices();
2802  }
2803 }
2804 
2805 
2806 
2807 template <int dim, int spacedim>
2808 void
2810 {
2811  // // Normally only one level is added, but if this Triangulation
2812  // // is created by copy_triangulation, it can be more than one level.
2813  // while (this->levels_hp.size() < this->tria->n_levels())
2814  // this->levels_hp.emplace_back(new ::internal::hp::DoFLevel);
2815  //
2816  // // Coarsening can lead to the loss of levels. Hence remove them.
2817  // while (this->levels_hp.size() > this->tria->n_levels())
2818  // {
2819  // // drop the last element. that also releases the memory pointed to
2820  // this->levels_hp.pop_back();
2821  // }
2822 
2823  this->hp_cell_active_fe_indices.resize(this->tria->n_levels());
2824  this->hp_cell_active_fe_indices.shrink_to_fit();
2825 
2826  this->hp_cell_future_fe_indices.resize(this->tria->n_levels());
2827  this->hp_cell_future_fe_indices.shrink_to_fit();
2828 
2829  for (unsigned int i = 0; i < this->hp_cell_future_fe_indices.size(); ++i)
2830  {
2831  // Resize active FE indices vectors. Use zero indicator to extend.
2832  this->hp_cell_active_fe_indices[i].resize(this->tria->n_raw_cells(i), 0);
2833 
2834  // Resize future FE indices vectors. Make sure that all
2835  // future FE indices have been cleared after refinement happened.
2836  //
2837  // We have used future FE indices to update all active FE indices
2838  // before refinement happened, thus we are safe to clear them now.
2839  this->hp_cell_future_fe_indices[i].assign(this->tria->n_raw_cells(i),
2840  invalid_active_fe_index);
2841  }
2842 }
2843 
2844 
2845 template <int dim, int spacedim>
2846 void
2848 {
2849  Assert(this->active_fe_index_transfer == nullptr, ExcInternalError());
2850 
2851  this->active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
2852 
2855 }
2856 
2857 
2858 
2859 template <int dim, int spacedim>
2860 void
2862 {
2863 #ifndef DEAL_II_WITH_P4EST
2864  Assert(false,
2865  ExcMessage(
2866  "You are attempting to use a functionality that is only available "
2867  "if deal.II was configured to use p4est, but cmake did not find a "
2868  "valid p4est library."));
2869 #else
2870  // the implementation below requires a p:d:T currently
2871  Assert(
2873  &this->get_triangulation()) != nullptr),
2874  ExcNotImplemented());
2875 
2876  Assert(active_fe_index_transfer == nullptr, ExcInternalError());
2877 
2878  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
2879 
2880  // If we work on a p::d::Triangulation, we have to transfer all
2881  // active FE indices since ownership of cells may change. We will
2882  // use our p::d::CellDataTransfer member to achieve this. Further,
2883  // we prepare the values in such a way that they will correspond to
2884  // the active FE indices on the new mesh.
2885 
2886  // Gather all current future FE indices.
2887  active_fe_index_transfer->active_fe_indices.resize(
2888  get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
2889 
2890  for (const auto &cell : active_cell_iterators())
2891  if (cell->is_locally_owned())
2892  active_fe_index_transfer->active_fe_indices[cell->active_cell_index()] =
2893  cell->future_fe_index();
2894 
2895  // Create transfer object and attach to it.
2896  const auto *distributed_tria =
2898  &this->get_triangulation());
2899 
2900  active_fe_index_transfer->cell_data_transfer = std::make_unique<
2901  parallel::distributed::
2902  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
2903  *distributed_tria,
2904  /*transfer_variable_size_data=*/false,
2905  /*refinement_strategy=*/
2906  &::AdaptationStrategies::Refinement::
2907  preserve<dim, spacedim, unsigned int>,
2908  /*coarsening_strategy=*/
2909  [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
2910  const std::vector<unsigned int> &children_fe_indices)
2911  -> unsigned int {
2912  return ::internal::hp::DoFHandlerImplementation::Implementation::
2913  determine_fe_from_children<dim, spacedim>(parent,
2914  children_fe_indices,
2915  fe_collection);
2916  });
2917 
2918  active_fe_index_transfer->cell_data_transfer
2919  ->prepare_for_coarsening_and_refinement(
2920  active_fe_index_transfer->active_fe_indices);
2921 #endif
2922 }
2923 
2924 
2925 
2926 template <int dim, int spacedim>
2927 void
2929 {
2930  update_active_fe_table();
2931 
2932  Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
2933 
2936 
2937  // We have to distribute the information about active FE indices
2938  // of all cells (including the artificial ones) on all processors,
2939  // if a parallel::shared::Triangulation has been used.
2942 
2943  // Free memory.
2944  this->active_fe_index_transfer.reset();
2945 }
2946 
2947 
2948 
2949 template <int dim, int spacedim>
2950 void
2952 {
2953 #ifndef DEAL_II_WITH_P4EST
2954  Assert(false, ExcInternalError());
2955 #else
2956  update_active_fe_table();
2957 
2958  Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
2959 
2960  // Unpack active FE indices.
2961  this->active_fe_index_transfer->active_fe_indices.resize(
2962  this->get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
2963  this->active_fe_index_transfer->cell_data_transfer->unpack(
2964  this->active_fe_index_transfer->active_fe_indices);
2965 
2966  // Update all locally owned active FE indices.
2967  this->set_active_fe_indices(
2968  this->active_fe_index_transfer->active_fe_indices);
2969 
2970  // Update active FE indices on ghost cells.
2973 
2974  // Free memory.
2975  this->active_fe_index_transfer.reset();
2976 #endif
2977 }
2978 
2979 
2980 
2981 template <int dim, int spacedim>
2982 void
2984 {
2985 #ifndef DEAL_II_WITH_P4EST
2986  Assert(false,
2987  ExcMessage(
2988  "You are attempting to use a functionality that is only available "
2989  "if deal.II was configured to use p4est, but cmake did not find a "
2990  "valid p4est library."));
2991 #else
2992  // the implementation below requires a p:d:T currently
2993  Assert(
2995  &this->get_triangulation()) != nullptr),
2996  ExcNotImplemented());
2997 
2998  Assert(active_fe_index_transfer == nullptr, ExcInternalError());
2999 
3000  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3001 
3002  // Create transfer object and attach to it.
3003  const auto *distributed_tria =
3005  &this->get_triangulation());
3006 
3007  active_fe_index_transfer->cell_data_transfer = std::make_unique<
3008  parallel::distributed::
3009  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3010  *distributed_tria,
3011  /*transfer_variable_size_data=*/false,
3012  /*refinement_strategy=*/
3013  &::AdaptationStrategies::Refinement::
3014  preserve<dim, spacedim, unsigned int>,
3015  /*coarsening_strategy=*/
3016  [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3017  const std::vector<unsigned int> &children_fe_indices)
3018  -> unsigned int {
3019  return ::internal::hp::DoFHandlerImplementation::Implementation::
3020  determine_fe_from_children<dim, spacedim>(parent,
3021  children_fe_indices,
3022  fe_collection);
3023  });
3024 
3025  // If we work on a p::d::Triangulation, we have to transfer all
3026  // active FE indices since ownership of cells may change.
3027 
3028  // Gather all current active FE indices
3029  get_active_fe_indices(active_fe_index_transfer->active_fe_indices);
3030 
3031  // Attach to transfer object
3032  active_fe_index_transfer->cell_data_transfer->prepare_for_serialization(
3033  active_fe_index_transfer->active_fe_indices);
3034 #endif
3035 }
3036 
3037 
3038 
3039 template <int dim, int spacedim>
3040 void
3042 {
3043 #ifndef DEAL_II_WITH_P4EST
3044  Assert(false,
3045  ExcMessage(
3046  "You are attempting to use a functionality that is only available "
3047  "if deal.II was configured to use p4est, but cmake did not find a "
3048  "valid p4est library."));
3049 #else
3050  // the implementation below requires a p:d:T currently
3051  Assert(
3053  &this->get_triangulation()) != nullptr),
3054  ExcNotImplemented());
3055 
3056  Assert(active_fe_index_transfer == nullptr, ExcInternalError());
3057 
3058  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3059 
3060  // Create transfer object and attach to it.
3061  const auto *distributed_tria =
3063  &this->get_triangulation());
3064 
3065  active_fe_index_transfer->cell_data_transfer = std::make_unique<
3066  parallel::distributed::
3067  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3068  *distributed_tria,
3069  /*transfer_variable_size_data=*/false,
3070  /*refinement_strategy=*/
3071  &::AdaptationStrategies::Refinement::
3072  preserve<dim, spacedim, unsigned int>,
3073  /*coarsening_strategy=*/
3074  [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3075  const std::vector<unsigned int> &children_fe_indices)
3076  -> unsigned int {
3077  return ::internal::hp::DoFHandlerImplementation::Implementation::
3078  determine_fe_from_children<dim, spacedim>(parent,
3079  children_fe_indices,
3080  fe_collection);
3081  });
3082 
3083  // Unpack active FE indices.
3084  active_fe_index_transfer->active_fe_indices.resize(
3085  get_triangulation().n_active_cells(), numbers::invalid_unsigned_int);
3086  active_fe_index_transfer->cell_data_transfer->deserialize(
3087  active_fe_index_transfer->active_fe_indices);
3088 
3089  // Update all locally owned active FE indices.
3090  set_active_fe_indices(active_fe_index_transfer->active_fe_indices);
3091 
3092  // Update active FE indices on ghost cells.
3095 
3096  // Free memory.
3097  active_fe_index_transfer.reset();
3098 #endif
3099 }
3100 
3101 
3102 
3103 template <int dim, int spacedim>
3105  : coarsest_level(numbers::invalid_unsigned_int)
3106  , finest_level(0)
3107 {}
3108 
3109 
3110 
3111 template <int dim, int spacedim>
3112 void
3114  const unsigned int cl,
3115  const unsigned int fl,
3116  const unsigned int dofs_per_vertex)
3117 {
3118  coarsest_level = cl;
3119  finest_level = fl;
3120 
3121  if (coarsest_level <= finest_level)
3122  {
3123  const unsigned int n_levels = finest_level - coarsest_level + 1;
3124  const unsigned int n_indices = n_levels * dofs_per_vertex;
3125 
3126  indices = std::make_unique<types::global_dof_index[]>(n_indices);
3127  std::fill(indices.get(),
3128  indices.get() + n_indices,
3130  }
3131  else
3132  indices.reset();
3133 }
3134 
3135 
3136 
3137 template <int dim, int spacedim>
3138 unsigned int
3140 {
3141  return coarsest_level;
3142 }
3143 
3144 
3145 
3146 template <int dim, int spacedim>
3147 unsigned int
3149 {
3150  return finest_level;
3151 }
3152 
3153 /*-------------- Explicit Instantiations -------------------------------*/
3154 #include "dof_handler.inst"
3155 
3156 
3157 
void init(const unsigned int coarsest_level, const unsigned int finest_level, const unsigned int dofs_per_vertex)
unsigned int get_finest_level() const
unsigned int get_coarsest_level() const
cell_iterator end() const
void pre_distributed_transfer_action()
void post_transfer_action()
virtual std::size_t memory_consumption() const
std::vector< std::vector< active_fe_index_type > > hp_cell_active_fe_indices
Definition: dof_handler.h:1505
unsigned int max_couplings_between_dofs() const
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition: dof_handler.h:1525
std::unique_ptr< ActiveFEIndexTransfer > active_fe_index_transfer
Definition: dof_handler.h:1537
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1429
void create_active_fe_table()
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:1422
void set_active_fe_indices(const std::vector< unsigned int > &active_fe_indices)
level_cell_iterator end_mg() const
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
level_cell_iterator begin_mg(const unsigned int level=0) const
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
types::global_dof_index n_boundary_dofs() const
void clear_mg_space()
void connect_to_triangulation_signals()
std::vector< std::vector< types::global_dof_index > > cell_dof_cache_indices
Definition: dof_handler.h:1460
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1518
std::vector< std::vector< active_fe_index_type > > hp_cell_future_fe_indices
Definition: dof_handler.h:1512
void post_distributed_transfer_action()
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition: dof_handler.h:1474
const Triangulation< dim, spacedim > & get_triangulation() const
void pre_transfer_action()
void clear_space()
std::array< std::vector< active_fe_index_type >, dim+1 > hp_object_fe_indices
Definition: dof_handler.h:1493
void reinit(const Triangulation< dim, spacedim > &tria)
std::vector< std::array< std::vector< offset_type >, dim+1 > > object_dof_ptr
Definition: dof_handler.h:1485
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
Definition: dof_handler.h:1531
active_cell_iterator begin_active(const unsigned int level=0) const
void distribute_mg_dofs()
active_cell_iterator end_active(const unsigned int level) const
bool hp_capability_enabled
Definition: dof_handler.h:1416
void initialize_local_block_info()
types::global_dof_index n_dofs() const
void update_active_fe_table()
std::vector< std::vector< offset_type > > cell_dof_cache_ptr
Definition: dof_handler.h:1466
void prepare_for_serialization_of_active_fe_indices()
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
cell_iterator begin(const unsigned int level=0) const
std::array< std::vector< offset_type >, dim+1 > hp_object_fe_ptr
Definition: dof_handler.h:1498
void clear()
void get_active_fe_indices(std::vector< unsigned int > &active_fe_indices) const
unsigned int max_couplings_between_boundary_dofs() const
void setup_policy()
void deserialize_active_fe_indices()
virtual ~DoFHandler() override
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
bool is_element(const size_type index) const
Definition: index_set.h:1734
virtual const MeshSmoothing & get_mesh_smoothing() const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_raw_lines() const
unsigned int n_raw_faces() const
unsigned int n_levels() const
cell_iterator end() const
unsigned int n_raw_cells(const unsigned int level) const
unsigned int max_adjacent_cells() const
bool vertex_used(const unsigned int index) const
active_cell_iterator end_active(const unsigned int level) const
unsigned int n_raw_quads() const
unsigned int n_cells() const
Signals signals
Definition: tria.h:2448
unsigned int n_vertices() const
unsigned int n_raw_hexs(const unsigned int level) const
unsigned int size() const
Definition: collection.h:264
unsigned int max_dofs_per_line() const
unsigned int max_dofs_per_hex() const
unsigned int max_dofs_per_vertex() const
unsigned int max_dofs_per_quad() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
unsigned int level
Definition: grid_out.cc:4608
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< level_cell_iterator > mg_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators() const
IteratorRange< cell_iterator > cell_iterators() const
IteratorRange< level_cell_iterator > mg_cell_iterators_on_level(const unsigned int level) const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNoFESelected()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcNoDominatedFiniteElementOnChildren()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
Task< RT > new_task(const std::function< RT()> &function)
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
Definition: dof_tools.cc:1370
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
@ valid
Iterator points to a valid object.
static const char T
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1483
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1648
Definition: hp.h:118
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13794
unsigned int dominated_future_fe_on_children(const typename DoFHandler< dim, spacedim >::cell_iterator &parent)
void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
std::string policy_to_string(const ::internal::DoFHandlerImplementation::Policy::PolicyBase< dim, spacedim > &policy)
Definition: dof_handler.cc:58
const types::boundary_id internal_face_boundary_id
Definition: types.h:260
static const unsigned int invalid_unsigned_int
Definition: types.h:201
const types::global_dof_index invalid_dof_index
Definition: types.h:216
unsigned int global_dof_index
Definition: types.h:76
unsigned int boundary_id
Definition: types.h:129
static void reserve_space_mg(DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:597
static void reserve_space(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:400
static void reserve_subentities(DoFHandler< dim, spacedim > &dof_handler, const unsigned int structdim, const unsigned int n_raw_entities, const T &cell_process)
Definition: dof_handler.cc:341
static unsigned int max_couplings_between_dofs(const DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:113
static void reset_to_empty_objects(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:268
static unsigned int max_couplings_between_dofs(const DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:232
static void reserve_space_mg(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:450
static unsigned int max_couplings_between_dofs(const DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:103
static void reserve_space_mg(DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:520
static void reserve_cells(DoFHandler< dim, spacedim > &dof_handler, const unsigned int n_inner_dofs_per_cell)
Definition: dof_handler.cc:292
static void collect_fe_indices_on_cells_to_be_refined(DoFHandler< dim, spacedim > &dof_handler)
static unsigned int determine_fe_from_children(const typename Triangulation< dim, spacedim >::cell_iterator &, const std::vector< unsigned int > &children_fe_indices, const ::hp::FECollection< dim, spacedim > &fe_collection)
static void reserve_space_cells(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:852
static unsigned int dominated_future_fe_on_children(const typename DoFHandler< dim, spacedim >::cell_iterator &parent)
static void distribute_fe_indices_on_refined_cells(DoFHandler< dim, spacedim > &dof_handler)
static void ensure_absence_of_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:699
static void reserve_space(::DoFHandler< 1, spacedim > &dof_handler)
static void reserve_space_faces(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:911
static void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
static void communicate_active_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
static void reserve_space_vertices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:719
static void reserve_space(::DoFHandler< 2, spacedim > &dof_handler)
static void reserve_space(::DoFHandler< 3, spacedim > &dof_handler)
const ::Triangulation< dim, spacedim > & tria