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