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