Reference documentation for deal.II version GIT 112f7bbc69 2023-05-28 21:25:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_handler.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/config.h>
17 
20 #include <deal.II/base/mpi.templates.h>
22 
23 #include <deal.II/distributed/cell_data_transfer.templates.h>
27 
30 
33 #include <deal.II/grid/tria.h>
37 
38 #include <algorithm>
39 #include <memory>
40 #include <set>
41 #include <unordered_set>
42 
44 
45 template <int dim, int spacedim>
46 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
48 
49 
50 namespace internal
51 {
52  template <int dim, int spacedim>
53  std::string
54  policy_to_string(const ::internal::DoFHandlerImplementation::Policy::
55  PolicyBase<dim, spacedim> &policy)
56  {
57  std::string policy_name;
58  if (dynamic_cast<const typename ::internal::DoFHandlerImplementation::
59  Policy::Sequential<dim, spacedim> *>(&policy) ||
60  dynamic_cast<const typename ::internal::DoFHandlerImplementation::
61  Policy::Sequential<dim, spacedim> *>(&policy))
62  policy_name = "Policy::Sequential<";
63  else if (dynamic_cast<
64  const typename ::internal::DoFHandlerImplementation::
65  Policy::ParallelDistributed<dim, spacedim> *>(&policy) ||
66  dynamic_cast<
67  const typename ::internal::DoFHandlerImplementation::
68  Policy::ParallelDistributed<dim, spacedim> *>(&policy))
69  policy_name = "Policy::ParallelDistributed<";
70  else if (dynamic_cast<
71  const typename ::internal::DoFHandlerImplementation::
72  Policy::ParallelShared<dim, spacedim> *>(&policy) ||
73  dynamic_cast<
74  const typename ::internal::DoFHandlerImplementation::
75  Policy::ParallelShared<dim, spacedim> *>(&policy))
76  policy_name = "Policy::ParallelShared<";
77  else
79  policy_name += Utilities::int_to_string(dim) + "," +
80  Utilities::int_to_string(spacedim) + ">";
81  return policy_name;
82  }
83 
84 
85  namespace DoFHandlerImplementation
86  {
92  {
97  template <int spacedim>
98  static unsigned int
100  {
101  return std::min(static_cast<types::global_dof_index>(
102  3 * dof_handler.fe_collection.max_dofs_per_vertex() +
103  2 * dof_handler.fe_collection.max_dofs_per_line()),
104  dof_handler.n_dofs());
105  }
106 
107  template <int spacedim>
108  static unsigned int
110  {
111  // get these numbers by drawing pictures
112  // and counting...
113  // example:
114  // | | |
115  // --x-----x--x--X--
116  // | | | |
117  // | x--x--x
118  // | | | |
119  // --x--x--*--x--x--
120  // | | | |
121  // x--x--x |
122  // | | | |
123  // --X--x--x-----x--
124  // | | |
125  // x = vertices connected with center vertex *;
126  // = total of 19
127  // (the X vertices are connected with * if
128  // the vertices adjacent to X are hanging
129  // nodes)
130  // count lines -> 28 (don't forget to count
131  // mother and children separately!)
132  types::global_dof_index max_couplings;
133  switch (dof_handler.tria->max_adjacent_cells())
134  {
135  case 4:
136  max_couplings =
137  19 * dof_handler.fe_collection.max_dofs_per_vertex() +
138  28 * dof_handler.fe_collection.max_dofs_per_line() +
139  8 * dof_handler.fe_collection.max_dofs_per_quad();
140  break;
141  case 5:
142  max_couplings =
143  21 * dof_handler.fe_collection.max_dofs_per_vertex() +
144  31 * dof_handler.fe_collection.max_dofs_per_line() +
145  9 * dof_handler.fe_collection.max_dofs_per_quad();
146  break;
147  case 6:
148  max_couplings =
149  28 * dof_handler.fe_collection.max_dofs_per_vertex() +
150  42 * dof_handler.fe_collection.max_dofs_per_line() +
151  12 * dof_handler.fe_collection.max_dofs_per_quad();
152  break;
153  case 7:
154  max_couplings =
155  30 * dof_handler.fe_collection.max_dofs_per_vertex() +
156  45 * dof_handler.fe_collection.max_dofs_per_line() +
157  13 * dof_handler.fe_collection.max_dofs_per_quad();
158  break;
159  case 8:
160  max_couplings =
161  37 * dof_handler.fe_collection.max_dofs_per_vertex() +
162  56 * dof_handler.fe_collection.max_dofs_per_line() +
163  16 * dof_handler.fe_collection.max_dofs_per_quad();
164  break;
165 
166  // the following numbers are not based on actual counting but by
167  // extrapolating the number sequences from the previous ones (for
168  // example, for n_dofs_per_vertex(), the sequence above is 19, 21,
169  // 28, 30, 37, and is continued as follows):
170  case 9:
171  max_couplings =
172  39 * dof_handler.fe_collection.max_dofs_per_vertex() +
173  59 * dof_handler.fe_collection.max_dofs_per_line() +
174  17 * dof_handler.fe_collection.max_dofs_per_quad();
175  break;
176  case 10:
177  max_couplings =
178  46 * dof_handler.fe_collection.max_dofs_per_vertex() +
179  70 * dof_handler.fe_collection.max_dofs_per_line() +
180  20 * dof_handler.fe_collection.max_dofs_per_quad();
181  break;
182  case 11:
183  max_couplings =
184  48 * dof_handler.fe_collection.max_dofs_per_vertex() +
185  73 * dof_handler.fe_collection.max_dofs_per_line() +
186  21 * dof_handler.fe_collection.max_dofs_per_quad();
187  break;
188  case 12:
189  max_couplings =
190  55 * dof_handler.fe_collection.max_dofs_per_vertex() +
191  84 * dof_handler.fe_collection.max_dofs_per_line() +
192  24 * dof_handler.fe_collection.max_dofs_per_quad();
193  break;
194  case 13:
195  max_couplings =
196  57 * dof_handler.fe_collection.max_dofs_per_vertex() +
197  87 * dof_handler.fe_collection.max_dofs_per_line() +
198  25 * dof_handler.fe_collection.max_dofs_per_quad();
199  break;
200  case 14:
201  max_couplings =
202  63 * dof_handler.fe_collection.max_dofs_per_vertex() +
203  98 * dof_handler.fe_collection.max_dofs_per_line() +
204  28 * dof_handler.fe_collection.max_dofs_per_quad();
205  break;
206  case 15:
207  max_couplings =
208  65 * dof_handler.fe_collection.max_dofs_per_vertex() +
209  103 * dof_handler.fe_collection.max_dofs_per_line() +
210  29 * dof_handler.fe_collection.max_dofs_per_quad();
211  break;
212  case 16:
213  max_couplings =
214  72 * dof_handler.fe_collection.max_dofs_per_vertex() +
215  114 * dof_handler.fe_collection.max_dofs_per_line() +
216  32 * dof_handler.fe_collection.max_dofs_per_quad();
217  break;
218 
219  default:
220  Assert(false, ExcNotImplemented());
221  max_couplings = 0;
222  }
223  return std::min(max_couplings, dof_handler.n_dofs());
224  }
225 
226  template <int spacedim>
227  static unsigned int
229  {
230  // TODO:[?] Invent significantly better estimates than the ones in this
231  // function
232 
233  // doing the same thing here is a rather complicated thing, compared
234  // to the 2d case, since it is hard to draw pictures with several
235  // refined hexahedra :-) so I presently only give a coarse
236  // estimate for the case that at most 8 hexes meet at each vertex
237  //
238  // can anyone give better estimate here?
239  const unsigned int max_adjacent_cells =
240  dof_handler.tria->max_adjacent_cells();
241 
242  types::global_dof_index max_couplings;
243  if (max_adjacent_cells <= 8)
244  max_couplings =
245  7 * 7 * 7 * dof_handler.fe_collection.max_dofs_per_vertex() +
246  7 * 6 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_line() +
247  9 * 4 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_quad() +
248  27 * dof_handler.fe_collection.max_dofs_per_hex();
249  else
250  {
251  Assert(false, ExcNotImplemented());
252  max_couplings = 0;
253  }
254 
255  return std::min(max_couplings, dof_handler.n_dofs());
256  }
257 
262  template <int dim, int spacedim>
263  static void
265  {
266  dof_handler.object_dof_indices.clear();
267  dof_handler.object_dof_indices.resize(dof_handler.tria->n_levels());
268  dof_handler.object_dof_indices.shrink_to_fit();
269 
270  dof_handler.object_dof_ptr.clear();
271  dof_handler.object_dof_ptr.resize(dof_handler.tria->n_levels());
272  dof_handler.object_dof_ptr.shrink_to_fit();
273  }
274 
278  template <int dim, int spacedim>
279  static void
281  const unsigned int n_inner_dofs_per_cell)
282  {
283  for (unsigned int i = 0; i < dof_handler.tria->n_levels(); ++i)
284  {
285  dof_handler.object_dof_ptr[i][dim].assign(
286  dof_handler.tria->n_raw_cells(i) + 1, 0);
287 
288  for (const auto &cell :
289  dof_handler.tria->cell_iterators_on_level(i))
290  if (cell->is_active() && !cell->is_artificial())
291  dof_handler.object_dof_ptr[i][dim][cell->index() + 1] =
292  n_inner_dofs_per_cell;
293 
294  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i); ++j)
295  dof_handler.object_dof_ptr[i][dim][j + 1] +=
296  dof_handler.object_dof_ptr[i][dim][j];
297 
298  dof_handler.object_dof_indices[i][dim].resize(
299  dof_handler.object_dof_ptr[i][dim].back(),
301  }
302  }
303 
308  template <int dim, int spacedim, typename T>
309  static void
311  const unsigned int structdim,
312  const unsigned int n_raw_entities,
313  const T & cell_process)
314  {
315  if (dof_handler.tria->n_cells() == 0)
316  return;
317 
318  dof_handler.object_dof_ptr[0][structdim].assign(n_raw_entities + 1, -1);
319  // determine for each entity the number of dofs
320  for (const auto &cell : dof_handler.tria->cell_iterators())
321  if (cell->is_active() && !cell->is_artificial())
322  cell_process(
323  cell,
324  [&](const unsigned int n_dofs_per_entity,
325  const unsigned int index) {
326  auto &n_dofs_per_entity_target =
327  dof_handler.object_dof_ptr[0][structdim][index + 1];
328 
329  // make sure that either the entity has not been visited or
330  // the entity has the same number of dofs assigned
331  Assert((n_dofs_per_entity_target ==
332  static_cast<
334  -1) ||
335  n_dofs_per_entity_target == n_dofs_per_entity),
337 
338  n_dofs_per_entity_target = n_dofs_per_entity;
339  });
340 
341  // convert the absolute numbers to CRS
342  dof_handler.object_dof_ptr[0][structdim][0] = 0;
343  for (unsigned int i = 1; i < n_raw_entities + 1; ++i)
344  {
345  if (dof_handler.object_dof_ptr[0][structdim][i] ==
346  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
347  -1))
348  dof_handler.object_dof_ptr[0][structdim][i] =
349  dof_handler.object_dof_ptr[0][structdim][i - 1];
350  else
351  dof_handler.object_dof_ptr[0][structdim][i] +=
352  dof_handler.object_dof_ptr[0][structdim][i - 1];
353  }
354 
355  // allocate memory for indices
356  dof_handler.object_dof_indices[0][structdim].resize(
357  dof_handler.object_dof_ptr[0][structdim].back(),
359  }
360 
367  template <int dim, int spacedim>
368  static void
370  {
371  reset_to_empty_objects(dof_handler);
372 
373  const auto &fe = dof_handler.get_fe();
374 
375  // cell
376  reserve_cells(dof_handler,
377  dim == 1 ? fe.n_dofs_per_line() :
378  (dim == 2 ? fe.n_dofs_per_quad(0) :
379  fe.n_dofs_per_hex()));
380 
381  // vertices
382  reserve_subentities(dof_handler,
383  0,
384  dof_handler.tria->n_vertices(),
385  [&](const auto &cell, const auto &process) {
386  for (const auto vertex_index :
387  cell->vertex_indices())
388  process(fe.n_dofs_per_vertex(),
389  cell->vertex_index(vertex_index));
390  });
391 
392  // lines
393  if (dim == 2 || dim == 3)
394  reserve_subentities(dof_handler,
395  1,
396  dof_handler.tria->n_raw_lines(),
397  [&](const auto &cell, const auto &process) {
398  for (const auto line_index :
399  cell->line_indices())
400  process(fe.n_dofs_per_line(),
401  cell->line(line_index)->index());
402  });
403 
404  // quads
405  if (dim == 3)
406  reserve_subentities(dof_handler,
407  2,
408  dof_handler.tria->n_raw_quads(),
409  [&](const auto &cell, const auto &process) {
410  for (const auto face_index :
411  cell->face_indices())
412  process(fe.n_dofs_per_quad(face_index),
413  cell->face(face_index)->index());
414  });
415  }
416 
417  template <int spacedim>
418  static void
420  {
421  Assert(dof_handler.get_triangulation().n_levels() > 0,
422  ExcMessage("Invalid triangulation"));
423  dof_handler.clear_mg_space();
424 
425  const ::Triangulation<1, spacedim> &tria =
426  dof_handler.get_triangulation();
427  const unsigned int dofs_per_line =
428  dof_handler.get_fe().n_dofs_per_line();
429  const unsigned int n_levels = tria.n_levels();
430 
431  for (unsigned int i = 0; i < n_levels; ++i)
432  {
433  dof_handler.mg_levels.emplace_back(
435  dof_handler.mg_levels.back()->dof_object.dofs =
436  std::vector<types::global_dof_index>(tria.n_raw_lines(i) *
437  dofs_per_line,
439  }
440 
441  const unsigned int n_vertices = tria.n_vertices();
442 
443  dof_handler.mg_vertex_dofs.resize(n_vertices);
444 
445  std::vector<unsigned int> max_level(n_vertices, 0);
446  std::vector<unsigned int> min_level(n_vertices, n_levels);
447 
448  for (typename ::Triangulation<1, spacedim>::cell_iterator cell =
449  tria.begin();
450  cell != tria.end();
451  ++cell)
452  {
453  const unsigned int level = cell->level();
454 
455  for (const auto vertex : cell->vertex_indices())
456  {
457  const unsigned int vertex_index = cell->vertex_index(vertex);
458 
459  if (min_level[vertex_index] > level)
460  min_level[vertex_index] = level;
461 
462  if (max_level[vertex_index] < level)
463  max_level[vertex_index] = level;
464  }
465  }
466 
467  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
468  if (tria.vertex_used(vertex))
469  {
470  Assert(min_level[vertex] < n_levels, ExcInternalError());
471  Assert(max_level[vertex] >= min_level[vertex],
472  ExcInternalError());
473  dof_handler.mg_vertex_dofs[vertex].init(
474  min_level[vertex],
475  max_level[vertex],
476  dof_handler.get_fe().n_dofs_per_vertex());
477  }
478 
479  else
480  {
481  Assert(min_level[vertex] == n_levels, ExcInternalError());
482  Assert(max_level[vertex] == 0, ExcInternalError());
483  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
484  }
485  }
486 
487  template <int spacedim>
488  static void
490  {
491  Assert(dof_handler.get_triangulation().n_levels() > 0,
492  ExcMessage("Invalid triangulation"));
493  dof_handler.clear_mg_space();
494 
495  const ::FiniteElement<2, spacedim> &fe = dof_handler.get_fe();
496  const ::Triangulation<2, spacedim> &tria =
497  dof_handler.get_triangulation();
498  const unsigned int n_levels = tria.n_levels();
499 
500  for (unsigned int i = 0; i < n_levels; ++i)
501  {
502  dof_handler.mg_levels.emplace_back(
503  std::make_unique<
505  dof_handler.mg_levels.back()->dof_object.dofs =
506  std::vector<types::global_dof_index>(
507  tria.n_raw_quads(i) *
508  fe.n_dofs_per_quad(0 /*note: in 2d there is only one quad*/),
510  }
511 
512  dof_handler.mg_faces =
513  std::make_unique<internal::DoFHandlerImplementation::DoFFaces<2>>();
514  dof_handler.mg_faces->lines.dofs =
515  std::vector<types::global_dof_index>(tria.n_raw_lines() *
516  fe.n_dofs_per_line(),
518 
519  const unsigned int n_vertices = tria.n_vertices();
520 
521  dof_handler.mg_vertex_dofs.resize(n_vertices);
522 
523  std::vector<unsigned int> max_level(n_vertices, 0);
524  std::vector<unsigned int> min_level(n_vertices, n_levels);
525 
526  for (typename ::Triangulation<2, spacedim>::cell_iterator cell =
527  tria.begin();
528  cell != tria.end();
529  ++cell)
530  {
531  const unsigned int level = cell->level();
532 
533  for (const auto vertex : cell->vertex_indices())
534  {
535  const unsigned int vertex_index = cell->vertex_index(vertex);
536 
537  if (min_level[vertex_index] > level)
538  min_level[vertex_index] = level;
539 
540  if (max_level[vertex_index] < level)
541  max_level[vertex_index] = level;
542  }
543  }
544 
545  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
546  if (tria.vertex_used(vertex))
547  {
548  Assert(min_level[vertex] < n_levels, ExcInternalError());
549  Assert(max_level[vertex] >= min_level[vertex],
550  ExcInternalError());
551  dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
552  max_level[vertex],
553  fe.n_dofs_per_vertex());
554  }
555 
556  else
557  {
558  Assert(min_level[vertex] == n_levels, ExcInternalError());
559  Assert(max_level[vertex] == 0, ExcInternalError());
560  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
561  }
562  }
563 
564  template <int spacedim>
565  static void
567  {
568  Assert(dof_handler.get_triangulation().n_levels() > 0,
569  ExcMessage("Invalid triangulation"));
570  dof_handler.clear_mg_space();
571 
572  const ::FiniteElement<3, spacedim> &fe = dof_handler.get_fe();
573  const ::Triangulation<3, spacedim> &tria =
574  dof_handler.get_triangulation();
575  const unsigned int n_levels = tria.n_levels();
576 
577  for (unsigned int i = 0; i < n_levels; ++i)
578  {
579  dof_handler.mg_levels.emplace_back(
580  std::make_unique<
582  dof_handler.mg_levels.back()->dof_object.dofs =
583  std::vector<types::global_dof_index>(tria.n_raw_hexs(i) *
584  fe.n_dofs_per_hex(),
586  }
587 
588  dof_handler.mg_faces =
589  std::make_unique<internal::DoFHandlerImplementation::DoFFaces<3>>();
590  dof_handler.mg_faces->lines.dofs =
591  std::vector<types::global_dof_index>(tria.n_raw_lines() *
592  fe.n_dofs_per_line(),
594 
595  // TODO: the implementation makes the assumption that all faces have the
596  // same number of dofs
597  AssertDimension(fe.n_unique_faces(), 1);
598  dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index>(
599  tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/),
601 
602  const unsigned int n_vertices = tria.n_vertices();
603 
604  dof_handler.mg_vertex_dofs.resize(n_vertices);
605 
606  std::vector<unsigned int> max_level(n_vertices, 0);
607  std::vector<unsigned int> min_level(n_vertices, n_levels);
608 
609  for (typename ::Triangulation<3, spacedim>::cell_iterator cell =
610  tria.begin();
611  cell != tria.end();
612  ++cell)
613  {
614  const unsigned int level = cell->level();
615 
616  for (const auto vertex : cell->vertex_indices())
617  {
618  const unsigned int vertex_index = cell->vertex_index(vertex);
619 
620  if (min_level[vertex_index] > level)
621  min_level[vertex_index] = level;
622 
623  if (max_level[vertex_index] < level)
624  max_level[vertex_index] = level;
625  }
626  }
627 
628  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
629  if (tria.vertex_used(vertex))
630  {
631  Assert(min_level[vertex] < n_levels, ExcInternalError());
632  Assert(max_level[vertex] >= min_level[vertex],
633  ExcInternalError());
634  dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
635  max_level[vertex],
636  fe.n_dofs_per_vertex());
637  }
638 
639  else
640  {
641  Assert(min_level[vertex] == n_levels, ExcInternalError());
642  Assert(max_level[vertex] == 0, ExcInternalError());
643  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
644  }
645  }
646  };
647  } // namespace DoFHandlerImplementation
648 
649 
650 
651  namespace hp
652  {
653  namespace DoFHandlerImplementation
654  {
660  {
666  template <int dim, int spacedim>
667  static void
669  DoFHandler<dim, spacedim> &dof_handler)
670  {
671  (void)dof_handler;
672  for (const auto &cell : dof_handler.active_cell_iterators())
673  if (cell->is_locally_owned())
674  Assert(
675  !cell->future_fe_index_set(),
676  ExcMessage(
677  "There shouldn't be any cells flagged for p-adaptation when partitioning."));
678  }
679 
680 
681 
686  template <int dim, int spacedim>
687  static void
689  {
690  // The final step in all of the reserve_space() functions is to set
691  // up vertex dof information. since vertices are sequentially
692  // numbered, what we do first is to set up an array in which
693  // we record whether a vertex is associated with any of the
694  // given fe's, by setting a bit. in a later step, we then
695  // actually allocate memory for the required dofs
696  //
697  // in the following, we only need to consider vertices that are
698  // adjacent to either a locally owned or a ghost cell; we never
699  // store anything on vertices that are only surrounded by
700  // artificial cells. so figure out that subset of vertices
701  // first
702  std::vector<bool> locally_used_vertices(
703  dof_handler.tria->n_vertices(), false);
704  for (const auto &cell : dof_handler.active_cell_iterators())
705  if (!cell->is_artificial())
706  for (const auto v : cell->vertex_indices())
707  locally_used_vertices[cell->vertex_index(v)] = true;
708 
709  std::vector<std::vector<bool>> vertex_fe_association(
710  dof_handler.fe_collection.size(),
711  std::vector<bool>(dof_handler.tria->n_vertices(), false));
712 
713  for (const auto &cell : dof_handler.active_cell_iterators())
714  if (!cell->is_artificial())
715  for (const auto v : cell->vertex_indices())
716  vertex_fe_association[cell->active_fe_index()]
717  [cell->vertex_index(v)] = true;
718 
719  // in debug mode, make sure that each vertex is associated
720  // with at least one FE (note that except for unused
721  // vertices, all vertices are actually active). this is of
722  // course only true for vertices that are part of either
723  // ghost or locally owned cells
724 #ifdef DEBUG
725  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
726  if (locally_used_vertices[v] == true)
727  if (dof_handler.tria->vertex_used(v) == true)
728  {
729  unsigned int fe = 0;
730  for (; fe < dof_handler.fe_collection.size(); ++fe)
731  if (vertex_fe_association[fe][v] == true)
732  break;
733  Assert(fe != dof_handler.fe_collection.size(),
734  ExcInternalError());
735  }
736 #endif
737 
738  const unsigned int d = 0;
739  const unsigned int l = 0;
740 
741  dof_handler.hp_object_fe_ptr[d].clear();
742  dof_handler.hp_object_fe_indices[d].clear();
743  dof_handler.object_dof_ptr[l][d].clear();
744  dof_handler.object_dof_indices[l][d].clear();
745 
746  dof_handler.hp_object_fe_ptr[d].reserve(
747  dof_handler.tria->n_vertices() + 1);
748 
749  unsigned int vertex_slots_needed = 0;
750  unsigned int fe_slots_needed = 0;
751 
752  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
753  {
754  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
755 
756  if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
757  {
758  for (unsigned int fe = 0;
759  fe < dof_handler.fe_collection.size();
760  ++fe)
761  if (vertex_fe_association[fe][v] == true)
762  {
763  fe_slots_needed++;
764  vertex_slots_needed +=
765  dof_handler.get_fe(fe).n_dofs_per_vertex();
766  }
767  }
768  }
769 
770  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
771 
772  dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
773  dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
774 
775  dof_handler.object_dof_indices[l][d].reserve(vertex_slots_needed);
776 
777  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
778  if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
779  {
780  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
781  ++fe)
782  if (vertex_fe_association[fe][v] == true)
783  {
784  dof_handler.hp_object_fe_indices[d].push_back(fe);
785  dof_handler.object_dof_ptr[l][d].push_back(
786  dof_handler.object_dof_indices[l][d].size());
787 
788  for (unsigned int i = 0;
789  i < dof_handler.get_fe(fe).n_dofs_per_vertex();
790  i++)
791  dof_handler.object_dof_indices[l][d].push_back(
793  }
794  }
795 
796 
797  dof_handler.object_dof_ptr[l][d].push_back(
798  dof_handler.object_dof_indices[l][d].size());
799 
800  AssertDimension(vertex_slots_needed,
801  dof_handler.object_dof_indices[l][d].size());
802  AssertDimension(fe_slots_needed,
803  dof_handler.hp_object_fe_indices[d].size());
804  AssertDimension(fe_slots_needed + 1,
805  dof_handler.object_dof_ptr[l][d].size());
806  AssertDimension(dof_handler.tria->n_vertices() + 1,
807  dof_handler.hp_object_fe_ptr[d].size());
808 
809  dof_handler.object_dof_indices[l][d].assign(
810  vertex_slots_needed, numbers::invalid_dof_index);
811  }
812 
813 
814 
819  template <int dim, int spacedim>
820  static void
822  {
823  (void)dof_handler;
824  // count how much space we need on each level for the cell
825  // dofs and set the dof_*_offsets data. initially set the
826  // latter to an invalid index, and only later set it to
827  // something reasonable for active dof_handler.cells
828  //
829  // note that for dof_handler.cells, the situation is simpler
830  // than for other (lower dimensional) objects since exactly
831  // one finite element is used for it
832  for (unsigned int level = 0; level < dof_handler.tria->n_levels();
833  ++level)
834  {
835  dof_handler.object_dof_ptr[level][dim] =
836  std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
837  dof_handler.tria->n_raw_cells(level),
838  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
839  -1));
840 
841  types::global_dof_index next_free_dof = 0;
842  for (auto cell :
844  if (cell->is_active() && !cell->is_artificial())
845  {
846  dof_handler.object_dof_ptr[level][dim][cell->index()] =
847  next_free_dof;
848  next_free_dof +=
849  cell->get_fe().template n_dofs_per_object<dim>();
850  }
851 
852  dof_handler.object_dof_indices[level][dim] =
853  std::vector<types::global_dof_index>(
854  next_free_dof, numbers::invalid_dof_index);
855  }
856  }
857 
858 
859 
864  template <int dim, int spacedim>
865  static void
867  {
868  // FACE DOFS
869  //
870  // Count face dofs, then allocate as much space
871  // as we need and prime the linked list for faces (see the
872  // description in hp::DoFLevel) with the indices we will
873  // need. Note that our task is more complicated than for the
874  // cell case above since two adjacent cells may have different
875  // active FE indices, in which case we need to allocate
876  // *two* sets of face dofs for the same face. But they don't
877  // *have* to be different, and so we need to prepare for this
878  // as well.
879  //
880  // The way we do things is that we loop over all active cells (these
881  // are the only ones that have DoFs anyway) and all their faces. We
882  // note in the vector face_touched whether we have previously
883  // visited a face and if so skip it
884  {
885  std::vector<bool> face_touched(dim == 2 ?
886  dof_handler.tria->n_raw_lines() :
887  dof_handler.tria->n_raw_quads());
888 
889  const unsigned int d = dim - 1;
890  const unsigned int l = 0;
891 
892  dof_handler.hp_object_fe_ptr[d].clear();
893  dof_handler.hp_object_fe_indices[d].clear();
894  dof_handler.object_dof_ptr[l][d].clear();
895  dof_handler.object_dof_indices[l][d].clear();
896 
897  dof_handler.hp_object_fe_ptr[d].resize(
898  dof_handler.tria->n_raw_faces() + 1);
899 
900  // An array to hold how many slots (see the hp::DoFLevel
901  // class) we will have to store on each level
902  unsigned int n_face_slots = 0;
903 
904  for (const auto &cell : dof_handler.active_cell_iterators())
905  if (!cell->is_artificial())
906  for (const auto face : cell->face_indices())
907  if (!face_touched[cell->face(face)->index()])
908  {
909  unsigned int fe_slots_needed = 0;
910 
911  if (cell->at_boundary(face) ||
912  cell->face(face)->has_children() ||
913  cell->neighbor_is_coarser(face) ||
914  (!cell->at_boundary(face) &&
915  cell->neighbor(face)->is_artificial()) ||
916  (!cell->at_boundary(face) &&
917  !cell->neighbor(face)->is_artificial() &&
918  (cell->active_fe_index() ==
919  cell->neighbor(face)->active_fe_index())))
920  {
921  fe_slots_needed = 1;
922  n_face_slots +=
923  dof_handler.get_fe(cell->active_fe_index())
924  .template n_dofs_per_object<dim - 1>(face);
925  }
926  else
927  {
928  fe_slots_needed = 2;
929  n_face_slots +=
930  dof_handler.get_fe(cell->active_fe_index())
931  .template n_dofs_per_object<dim - 1>(face) +
932  dof_handler
933  .get_fe(cell->neighbor(face)->active_fe_index())
934  .template n_dofs_per_object<dim - 1>(
935  cell->neighbor_face_no(face));
936  }
937 
938  // mark this face as visited
939  face_touched[cell->face(face)->index()] = true;
940 
941  dof_handler
942  .hp_object_fe_ptr[d][cell->face(face)->index() + 1] =
943  fe_slots_needed;
944  }
945 
946  for (unsigned int i = 1; i < dof_handler.hp_object_fe_ptr[d].size();
947  i++)
948  dof_handler.hp_object_fe_ptr[d][i] +=
949  dof_handler.hp_object_fe_ptr[d][i - 1];
950 
951 
952  dof_handler.hp_object_fe_indices[d].resize(
953  dof_handler.hp_object_fe_ptr[d].back());
954  dof_handler.object_dof_ptr[l][d].resize(
955  dof_handler.hp_object_fe_ptr[d].back() + 1);
956 
957  dof_handler.object_dof_indices[l][d].reserve(n_face_slots);
958 
959 
960  // With the memory now allocated, loop over the
961  // dof_handler cells again and prime the _offset values as
962  // well as the fe_index fields
963  face_touched = std::vector<bool>(face_touched.size());
964 
965  for (const auto &cell : dof_handler.active_cell_iterators())
966  if (!cell->is_artificial())
967  for (const auto face : cell->face_indices())
968  if (!face_touched[cell->face(face)->index()])
969  {
970  // Same decision tree as before
971  if (cell->at_boundary(face) ||
972  cell->face(face)->has_children() ||
973  cell->neighbor_is_coarser(face) ||
974  (!cell->at_boundary(face) &&
975  cell->neighbor(face)->is_artificial()) ||
976  (!cell->at_boundary(face) &&
977  !cell->neighbor(face)->is_artificial() &&
978  (cell->active_fe_index() ==
979  cell->neighbor(face)->active_fe_index())))
980  {
981  const types::fe_index fe = cell->active_fe_index();
982  const unsigned int n_dofs =
983  dof_handler.get_fe(fe)
984  .template n_dofs_per_object<dim - 1>(face);
985  const unsigned int offset =
986  dof_handler
987  .hp_object_fe_ptr[d][cell->face(face)->index()];
988 
989  dof_handler.hp_object_fe_indices[d][offset] = fe;
990  dof_handler.object_dof_ptr[l][d][offset + 1] = n_dofs;
991 
992  for (unsigned int i = 0; i < n_dofs; ++i)
993  dof_handler.object_dof_indices[l][d].push_back(
995  }
996  else
997  {
998  types::fe_index fe_1 = cell->active_fe_index();
999  unsigned int face_no_1 = face;
1000  types::fe_index fe_2 =
1001  cell->neighbor(face)->active_fe_index();
1002  unsigned int face_no_2 = cell->neighbor_face_no(face);
1003 
1004  if (fe_2 < fe_1)
1005  {
1006  std::swap(fe_1, fe_2);
1007  std::swap(face_no_1, face_no_2);
1008  }
1009 
1010  const unsigned int n_dofs_1 =
1011  dof_handler.get_fe(fe_1)
1012  .template n_dofs_per_object<dim - 1>(face_no_1);
1013 
1014  const unsigned int n_dofs_2 =
1015  dof_handler.get_fe(fe_2)
1016  .template n_dofs_per_object<dim - 1>(face_no_2);
1017 
1018  const unsigned int offset =
1019  dof_handler
1020  .hp_object_fe_ptr[d][cell->face(face)->index()];
1021 
1022  dof_handler.hp_object_fe_indices[d].push_back(
1023  cell->active_fe_index());
1024  dof_handler.object_dof_ptr[l][d].push_back(
1025  dof_handler.object_dof_indices[l][d].size());
1026 
1027  dof_handler.hp_object_fe_indices[d][offset + 0] =
1028  fe_1;
1029  dof_handler.hp_object_fe_indices[d][offset + 1] =
1030  fe_2;
1031  dof_handler.object_dof_ptr[l][d][offset + 1] =
1032  n_dofs_1;
1033  dof_handler.object_dof_ptr[l][d][offset + 2] =
1034  n_dofs_2;
1035 
1036 
1037  for (unsigned int i = 0; i < n_dofs_1 + n_dofs_2; ++i)
1038  dof_handler.object_dof_indices[l][d].push_back(
1040  }
1041 
1042  // mark this face as visited
1043  face_touched[cell->face(face)->index()] = true;
1044  }
1045 
1046  for (unsigned int i = 1;
1047  i < dof_handler.object_dof_ptr[l][d].size();
1048  i++)
1049  dof_handler.object_dof_ptr[l][d][i] +=
1050  dof_handler.object_dof_ptr[l][d][i - 1];
1051  }
1052  }
1053 
1054 
1055 
1062  template <int spacedim>
1063  static void
1065  {
1066  Assert(dof_handler.fe_collection.size() > 0,
1068  Assert(dof_handler.tria->n_levels() > 0,
1069  ExcMessage("The current Triangulation must not be empty."));
1070  Assert(dof_handler.tria->n_levels() ==
1071  dof_handler.hp_cell_future_fe_indices.size(),
1072  ExcInternalError());
1073 
1075  reset_to_empty_objects(dof_handler);
1076 
1077  Threads::TaskGroup<> tasks;
1078  tasks +=
1079  Threads::new_task(&reserve_space_cells<1, spacedim>, dof_handler);
1080  tasks += Threads::new_task(&reserve_space_vertices<1, spacedim>,
1081  dof_handler);
1082  tasks.join_all();
1083  }
1084 
1085 
1086 
1087  template <int spacedim>
1088  static void
1090  {
1091  Assert(dof_handler.fe_collection.size() > 0,
1093  Assert(dof_handler.tria->n_levels() > 0,
1094  ExcMessage("The current Triangulation must not be empty."));
1095  Assert(dof_handler.tria->n_levels() ==
1096  dof_handler.hp_cell_future_fe_indices.size(),
1097  ExcInternalError());
1098 
1100  reset_to_empty_objects(dof_handler);
1101 
1102  Threads::TaskGroup<> tasks;
1103  tasks +=
1104  Threads::new_task(&reserve_space_cells<2, spacedim>, dof_handler);
1105  tasks +=
1106  Threads::new_task(&reserve_space_faces<2, spacedim>, dof_handler);
1107  tasks += Threads::new_task(&reserve_space_vertices<2, spacedim>,
1108  dof_handler);
1109  tasks.join_all();
1110  }
1111 
1112 
1113 
1114  template <int spacedim>
1115  static void
1117  {
1118  Assert(dof_handler.fe_collection.size() > 0,
1120  Assert(dof_handler.tria->n_levels() > 0,
1121  ExcMessage("The current Triangulation must not be empty."));
1122  Assert(dof_handler.tria->n_levels() ==
1123  dof_handler.hp_cell_future_fe_indices.size(),
1124  ExcInternalError());
1125 
1127  reset_to_empty_objects(dof_handler);
1128 
1129  Threads::TaskGroup<> tasks;
1130  tasks +=
1131  Threads::new_task(&reserve_space_cells<3, spacedim>, dof_handler);
1132  tasks +=
1133  Threads::new_task(&reserve_space_faces<3, spacedim>, dof_handler);
1134  tasks += Threads::new_task(&reserve_space_vertices<3, spacedim>,
1135  dof_handler);
1136 
1137  // While the tasks above are running, we can turn to line dofs
1138 
1139  // the situation here is pretty much like with vertices:
1140  // there can be an arbitrary number of finite elements
1141  // associated with each line.
1142  //
1143  // the algorithm we use is somewhat similar to what we do in
1144  // reserve_space_vertices()
1145  {
1146  // what we do first is to set up an array in which we
1147  // record whether a line is associated with any of the
1148  // given fe's, by setting a bit. in a later step, we
1149  // then actually allocate memory for the required dofs
1150  std::vector<std::vector<bool>> line_fe_association(
1151  dof_handler.fe_collection.size(),
1152  std::vector<bool>(dof_handler.tria->n_raw_lines(), false));
1153 
1154  for (const auto &cell : dof_handler.active_cell_iterators())
1155  if (!cell->is_artificial())
1156  for (const auto l : cell->line_indices())
1157  line_fe_association[cell->active_fe_index()]
1158  [cell->line_index(l)] = true;
1159 
1160  // first check which of the lines is used at all,
1161  // i.e. is associated with a finite element. we do this
1162  // since not all lines may actually be used, in which
1163  // case we do not have to allocate any memory at all
1164  std::vector<bool> line_is_used(dof_handler.tria->n_raw_lines(),
1165  false);
1166  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1167  ++line)
1168  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
1169  ++fe)
1170  if (line_fe_association[fe][line] == true)
1171  {
1172  line_is_used[line] = true;
1173  break;
1174  }
1175 
1176 
1177 
1178  const unsigned int d = 1;
1179  const unsigned int l = 0;
1180 
1181  dof_handler.hp_object_fe_ptr[d].clear();
1182  dof_handler.hp_object_fe_indices[d].clear();
1183  dof_handler.object_dof_ptr[l][d].clear();
1184  dof_handler.object_dof_indices[l][d].clear();
1185 
1186  dof_handler.hp_object_fe_ptr[d].reserve(
1187  dof_handler.tria->n_raw_lines() + 1);
1188 
1189  unsigned int line_slots_needed = 0;
1190  unsigned int fe_slots_needed = 0;
1191 
1192  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1193  ++line)
1194  {
1195  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1196 
1197  if (line_is_used[line] == true)
1198  {
1199  for (unsigned int fe = 0;
1200  fe < dof_handler.fe_collection.size();
1201  ++fe)
1202  if (line_fe_association[fe][line] == true)
1203  {
1204  fe_slots_needed++;
1205  line_slots_needed +=
1206  dof_handler.get_fe(fe).n_dofs_per_line();
1207  }
1208  }
1209  }
1210 
1211  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1212 
1213  // make sure that all entries have been set
1214  AssertDimension(dof_handler.hp_object_fe_ptr[d].size(),
1215  dof_handler.tria->n_raw_lines() + 1);
1216 
1217  dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
1218  dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
1219 
1220  dof_handler.object_dof_indices[l][d].reserve(line_slots_needed);
1221 
1222  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1223  ++line)
1224  if (line_is_used[line] == true)
1225  {
1226  for (unsigned int fe = 0;
1227  fe < dof_handler.fe_collection.size();
1228  ++fe)
1229  if (line_fe_association[fe][line] == true)
1230  {
1231  dof_handler.hp_object_fe_indices[d].push_back(fe);
1232  dof_handler.object_dof_ptr[l][d].push_back(
1233  dof_handler.object_dof_indices[l][d].size());
1234 
1235  for (unsigned int i = 0;
1236  i < dof_handler.get_fe(fe).n_dofs_per_line();
1237  i++)
1238  dof_handler.object_dof_indices[l][d].push_back(
1240  }
1241  }
1242 
1243  dof_handler.object_dof_ptr[l][d].push_back(
1244  dof_handler.object_dof_indices[l][d].size());
1245 
1246  // make sure that all entries have been set
1247  AssertDimension(dof_handler.hp_object_fe_indices[d].size(),
1248  fe_slots_needed);
1249  AssertDimension(dof_handler.object_dof_ptr[l][d].size(),
1250  fe_slots_needed + 1);
1251  AssertDimension(dof_handler.object_dof_indices[l][d].size(),
1252  line_slots_needed);
1253  }
1254 
1255  // Ensure that everything is done at this point.
1256  tasks.join_all();
1257  }
1258 
1259 
1260 
1272  template <int dim, int spacedim>
1273  static void
1275  {
1276  Assert(
1277  dof_handler.hp_capability_enabled == true,
1279 
1280  if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1281  dynamic_cast<
1282  const ::parallel::shared::Triangulation<dim, spacedim>
1283  *>(&dof_handler.get_triangulation()))
1284  {
1285  // we have a shared triangulation. in this case, every processor
1286  // knows about all cells, but every processor only has knowledge
1287  // about the active FE index on the cells it owns.
1288  //
1289  // we can create a complete set of active FE indices by letting
1290  // every processor create a vector of indices for all cells,
1291  // filling only those on the cells it owns and setting the indices
1292  // on the other cells to zero. then we add all of these vectors
1293  // up, and because every vector entry has exactly one processor
1294  // that owns it, the sum is correct
1295  std::vector<types::fe_index> active_fe_indices(
1296  tr->n_active_cells(), 0u);
1297  for (const auto &cell : dof_handler.active_cell_iterators())
1298  if (cell->is_locally_owned())
1299  active_fe_indices[cell->active_cell_index()] =
1300  cell->active_fe_index();
1301 
1302  Utilities::MPI::sum(active_fe_indices,
1303  tr->get_communicator(),
1304  active_fe_indices);
1305 
1306  // now go back and fill the active FE index on all other
1307  // cells. we would like to call cell->set_active_fe_index(),
1308  // but that function does not allow setting these indices on
1309  // non-locally_owned cells. so we have to work around the
1310  // issue a little bit by accessing the underlying data
1311  // structures directly
1312  for (const auto &cell : dof_handler.active_cell_iterators())
1313  if (!cell->is_locally_owned())
1314  dof_handler
1315  .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1316  active_fe_indices[cell->active_cell_index()];
1317  }
1318  else if (const ::parallel::
1319  DistributedTriangulationBase<dim, spacedim> *tr =
1320  dynamic_cast<
1321  const ::parallel::
1322  DistributedTriangulationBase<dim, spacedim> *>(
1323  &dof_handler.get_triangulation()))
1324  {
1325  // For completely distributed meshes, use the function that is
1326  // able to move data from locally owned cells on one processor to
1327  // the corresponding ghost cells on others. To this end, we need
1328  // to have functions that can pack and unpack the data we want to
1329  // transport -- namely, the single unsigned int active_fe_index
1330  // objects
1331  auto pack =
1332  [](
1334  &cell) -> types::fe_index {
1335  return cell->active_fe_index();
1336  };
1337 
1338  auto unpack =
1339  [&dof_handler](
1341  & cell,
1342  const types::fe_index active_fe_index) -> void {
1343  // we would like to say
1344  // cell->set_active_fe_index(active_fe_index);
1345  // but this is not allowed on cells that are not
1346  // locally owned, and we are on a ghost cell
1347  dof_handler
1348  .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1349  active_fe_index;
1350  };
1351 
1354  DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1355  }
1356  else
1357  {
1358  // a sequential triangulation. there is nothing we need to do here
1359  Assert(
1360  (dynamic_cast<
1361  const ::parallel::TriangulationBase<dim, spacedim> *>(
1362  &dof_handler.get_triangulation()) == nullptr),
1363  ExcInternalError());
1364  }
1365  }
1366 
1367 
1368 
1382  template <int dim, int spacedim>
1383  static void
1385  {
1386  Assert(
1387  dof_handler.hp_capability_enabled == true,
1389 
1390  if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1391  dynamic_cast<
1392  const ::parallel::shared::Triangulation<dim, spacedim>
1393  *>(&dof_handler.get_triangulation()))
1394  {
1395  std::vector<types::fe_index> future_fe_indices(
1396  tr->n_active_cells(), 0u);
1397  for (const auto &cell : dof_handler.active_cell_iterators() |
1399  future_fe_indices[cell->active_cell_index()] =
1400  dof_handler
1401  .hp_cell_future_fe_indices[cell->level()][cell->index()];
1402 
1403  Utilities::MPI::sum(future_fe_indices,
1404  tr->get_communicator(),
1405  future_fe_indices);
1406 
1407  for (const auto &cell : dof_handler.active_cell_iterators())
1408  if (!cell->is_locally_owned())
1409  dof_handler
1410  .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1411  future_fe_indices[cell->active_cell_index()];
1412  }
1413  else if (const ::parallel::
1414  DistributedTriangulationBase<dim, spacedim> *tr =
1415  dynamic_cast<
1416  const ::parallel::
1417  DistributedTriangulationBase<dim, spacedim> *>(
1418  &dof_handler.get_triangulation()))
1419  {
1420  auto pack =
1421  [&dof_handler](
1423  &cell) -> types::fe_index {
1424  return dof_handler
1425  .hp_cell_future_fe_indices[cell->level()][cell->index()];
1426  };
1427 
1428  auto unpack =
1429  [&dof_handler](
1431  & cell,
1432  const types::fe_index future_fe_index) -> void {
1433  dof_handler
1434  .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1435  future_fe_index;
1436  };
1437 
1440  DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1441  }
1442  else
1443  {
1444  Assert(
1445  (dynamic_cast<
1446  const ::parallel::TriangulationBase<dim, spacedim> *>(
1447  &dof_handler.get_triangulation()) == nullptr),
1448  ExcInternalError());
1449  }
1450  }
1451 
1452 
1453 
1474  template <int dim, int spacedim>
1475  static void
1477  DoFHandler<dim, spacedim> &dof_handler)
1478  {
1479  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1480 
1481  for (const auto &cell : dof_handler.active_cell_iterators())
1482  if (cell->is_locally_owned())
1483  {
1484  if (cell->refine_flag_set())
1485  {
1486  // Store the active FE index of each cell that will be
1487  // refined to and distribute it later on its children.
1488  // Pick their future index if flagged for p-refinement.
1489  fe_transfer->refined_cells_fe_index.insert(
1490  {cell, cell->future_fe_index()});
1491  }
1492  else if (cell->coarsen_flag_set())
1493  {
1494  // From all cells that will be coarsened, determine their
1495  // parent and calculate its proper active FE index, so that
1496  // it can be set after refinement. But first, check if that
1497  // particular cell has a parent at all.
1498  Assert(cell->level() > 0, ExcInternalError());
1499  const auto &parent = cell->parent();
1500 
1501  // Check if the active FE index for the current cell has
1502  // been determined already.
1503  if (fe_transfer->coarsened_cells_fe_index.find(parent) ==
1504  fe_transfer->coarsened_cells_fe_index.end())
1505  {
1506  // Find a suitable active FE index for the parent cell
1507  // based on the 'least dominant finite element' of its
1508  // children. Consider the childrens' hypothetical future
1509  // index when they have been flagged for p-refinement.
1510 #ifdef DEBUG
1511  for (const auto &child : parent->child_iterators())
1512  Assert(child->is_active() &&
1513  child->coarsen_flag_set(),
1514  typename ::Triangulation<
1515  dim>::ExcInconsistentCoarseningFlags());
1516 #endif
1517 
1518  const types::fe_index fe_index = ::internal::hp::
1519  DoFHandlerImplementation::Implementation::
1520  dominated_future_fe_on_children<dim, spacedim>(
1521  parent);
1522 
1523  fe_transfer->coarsened_cells_fe_index.insert(
1524  {parent, fe_index});
1525  }
1526  }
1527  else
1528  {
1529  // No h-refinement is scheduled for this cell.
1530  // However, it may have p-refinement indicators, so we
1531  // choose a new active FE index based on its flags.
1532  if (cell->future_fe_index_set() == true)
1533  fe_transfer->persisting_cells_fe_index.insert(
1534  {cell, cell->future_fe_index()});
1535  }
1536  }
1537  }
1538 
1539 
1540 
1545  template <int dim, int spacedim>
1546  static void
1548  DoFHandler<dim, spacedim> &dof_handler)
1549  {
1550  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1551 
1552  // Set active FE indices on persisting cells.
1553  for (const auto &persist : fe_transfer->persisting_cells_fe_index)
1554  {
1555  const auto &cell = persist.first;
1556 
1557  if (cell->is_locally_owned())
1558  {
1559  Assert(cell->is_active(), ExcInternalError());
1560  cell->set_active_fe_index(persist.second);
1561  }
1562  }
1563 
1564  // Distribute active FE indices from all refined cells on their
1565  // respective children.
1566  for (const auto &refine : fe_transfer->refined_cells_fe_index)
1567  {
1568  const auto &parent = refine.first;
1569 
1570  for (const auto &child : parent->child_iterators())
1571  if (child->is_locally_owned())
1572  {
1573  Assert(child->is_active(), ExcInternalError());
1574  child->set_active_fe_index(refine.second);
1575  }
1576  }
1577 
1578  // Set active FE indices on coarsened cells that have been determined
1579  // before the actual coarsening happened.
1580  for (const auto &coarsen : fe_transfer->coarsened_cells_fe_index)
1581  {
1582  const auto &cell = coarsen.first;
1583 
1584  if (cell->is_locally_owned())
1585  {
1586  Assert(cell->is_active(), ExcInternalError());
1587  cell->set_active_fe_index(coarsen.second);
1588  }
1589  }
1590  }
1591 
1592 
1603  template <int dim, int spacedim>
1604  static types::fe_index
1607  const std::vector<types::fe_index> & children_fe_indices,
1608  const ::hp::FECollection<dim, spacedim> &fe_collection)
1609  {
1610  Assert(!children_fe_indices.empty(), ExcInternalError());
1611 
1612  // convert vector to set
1613  // TODO: Change set to types::fe_index
1614  const std::set<unsigned int> children_fe_indices_set(
1615  children_fe_indices.begin(), children_fe_indices.end());
1616 
1617  const types::fe_index dominated_fe_index =
1618  fe_collection.find_dominated_fe_extended(children_fe_indices_set,
1619  /*codim=*/0);
1620 
1621  Assert(dominated_fe_index != numbers::invalid_fe_index,
1623 
1624  return dominated_fe_index;
1625  }
1626 
1627 
1635  template <int dim, int spacedim>
1636  static types::fe_index
1638  const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
1639  {
1640  Assert(
1641  !parent->is_active(),
1642  ExcMessage(
1643  "You ask for information on children of this cell which is only "
1644  "available for active cells. This cell has no children."));
1645 
1646  const auto &dof_handler = parent->get_dof_handler();
1647  Assert(
1648  dof_handler.has_hp_capabilities(),
1650 
1651  // TODO: Change set to types::fe_index
1652  std::set<unsigned int> future_fe_indices_children;
1653  for (const auto &child : parent->child_iterators())
1654  {
1655  Assert(
1656  child->is_active(),
1657  ExcMessage(
1658  "You ask for information on children of this cell which is only "
1659  "available for active cells. One of its children is not active."));
1660 
1661  // Ghost siblings might occur on parallel::shared::Triangulation
1662  // objects. The public interface does not allow to access future
1663  // FE indices on ghost cells. However, we need this information
1664  // here and thus call the internal function that does not check
1665  // for cell ownership. This requires that future FE indices have
1666  // been communicated prior to calling this function.
1667  const types::fe_index future_fe_index_child =
1668  ::internal::DoFCellAccessorImplementation::
1669  Implementation::future_fe_index<dim, spacedim, false>(*child);
1670 
1671  future_fe_indices_children.insert(future_fe_index_child);
1672  }
1673  Assert(!future_fe_indices_children.empty(), ExcInternalError());
1674 
1675  const types::fe_index future_fe_index =
1676  dof_handler.fe_collection.find_dominated_fe_extended(
1677  future_fe_indices_children,
1678  /*codim=*/0);
1679 
1680  Assert(future_fe_index != numbers::invalid_fe_index,
1682 
1683  return future_fe_index;
1684  }
1685  };
1686 
1687 
1688 
1692  template <int dim, int spacedim>
1693  void
1695  {
1696  Implementation::communicate_future_fe_indices<dim, spacedim>(
1697  dof_handler);
1698  }
1699 
1700 
1701 
1705  template <int dim, int spacedim>
1706  unsigned int
1708  const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
1709  {
1710  return Implementation::dominated_future_fe_on_children<dim, spacedim>(
1711  parent);
1712  }
1713  } // namespace DoFHandlerImplementation
1714  } // namespace hp
1715 } // namespace internal
1716 
1717 
1718 
1719 template <int dim, int spacedim>
1720 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
1722  : hp_capability_enabled(true)
1723  , tria(nullptr, typeid(*this).name())
1724  , mg_faces(nullptr)
1725 {}
1726 
1727 
1728 
1729 template <int dim, int spacedim>
1730 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
1732  : DoFHandler()
1733 {
1734  reinit(tria);
1735 }
1736 
1737 
1738 
1739 template <int dim, int spacedim>
1740 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
1742 {
1743  // unsubscribe all attachments to signals of the underlying triangulation
1744  for (auto &connection : this->tria_listeners)
1745  connection.disconnect();
1746  this->tria_listeners.clear();
1747 
1748  for (auto &connection : this->tria_listeners_for_transfer)
1749  connection.disconnect();
1750  this->tria_listeners_for_transfer.clear();
1751 
1752  // release allocated memory
1753  // virtual functions called in constructors and destructors never use the
1754  // override in a derived class
1755  // for clarity be explicit on which function is called
1757 
1758  // also release the policy. this needs to happen before the
1759  // current object disappears because the policy objects
1760  // store references to the DoFhandler object they work on
1761  this->policy.reset();
1762 }
1763 
1764 
1765 
1766 template <int dim, int spacedim>
1767 DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
1769 {
1770  //
1771  // call destructor
1772  //
1773  // remove association with old triangulation
1774  for (auto &connection : this->tria_listeners)
1775  connection.disconnect();
1776  this->tria_listeners.clear();
1777 
1778  for (auto &connection : this->tria_listeners_for_transfer)
1779  connection.disconnect();
1780  this->tria_listeners_for_transfer.clear();
1781 
1782  // release allocated memory and policy
1784  this->policy.reset();
1785 
1786  // reset the finite element collection
1787  this->fe_collection = hp::FECollection<dim, spacedim>();
1788 
1789  //
1790  // call constructor
1791  //
1792  // establish connection to new triangulation
1793  this->tria = &tria;
1794  this->setup_policy();
1795 
1796  // start in hp-mode and let distribute_dofs toggle it if necessary
1797  hp_capability_enabled = true;
1798  this->connect_to_triangulation_signals();
1799  this->create_active_fe_table();
1800 }
1801 
1802 
1803 
1804 /*------------------------ Cell iterator functions ------------------------*/
1805 
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].size() == 0 &&
2820  this->hp_cell_future_fe_indices[level].size() == 0)
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 
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:1599
std::vector< types::fe_index > get_active_fe_indices() const
std::unique_ptr< ActiveFEIndexTransfer > active_fe_index_transfer
Definition: dof_handler.h:1611
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1519
void create_active_fe_table()
std::vector< std::vector< types::fe_index > > hp_cell_future_fe_indices
Definition: dof_handler.h:1586
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:1512
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:1580
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1592
void post_distributed_transfer_action()
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition: dof_handler.h:1550
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:1561
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
Definition: dof_handler.h:1605
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:1506
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:1574
void clear()
std::array< std::vector< types::fe_index >, dim+1 > hp_object_fe_indices
Definition: dof_handler.h:1569
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:1776
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:2480
unsigned int n_vertices() const
unsigned int n_raw_hexs(const unsigned int level) const
unsigned int size() const
Definition: collection.h:264
unsigned int max_dofs_per_line() const
unsigned int max_dofs_per_hex() const
unsigned int max_dofs_per_vertex() const
unsigned int max_dofs_per_quad() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_CXX20_REQUIRES(condition)
Definition: config.h:162
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
unsigned int level
Definition: grid_out.cc:4618
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:1614
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNoFESelected()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
static ::ExceptionBase & ExcNoDominatedFiniteElementOnChildren()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
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:1363
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
@ valid
Iterator points to a valid object.
static const char T
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
void swap(MemorySpaceData< T, MemorySpace > &u, MemorySpaceData< T, MemorySpace > &v)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1351
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:1510
Definition: hp.h:118
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13833
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:54
const types::boundary_id internal_face_boundary_id
Definition: types.h:277
const types::fe_index invalid_fe_index
Definition: types.h:228
static const unsigned int invalid_unsigned_int
Definition: types.h:213
const types::global_dof_index invalid_dof_index
Definition: types.h:233
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:141
static void reserve_space_mg(DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:566
static void reserve_space(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:369
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:310
static unsigned int max_couplings_between_dofs(const DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:109
static void reset_to_empty_objects(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:264
static unsigned int max_couplings_between_dofs(const DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:228
static void reserve_space_mg(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:419
static unsigned int max_couplings_between_dofs(const DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:99
static void reserve_space_mg(DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:489
static void reserve_cells(DoFHandler< dim, spacedim > &dof_handler, const unsigned int n_inner_dofs_per_cell)
Definition: dof_handler.cc:280
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:821
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:668
static void reserve_space_faces(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:866
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:688
static void reserve_space(DoFHandler< 1, spacedim > &dof_handler)
const ::Triangulation< dim, spacedim > & tria