Reference documentation for deal.II version Git e87e13460a 2021-10-28 06:28:57 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_handler.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/config.h>
17 
20 #include <deal.II/base/mpi.templates.h>
21 
22 #include <deal.II/distributed/cell_data_transfer.templates.h>
26 
29 
31 #include <deal.II/grid/tria.h>
35 
36 #include <algorithm>
37 #include <memory>
38 #include <set>
39 #include <unordered_set>
40 
42 
43 template <int dim, int spacedim>
45 
46 
47 template <int dim, int spacedim>
50 
51 
52 namespace internal
53 {
54  template <int dim, int spacedim>
55  std::string
56  policy_to_string(const ::internal::DoFHandlerImplementation::Policy::
57  PolicyBase<dim, spacedim> &policy)
58  {
59  std::string policy_name;
60  if (dynamic_cast<const typename ::internal::DoFHandlerImplementation::
61  Policy::Sequential<dim, spacedim> *>(&policy) ||
62  dynamic_cast<const typename ::internal::DoFHandlerImplementation::
63  Policy::Sequential<dim, spacedim> *>(&policy))
64  policy_name = "Policy::Sequential<";
65  else if (dynamic_cast<
66  const typename ::internal::DoFHandlerImplementation::
67  Policy::ParallelDistributed<dim, spacedim> *>(&policy) ||
68  dynamic_cast<
69  const typename ::internal::DoFHandlerImplementation::
70  Policy::ParallelDistributed<dim, spacedim> *>(&policy))
71  policy_name = "Policy::ParallelDistributed<";
72  else if (dynamic_cast<
73  const typename ::internal::DoFHandlerImplementation::
74  Policy::ParallelShared<dim, spacedim> *>(&policy) ||
75  dynamic_cast<
76  const typename ::internal::DoFHandlerImplementation::
77  Policy::ParallelShared<dim, spacedim> *>(&policy))
78  policy_name = "Policy::ParallelShared<";
79  else
81  policy_name += Utilities::int_to_string(dim) + "," +
82  Utilities::int_to_string(spacedim) + ">";
83  return policy_name;
84  }
85 
86 
87  namespace DoFHandlerImplementation
88  {
94  {
99  template <int spacedim>
100  static unsigned int
102  {
103  return std::min(static_cast<types::global_dof_index>(
104  3 * dof_handler.fe_collection.max_dofs_per_vertex() +
105  2 * dof_handler.fe_collection.max_dofs_per_line()),
106  dof_handler.n_dofs());
107  }
108 
109  template <int spacedim>
110  static unsigned int
112  {
113  // get these numbers by drawing pictures
114  // and counting...
115  // example:
116  // | | |
117  // --x-----x--x--X--
118  // | | | |
119  // | x--x--x
120  // | | | |
121  // --x--x--*--x--x--
122  // | | | |
123  // x--x--x |
124  // | | | |
125  // --X--x--x-----x--
126  // | | |
127  // x = vertices connected with center vertex *;
128  // = total of 19
129  // (the X vertices are connected with * if
130  // the vertices adjacent to X are hanging
131  // nodes)
132  // count lines -> 28 (don't forget to count
133  // mother and children separately!)
134  types::global_dof_index max_couplings;
135  switch (dof_handler.tria->max_adjacent_cells())
136  {
137  case 4:
138  max_couplings =
139  19 * dof_handler.fe_collection.max_dofs_per_vertex() +
140  28 * dof_handler.fe_collection.max_dofs_per_line() +
141  8 * dof_handler.fe_collection.max_dofs_per_quad();
142  break;
143  case 5:
144  max_couplings =
145  21 * dof_handler.fe_collection.max_dofs_per_vertex() +
146  31 * dof_handler.fe_collection.max_dofs_per_line() +
147  9 * dof_handler.fe_collection.max_dofs_per_quad();
148  break;
149  case 6:
150  max_couplings =
151  28 * dof_handler.fe_collection.max_dofs_per_vertex() +
152  42 * dof_handler.fe_collection.max_dofs_per_line() +
153  12 * dof_handler.fe_collection.max_dofs_per_quad();
154  break;
155  case 7:
156  max_couplings =
157  30 * dof_handler.fe_collection.max_dofs_per_vertex() +
158  45 * dof_handler.fe_collection.max_dofs_per_line() +
159  13 * dof_handler.fe_collection.max_dofs_per_quad();
160  break;
161  case 8:
162  max_couplings =
163  37 * dof_handler.fe_collection.max_dofs_per_vertex() +
164  56 * dof_handler.fe_collection.max_dofs_per_line() +
165  16 * dof_handler.fe_collection.max_dofs_per_quad();
166  break;
167 
168  // the following numbers are not based on actual counting but by
169  // extrapolating the number sequences from the previous ones (for
170  // example, for n_dofs_per_vertex(), the sequence above is 19, 21,
171  // 28, 30, 37, and is continued as follows):
172  case 9:
173  max_couplings =
174  39 * dof_handler.fe_collection.max_dofs_per_vertex() +
175  59 * dof_handler.fe_collection.max_dofs_per_line() +
176  17 * dof_handler.fe_collection.max_dofs_per_quad();
177  break;
178  case 10:
179  max_couplings =
180  46 * dof_handler.fe_collection.max_dofs_per_vertex() +
181  70 * dof_handler.fe_collection.max_dofs_per_line() +
182  20 * dof_handler.fe_collection.max_dofs_per_quad();
183  break;
184  case 11:
185  max_couplings =
186  48 * dof_handler.fe_collection.max_dofs_per_vertex() +
187  73 * dof_handler.fe_collection.max_dofs_per_line() +
188  21 * dof_handler.fe_collection.max_dofs_per_quad();
189  break;
190  case 12:
191  max_couplings =
192  55 * dof_handler.fe_collection.max_dofs_per_vertex() +
193  84 * dof_handler.fe_collection.max_dofs_per_line() +
194  24 * dof_handler.fe_collection.max_dofs_per_quad();
195  break;
196  case 13:
197  max_couplings =
198  57 * dof_handler.fe_collection.max_dofs_per_vertex() +
199  87 * dof_handler.fe_collection.max_dofs_per_line() +
200  25 * dof_handler.fe_collection.max_dofs_per_quad();
201  break;
202  case 14:
203  max_couplings =
204  63 * dof_handler.fe_collection.max_dofs_per_vertex() +
205  98 * dof_handler.fe_collection.max_dofs_per_line() +
206  28 * dof_handler.fe_collection.max_dofs_per_quad();
207  break;
208  case 15:
209  max_couplings =
210  65 * dof_handler.fe_collection.max_dofs_per_vertex() +
211  103 * dof_handler.fe_collection.max_dofs_per_line() +
212  29 * dof_handler.fe_collection.max_dofs_per_quad();
213  break;
214  case 16:
215  max_couplings =
216  72 * dof_handler.fe_collection.max_dofs_per_vertex() +
217  114 * dof_handler.fe_collection.max_dofs_per_line() +
218  32 * dof_handler.fe_collection.max_dofs_per_quad();
219  break;
220 
221  default:
222  Assert(false, ExcNotImplemented());
223  max_couplings = 0;
224  }
225  return std::min(max_couplings, dof_handler.n_dofs());
226  }
227 
228  template <int spacedim>
229  static unsigned int
231  {
232  // TODO:[?] Invent significantly better estimates than the ones in this
233  // function
234 
235  // doing the same thing here is a rather complicated thing, compared
236  // to the 2d case, since it is hard to draw pictures with several
237  // refined hexahedra :-) so I presently only give a coarse
238  // estimate for the case that at most 8 hexes meet at each vertex
239  //
240  // can anyone give better estimate here?
241  const unsigned int max_adjacent_cells =
242  dof_handler.tria->max_adjacent_cells();
243 
244  types::global_dof_index max_couplings;
245  if (max_adjacent_cells <= 8)
246  max_couplings =
247  7 * 7 * 7 * dof_handler.fe_collection.max_dofs_per_vertex() +
248  7 * 6 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_line() +
249  9 * 4 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_quad() +
250  27 * dof_handler.fe_collection.max_dofs_per_hex();
251  else
252  {
253  Assert(false, ExcNotImplemented());
254  max_couplings = 0;
255  }
256 
257  return std::min(max_couplings, dof_handler.n_dofs());
258  }
259 
264  template <int dim, int spacedim>
265  static void
267  {
268  dof_handler.object_dof_indices.clear();
269  dof_handler.object_dof_indices.resize(dof_handler.tria->n_levels());
270  dof_handler.object_dof_indices.shrink_to_fit();
271 
272  dof_handler.object_dof_ptr.clear();
273  dof_handler.object_dof_ptr.resize(dof_handler.tria->n_levels());
274  dof_handler.object_dof_ptr.shrink_to_fit();
275 
276  dof_handler.cell_dof_cache_indices.clear();
277  dof_handler.cell_dof_cache_indices.resize(dof_handler.tria->n_levels());
278  dof_handler.cell_dof_cache_indices.shrink_to_fit();
279 
280  dof_handler.cell_dof_cache_ptr.clear();
281  dof_handler.cell_dof_cache_ptr.resize(dof_handler.tria->n_levels());
282  dof_handler.cell_dof_cache_ptr.shrink_to_fit();
283  }
284 
288  template <int dim, int spacedim>
289  static void
291  const unsigned int n_inner_dofs_per_cell)
292  {
293  for (unsigned int i = 0; i < dof_handler.tria->n_levels(); ++i)
294  {
295  // 1) object_dof_indices
296  dof_handler.object_dof_ptr[i][dim].assign(
297  dof_handler.tria->n_raw_cells(i) + 1, 0);
298 
299  for (const auto &cell :
300  dof_handler.tria->cell_iterators_on_level(i))
301  if (cell->is_active() && !cell->is_artificial())
302  dof_handler.object_dof_ptr[i][dim][cell->index() + 1] =
303  n_inner_dofs_per_cell;
304 
305  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i); ++j)
306  dof_handler.object_dof_ptr[i][dim][j + 1] +=
307  dof_handler.object_dof_ptr[i][dim][j];
308 
309  dof_handler.object_dof_indices[i][dim].resize(
310  dof_handler.object_dof_ptr[i][dim].back(),
312 
313  // 2) cell_dof_cache_indices
314  dof_handler.cell_dof_cache_ptr[i].assign(
315  dof_handler.tria->n_raw_cells(i) + 1, 0);
316 
317  for (const auto &cell :
318  dof_handler.tria->cell_iterators_on_level(i))
319  if (cell->is_active() && !cell->is_artificial())
320  dof_handler.cell_dof_cache_ptr[i][cell->index() + 1] =
321  dof_handler.get_fe().n_dofs_per_cell();
322 
323  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i); ++j)
324  dof_handler.cell_dof_cache_ptr[i][j + 1] +=
325  dof_handler.cell_dof_cache_ptr[i][j];
326 
327  dof_handler.cell_dof_cache_indices[i].resize(
328  dof_handler.cell_dof_cache_ptr[i].back(),
330  }
331  }
332 
337  template <int dim, int spacedim, typename T>
338  static void
340  const unsigned int structdim,
341  const unsigned int n_raw_entities,
342  const T & cell_process)
343  {
344  if (dof_handler.tria->n_cells() == 0)
345  return;
346 
347  dof_handler.object_dof_ptr[0][structdim].assign(n_raw_entities + 1, -1);
348  // determine for each entity the number of dofs
349  for (const auto &cell : dof_handler.tria->cell_iterators())
350  if (cell->is_active() && !cell->is_artificial())
351  cell_process(
352  cell,
353  [&](const unsigned int n_dofs_per_entity,
354  const unsigned int index) {
355  auto &n_dofs_per_entity_target =
356  dof_handler.object_dof_ptr[0][structdim][index + 1];
357 
358  // make sure that either the entity has not been visited or
359  // the entity has the same number of dofs assigned
360  Assert((n_dofs_per_entity_target ==
361  static_cast<
363  -1) ||
364  n_dofs_per_entity_target == n_dofs_per_entity),
366 
367  n_dofs_per_entity_target = n_dofs_per_entity;
368  });
369 
370  // convert the absolute numbers to CRS
371  dof_handler.object_dof_ptr[0][structdim][0] = 0;
372  for (unsigned int i = 1; i < n_raw_entities + 1; ++i)
373  {
374  if (dof_handler.object_dof_ptr[0][structdim][i] ==
375  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
376  -1))
377  dof_handler.object_dof_ptr[0][structdim][i] =
378  dof_handler.object_dof_ptr[0][structdim][i - 1];
379  else
380  dof_handler.object_dof_ptr[0][structdim][i] +=
381  dof_handler.object_dof_ptr[0][structdim][i - 1];
382  }
383 
384  // allocate memory for indices
385  dof_handler.object_dof_indices[0][structdim].resize(
386  dof_handler.object_dof_ptr[0][structdim].back(),
388  }
389 
396  template <int dim, int spacedim>
397  static void
399  {
400  reset_to_empty_objects(dof_handler);
401 
402  const auto &fe = dof_handler.get_fe();
403 
404  // cell
405  reserve_cells(dof_handler,
406  dim == 1 ? fe.n_dofs_per_line() :
407  (dim == 2 ? fe.n_dofs_per_quad(0) :
408  fe.n_dofs_per_hex()));
409 
410  // vertices
411  reserve_subentities(dof_handler,
412  0,
413  dof_handler.tria->n_vertices(),
414  [&](const auto &cell, const auto &process) {
415  for (const auto vertex_index :
416  cell->vertex_indices())
417  process(fe.n_dofs_per_vertex(),
418  cell->vertex_index(vertex_index));
419  });
420 
421  // lines
422  if (dim == 2 || dim == 3)
423  reserve_subentities(dof_handler,
424  1,
425  dof_handler.tria->n_raw_lines(),
426  [&](const auto &cell, const auto &process) {
427  for (const auto line_index :
428  cell->line_indices())
429  process(fe.n_dofs_per_line(),
430  cell->line(line_index)->index());
431  });
432 
433  // quads
434  if (dim == 3)
435  reserve_subentities(dof_handler,
436  2,
437  dof_handler.tria->n_raw_quads(),
438  [&](const auto &cell, const auto &process) {
439  for (const auto face_index :
440  cell->face_indices())
441  process(fe.n_dofs_per_quad(face_index),
442  cell->face(face_index)->index());
443  });
444  }
445 
446  template <int spacedim>
447  static void
449  {
450  Assert(dof_handler.get_triangulation().n_levels() > 0,
451  ExcMessage("Invalid triangulation"));
452  dof_handler.clear_mg_space();
453 
454  const ::Triangulation<1, spacedim> &tria =
455  dof_handler.get_triangulation();
456  const unsigned int dofs_per_line =
457  dof_handler.get_fe().n_dofs_per_line();
458  const unsigned int n_levels = tria.n_levels();
459 
460  for (unsigned int i = 0; i < n_levels; ++i)
461  {
462  dof_handler.mg_levels.emplace_back(
464  dof_handler.mg_levels.back()->dof_object.dofs =
465  std::vector<types::global_dof_index>(tria.n_raw_lines(i) *
466  dofs_per_line,
468  }
469 
470  const unsigned int n_vertices = tria.n_vertices();
471 
472  dof_handler.mg_vertex_dofs.resize(n_vertices);
473 
474  std::vector<unsigned int> max_level(n_vertices, 0);
475  std::vector<unsigned int> min_level(n_vertices, n_levels);
476 
477  for (typename ::Triangulation<1, spacedim>::cell_iterator cell =
478  tria.begin();
479  cell != tria.end();
480  ++cell)
481  {
482  const unsigned int level = cell->level();
483 
484  for (const auto vertex : cell->vertex_indices())
485  {
486  const unsigned int vertex_index = cell->vertex_index(vertex);
487 
488  if (min_level[vertex_index] > level)
489  min_level[vertex_index] = level;
490 
491  if (max_level[vertex_index] < level)
492  max_level[vertex_index] = level;
493  }
494  }
495 
496  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
497  if (tria.vertex_used(vertex))
498  {
499  Assert(min_level[vertex] < n_levels, ExcInternalError());
500  Assert(max_level[vertex] >= min_level[vertex],
501  ExcInternalError());
502  dof_handler.mg_vertex_dofs[vertex].init(
503  min_level[vertex],
504  max_level[vertex],
505  dof_handler.get_fe().n_dofs_per_vertex());
506  }
507 
508  else
509  {
510  Assert(min_level[vertex] == n_levels, ExcInternalError());
511  Assert(max_level[vertex] == 0, ExcInternalError());
512  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
513  }
514  }
515 
516  template <int spacedim>
517  static void
519  {
520  Assert(dof_handler.get_triangulation().n_levels() > 0,
521  ExcMessage("Invalid triangulation"));
522  dof_handler.clear_mg_space();
523 
524  const ::FiniteElement<2, spacedim> &fe = dof_handler.get_fe();
525  const ::Triangulation<2, spacedim> &tria =
526  dof_handler.get_triangulation();
527  const unsigned int n_levels = tria.n_levels();
528 
529  for (unsigned int i = 0; i < n_levels; ++i)
530  {
531  dof_handler.mg_levels.emplace_back(
532  std::make_unique<
534  dof_handler.mg_levels.back()->dof_object.dofs =
535  std::vector<types::global_dof_index>(
536  tria.n_raw_quads(i) *
537  fe.n_dofs_per_quad(0 /*note: in 2D there is only one quad*/),
539  }
540 
541  dof_handler.mg_faces =
542  std::make_unique<internal::DoFHandlerImplementation::DoFFaces<2>>();
543  dof_handler.mg_faces->lines.dofs =
544  std::vector<types::global_dof_index>(tria.n_raw_lines() *
545  fe.n_dofs_per_line(),
547 
548  const unsigned int n_vertices = tria.n_vertices();
549 
550  dof_handler.mg_vertex_dofs.resize(n_vertices);
551 
552  std::vector<unsigned int> max_level(n_vertices, 0);
553  std::vector<unsigned int> min_level(n_vertices, n_levels);
554 
555  for (typename ::Triangulation<2, spacedim>::cell_iterator cell =
556  tria.begin();
557  cell != tria.end();
558  ++cell)
559  {
560  const unsigned int level = cell->level();
561 
562  for (const auto vertex : cell->vertex_indices())
563  {
564  const unsigned int vertex_index = cell->vertex_index(vertex);
565 
566  if (min_level[vertex_index] > level)
567  min_level[vertex_index] = level;
568 
569  if (max_level[vertex_index] < level)
570  max_level[vertex_index] = level;
571  }
572  }
573 
574  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
575  if (tria.vertex_used(vertex))
576  {
577  Assert(min_level[vertex] < n_levels, ExcInternalError());
578  Assert(max_level[vertex] >= min_level[vertex],
579  ExcInternalError());
580  dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
581  max_level[vertex],
582  fe.n_dofs_per_vertex());
583  }
584 
585  else
586  {
587  Assert(min_level[vertex] == n_levels, ExcInternalError());
588  Assert(max_level[vertex] == 0, ExcInternalError());
589  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
590  }
591  }
592 
593  template <int spacedim>
594  static void
596  {
597  Assert(dof_handler.get_triangulation().n_levels() > 0,
598  ExcMessage("Invalid triangulation"));
599  dof_handler.clear_mg_space();
600 
601  const ::FiniteElement<3, spacedim> &fe = dof_handler.get_fe();
602  const ::Triangulation<3, spacedim> &tria =
603  dof_handler.get_triangulation();
604  const unsigned int n_levels = tria.n_levels();
605 
606  for (unsigned int i = 0; i < n_levels; ++i)
607  {
608  dof_handler.mg_levels.emplace_back(
609  std::make_unique<
611  dof_handler.mg_levels.back()->dof_object.dofs =
612  std::vector<types::global_dof_index>(tria.n_raw_hexs(i) *
613  fe.n_dofs_per_hex(),
615  }
616 
617  dof_handler.mg_faces =
618  std::make_unique<internal::DoFHandlerImplementation::DoFFaces<3>>();
619  dof_handler.mg_faces->lines.dofs =
620  std::vector<types::global_dof_index>(tria.n_raw_lines() *
621  fe.n_dofs_per_line(),
623 
624  // TODO: the implementation makes the assumption that all faces have the
625  // same number of dofs
626  AssertDimension(fe.n_unique_faces(), 1);
627  dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index>(
628  tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/),
630 
631  const unsigned int n_vertices = tria.n_vertices();
632 
633  dof_handler.mg_vertex_dofs.resize(n_vertices);
634 
635  std::vector<unsigned int> max_level(n_vertices, 0);
636  std::vector<unsigned int> min_level(n_vertices, n_levels);
637 
638  for (typename ::Triangulation<3, spacedim>::cell_iterator cell =
639  tria.begin();
640  cell != tria.end();
641  ++cell)
642  {
643  const unsigned int level = cell->level();
644 
645  for (const auto vertex : cell->vertex_indices())
646  {
647  const unsigned int vertex_index = cell->vertex_index(vertex);
648 
649  if (min_level[vertex_index] > level)
650  min_level[vertex_index] = level;
651 
652  if (max_level[vertex_index] < level)
653  max_level[vertex_index] = level;
654  }
655  }
656 
657  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
658  if (tria.vertex_used(vertex))
659  {
660  Assert(min_level[vertex] < n_levels, ExcInternalError());
661  Assert(max_level[vertex] >= min_level[vertex],
662  ExcInternalError());
663  dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
664  max_level[vertex],
665  fe.n_dofs_per_vertex());
666  }
667 
668  else
669  {
670  Assert(min_level[vertex] == n_levels, ExcInternalError());
671  Assert(max_level[vertex] == 0, ExcInternalError());
672  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
673  }
674  }
675 
676  template <int spacedim>
679  const DoFHandler<1, spacedim> &dof_handler,
681  &mg_level,
683  &,
684  const unsigned int obj_index,
685  const unsigned int fe_index,
686  const unsigned int local_index,
687  const std::integral_constant<int, 1>)
688  {
689  Assert(dof_handler.hp_capability_enabled == false,
691 
692  return mg_level->dof_object.get_dof_index(
693  static_cast<const DoFHandler<1, spacedim> &>(dof_handler),
694  obj_index,
695  fe_index,
696  local_index);
697  }
698 
699  template <int spacedim>
702  const DoFHandler<2, spacedim> &dof_handler,
704  &,
706  & mg_faces,
707  const unsigned int obj_index,
708  const unsigned int fe_index,
709  const unsigned int local_index,
710  const std::integral_constant<int, 1>)
711  {
712  return mg_faces->lines.get_dof_index(
713  static_cast<const DoFHandler<2, spacedim> &>(dof_handler),
714  obj_index,
715  fe_index,
716  local_index);
717  }
718 
719  template <int spacedim>
722  const DoFHandler<2, spacedim> &dof_handler,
724  &mg_level,
726  &,
727  const unsigned int obj_index,
728  const unsigned int fe_index,
729  const unsigned int local_index,
730  const std::integral_constant<int, 2>)
731  {
732  Assert(dof_handler.hp_capability_enabled == false,
734  return mg_level->dof_object.get_dof_index(
735  static_cast<const DoFHandler<2, spacedim> &>(dof_handler),
736  obj_index,
737  fe_index,
738  local_index);
739  }
740 
741  template <int spacedim>
744  const DoFHandler<3, spacedim> &dof_handler,
746  &,
748  & mg_faces,
749  const unsigned int obj_index,
750  const unsigned int fe_index,
751  const unsigned int local_index,
752  const std::integral_constant<int, 1>)
753  {
754  Assert(dof_handler.hp_capability_enabled == false,
756  return mg_faces->lines.get_dof_index(
757  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
758  obj_index,
759  fe_index,
760  local_index);
761  }
762 
763  template <int spacedim>
766  const DoFHandler<3, spacedim> &dof_handler,
768  &,
770  & mg_faces,
771  const unsigned int obj_index,
772  const unsigned int fe_index,
773  const unsigned int local_index,
774  const std::integral_constant<int, 2>)
775  {
776  Assert(dof_handler.hp_capability_enabled == false,
778  return mg_faces->quads.get_dof_index(
779  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
780  obj_index,
781  fe_index,
782  local_index);
783  }
784 
785  template <int spacedim>
788  const DoFHandler<3, spacedim> &dof_handler,
790  &mg_level,
792  &,
793  const unsigned int obj_index,
794  const unsigned int fe_index,
795  const unsigned int local_index,
796  const std::integral_constant<int, 3>)
797  {
798  Assert(dof_handler.hp_capability_enabled == false,
800  return mg_level->dof_object.get_dof_index(
801  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
802  obj_index,
803  fe_index,
804  local_index);
805  }
806 
807  template <int spacedim>
808  static void
810  const DoFHandler<1, spacedim> &dof_handler,
812  &mg_level,
814  &,
815  const unsigned int obj_index,
816  const unsigned int fe_index,
817  const unsigned int local_index,
819  const std::integral_constant<int, 1>)
820  {
821  Assert(dof_handler.hp_capability_enabled == false,
823  mg_level->dof_object.set_dof_index(
824  static_cast<const DoFHandler<1, spacedim> &>(dof_handler),
825  obj_index,
826  fe_index,
827  local_index,
828  global_index);
829  }
830 
831  template <int spacedim>
832  static void
834  const DoFHandler<2, spacedim> &dof_handler,
836  &,
838  & mg_faces,
839  const unsigned int obj_index,
840  const unsigned int fe_index,
841  const unsigned int local_index,
843  const std::integral_constant<int, 1>)
844  {
845  Assert(dof_handler.hp_capability_enabled == false,
847  mg_faces->lines.set_dof_index(
848  static_cast<const DoFHandler<2, spacedim> &>(dof_handler),
849  obj_index,
850  fe_index,
851  local_index,
852  global_index);
853  }
854 
855  template <int spacedim>
856  static void
858  const DoFHandler<2, spacedim> &dof_handler,
860  &mg_level,
862  &,
863  const unsigned int obj_index,
864  const unsigned int fe_index,
865  const unsigned int local_index,
867  const std::integral_constant<int, 2>)
868  {
869  Assert(dof_handler.hp_capability_enabled == false,
871  mg_level->dof_object.set_dof_index(
872  static_cast<const DoFHandler<2, spacedim> &>(dof_handler),
873  obj_index,
874  fe_index,
875  local_index,
876  global_index);
877  }
878 
879  template <int spacedim>
880  static void
882  const DoFHandler<3, spacedim> &dof_handler,
884  &,
886  & mg_faces,
887  const unsigned int obj_index,
888  const unsigned int fe_index,
889  const unsigned int local_index,
891  const std::integral_constant<int, 1>)
892  {
893  Assert(dof_handler.hp_capability_enabled == false,
895  mg_faces->lines.set_dof_index(
896  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
897  obj_index,
898  fe_index,
899  local_index,
900  global_index);
901  }
902 
903  template <int spacedim>
904  static void
906  const DoFHandler<3, spacedim> &dof_handler,
908  &,
910  & mg_faces,
911  const unsigned int obj_index,
912  const unsigned int fe_index,
913  const unsigned int local_index,
915  const std::integral_constant<int, 2>)
916  {
917  Assert(dof_handler.hp_capability_enabled == false,
919  mg_faces->quads.set_dof_index(
920  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
921  obj_index,
922  fe_index,
923  local_index,
924  global_index);
925  }
926 
927  template <int spacedim>
928  static void
930  const DoFHandler<3, spacedim> &dof_handler,
932  &mg_level,
934  &,
935  const unsigned int obj_index,
936  const unsigned int fe_index,
937  const unsigned int local_index,
939  const std::integral_constant<int, 3>)
940  {
941  Assert(dof_handler.hp_capability_enabled == false,
943  mg_level->dof_object.set_dof_index(
944  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
945  obj_index,
946  fe_index,
947  local_index,
948  global_index);
949  }
950  };
951  } // namespace DoFHandlerImplementation
952 
953 
954 
955  namespace hp
956  {
957  namespace DoFHandlerImplementation
958  {
964  {
970  template <int dim, int spacedim>
971  static void
973  DoFHandler<dim, spacedim> &dof_handler)
974  {
975  (void)dof_handler;
976  for (const auto &cell : dof_handler.active_cell_iterators())
977  if (cell->is_locally_owned())
978  Assert(
979  !cell->future_fe_index_set(),
980  ExcMessage(
981  "There shouldn't be any cells flagged for p-adaptation when partitioning."));
982  }
983 
984 
985 
990  template <int dim, int spacedim>
991  static void
993  {
994  // The final step in all of the reserve_space() functions is to set
995  // up vertex dof information. since vertices are sequentially
996  // numbered, what we do first is to set up an array in which
997  // we record whether a vertex is associated with any of the
998  // given fe's, by setting a bit. in a later step, we then
999  // actually allocate memory for the required dofs
1000  //
1001  // in the following, we only need to consider vertices that are
1002  // adjacent to either a locally owned or a ghost cell; we never
1003  // store anything on vertices that are only surrounded by
1004  // artificial cells. so figure out that subset of vertices
1005  // first
1006  std::vector<bool> locally_used_vertices(
1007  dof_handler.tria->n_vertices(), false);
1008  for (const auto &cell : dof_handler.active_cell_iterators())
1009  if (!cell->is_artificial())
1010  for (const auto v : cell->vertex_indices())
1011  locally_used_vertices[cell->vertex_index(v)] = true;
1012 
1013  std::vector<std::vector<bool>> vertex_fe_association(
1014  dof_handler.fe_collection.size(),
1015  std::vector<bool>(dof_handler.tria->n_vertices(), false));
1016 
1017  for (const auto &cell : dof_handler.active_cell_iterators())
1018  if (!cell->is_artificial())
1019  for (const auto v : cell->vertex_indices())
1020  vertex_fe_association[cell->active_fe_index()]
1021  [cell->vertex_index(v)] = true;
1022 
1023  // in debug mode, make sure that each vertex is associated
1024  // with at least one FE (note that except for unused
1025  // vertices, all vertices are actually active). this is of
1026  // course only true for vertices that are part of either
1027  // ghost or locally owned cells
1028 #ifdef DEBUG
1029  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
1030  if (locally_used_vertices[v] == true)
1031  if (dof_handler.tria->vertex_used(v) == true)
1032  {
1033  unsigned int fe = 0;
1034  for (; fe < dof_handler.fe_collection.size(); ++fe)
1035  if (vertex_fe_association[fe][v] == true)
1036  break;
1037  Assert(fe != dof_handler.fe_collection.size(),
1038  ExcInternalError());
1039  }
1040 #endif
1041 
1042  const unsigned int d = 0;
1043  const unsigned int l = 0;
1044 
1045  dof_handler.hp_object_fe_ptr[d].clear();
1046  dof_handler.hp_object_fe_indices[d].clear();
1047  dof_handler.object_dof_ptr[l][d].clear();
1048  dof_handler.object_dof_indices[l][d].clear();
1049 
1050  dof_handler.hp_object_fe_ptr[d].reserve(
1051  dof_handler.tria->n_vertices() + 1);
1052 
1053  unsigned int vertex_slots_needed = 0;
1054  unsigned int fe_slots_needed = 0;
1055 
1056  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
1057  {
1058  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1059 
1060  if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
1061  {
1062  for (unsigned int fe = 0;
1063  fe < dof_handler.fe_collection.size();
1064  ++fe)
1065  if (vertex_fe_association[fe][v] == true)
1066  {
1067  fe_slots_needed++;
1068  vertex_slots_needed +=
1069  dof_handler.get_fe(fe).n_dofs_per_vertex();
1070  }
1071  }
1072  }
1073 
1074  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1075 
1076  dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
1077  dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
1078 
1079  dof_handler.object_dof_indices[l][d].reserve(vertex_slots_needed);
1080 
1081  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
1082  if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
1083  {
1084  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
1085  ++fe)
1086  if (vertex_fe_association[fe][v] == true)
1087  {
1088  dof_handler.hp_object_fe_indices[d].push_back(fe);
1089  dof_handler.object_dof_ptr[l][d].push_back(
1090  dof_handler.object_dof_indices[l][d].size());
1091 
1092  for (unsigned int i = 0;
1093  i < dof_handler.get_fe(fe).n_dofs_per_vertex();
1094  i++)
1095  dof_handler.object_dof_indices[l][d].push_back(
1097  }
1098  }
1099 
1100 
1101  dof_handler.object_dof_ptr[l][d].push_back(
1102  dof_handler.object_dof_indices[l][d].size());
1103 
1104  AssertDimension(vertex_slots_needed,
1105  dof_handler.object_dof_indices[l][d].size());
1106  AssertDimension(fe_slots_needed,
1107  dof_handler.hp_object_fe_indices[d].size());
1108  AssertDimension(fe_slots_needed + 1,
1109  dof_handler.object_dof_ptr[l][d].size());
1110  AssertDimension(dof_handler.tria->n_vertices() + 1,
1111  dof_handler.hp_object_fe_ptr[d].size());
1112 
1113  dof_handler.object_dof_indices[l][d].assign(
1114  vertex_slots_needed, numbers::invalid_dof_index);
1115  }
1116 
1117 
1118 
1123  template <int dim, int spacedim>
1124  static void
1126  {
1127  (void)dof_handler;
1128  // count how much space we need on each level for the cell
1129  // dofs and set the dof_*_offsets data. initially set the
1130  // latter to an invalid index, and only later set it to
1131  // something reasonable for active dof_handler.cells
1132  //
1133  // note that for dof_handler.cells, the situation is simpler
1134  // than for other (lower dimensional) objects since exactly
1135  // one finite element is used for it
1136  for (unsigned int level = 0; level < dof_handler.tria->n_levels();
1137  ++level)
1138  {
1139  dof_handler.object_dof_ptr[level][dim] =
1140  std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
1141  dof_handler.tria->n_raw_cells(level),
1142  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
1143  -1));
1144  dof_handler.cell_dof_cache_ptr[level] =
1145  std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
1146  dof_handler.tria->n_raw_cells(level),
1147  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
1148  -1));
1149 
1150  types::global_dof_index next_free_dof = 0;
1151  types::global_dof_index cache_size = 0;
1152 
1153  for (auto cell :
1155  if (cell->is_active() && !cell->is_artificial())
1156  {
1157  dof_handler.object_dof_ptr[level][dim][cell->index()] =
1158  next_free_dof;
1159  next_free_dof +=
1160  cell->get_fe().template n_dofs_per_object<dim>();
1161 
1162  dof_handler.cell_dof_cache_ptr[level][cell->index()] =
1163  cache_size;
1164  cache_size += cell->get_fe().n_dofs_per_cell();
1165  }
1166 
1167  dof_handler.object_dof_indices[level][dim] =
1168  std::vector<types::global_dof_index>(
1169  next_free_dof, numbers::invalid_dof_index);
1170  dof_handler.cell_dof_cache_indices[level] =
1171  std::vector<types::global_dof_index>(
1172  cache_size, numbers::invalid_dof_index);
1173  }
1174  }
1175 
1176 
1177 
1182  template <int dim, int spacedim>
1183  static void
1185  {
1186  // FACE DOFS
1187  //
1188  // Count face dofs, then allocate as much space
1189  // as we need and prime the linked list for faces (see the
1190  // description in hp::DoFLevel) with the indices we will
1191  // need. Note that our task is more complicated than for the
1192  // cell case above since two adjacent cells may have different
1193  // active FE indices, in which case we need to allocate
1194  // *two* sets of face dofs for the same face. But they don't
1195  // *have* to be different, and so we need to prepare for this
1196  // as well.
1197  //
1198  // The way we do things is that we loop over all active
1199  // cells (these are the only ones that have DoFs
1200  // anyway) and all their faces. We note in the
1201  // user flags whether we have previously visited a face and
1202  // if so skip it (consequently, we have to save and later
1203  // restore the face flags)
1204  {
1205  std::vector<bool> saved_face_user_flags;
1206  switch (dim)
1207  {
1208  case 2:
1209  {
1210  const_cast<::Triangulation<dim, spacedim> &>(
1211  *dof_handler.tria)
1212  .save_user_flags_line(saved_face_user_flags);
1213  const_cast<::Triangulation<dim, spacedim> &>(
1214  *dof_handler.tria)
1215  .clear_user_flags_line();
1216 
1217  break;
1218  }
1219 
1220  case 3:
1221  {
1222  const_cast<::Triangulation<dim, spacedim> &>(
1223  *dof_handler.tria)
1224  .save_user_flags_quad(saved_face_user_flags);
1225  const_cast<::Triangulation<dim, spacedim> &>(
1226  *dof_handler.tria)
1227  .clear_user_flags_quad();
1228 
1229  break;
1230  }
1231 
1232  default:
1233  Assert(false, ExcNotImplemented());
1234  }
1235 
1236  const unsigned int d = dim - 1;
1237  const unsigned int l = 0;
1238 
1239  dof_handler.hp_object_fe_ptr[d].clear();
1240  dof_handler.hp_object_fe_indices[d].clear();
1241  dof_handler.object_dof_ptr[l][d].clear();
1242  dof_handler.object_dof_indices[l][d].clear();
1243 
1244  dof_handler.hp_object_fe_ptr[d].resize(
1245  dof_handler.tria->n_raw_faces() + 1);
1246 
1247  // An array to hold how many slots (see the hp::DoFLevel
1248  // class) we will have to store on each level
1249  unsigned int n_face_slots = 0;
1250 
1251  for (const auto &cell : dof_handler.active_cell_iterators())
1252  if (!cell->is_artificial())
1253  for (const auto face : cell->face_indices())
1254  if (cell->face(face)->user_flag_set() == false)
1255  {
1256  unsigned int fe_slots_needed = 0;
1257 
1258  if (cell->at_boundary(face) ||
1259  cell->face(face)->has_children() ||
1260  cell->neighbor_is_coarser(face) ||
1261  (!cell->at_boundary(face) &&
1262  cell->neighbor(face)->is_artificial()) ||
1263  (!cell->at_boundary(face) &&
1264  !cell->neighbor(face)->is_artificial() &&
1265  (cell->active_fe_index() ==
1266  cell->neighbor(face)->active_fe_index())))
1267  {
1268  fe_slots_needed = 1;
1269  n_face_slots +=
1270  dof_handler.get_fe(cell->active_fe_index())
1271  .template n_dofs_per_object<dim - 1>(face);
1272  }
1273  else
1274  {
1275  fe_slots_needed = 2;
1276  n_face_slots +=
1277  dof_handler.get_fe(cell->active_fe_index())
1278  .template n_dofs_per_object<dim - 1>(face) +
1279  dof_handler
1280  .get_fe(cell->neighbor(face)->active_fe_index())
1281  .template n_dofs_per_object<dim - 1>(
1282  cell->neighbor_face_no(face));
1283  }
1284 
1285  // mark this face as visited
1286  cell->face(face)->set_user_flag();
1287 
1288  dof_handler
1289  .hp_object_fe_ptr[d][cell->face(face)->index() + 1] =
1290  fe_slots_needed;
1291  }
1292 
1293  for (unsigned int i = 1; i < dof_handler.hp_object_fe_ptr[d].size();
1294  i++)
1295  dof_handler.hp_object_fe_ptr[d][i] +=
1296  dof_handler.hp_object_fe_ptr[d][i - 1];
1297 
1298 
1299  dof_handler.hp_object_fe_indices[d].resize(
1300  dof_handler.hp_object_fe_ptr[d].back());
1301  dof_handler.object_dof_ptr[l][d].resize(
1302  dof_handler.hp_object_fe_ptr[d].back() + 1);
1303 
1304  dof_handler.object_dof_indices[l][d].reserve(n_face_slots);
1305 
1306 
1307  // With the memory now allocated, loop over the
1308  // dof_handler cells again and prime the _offset values as
1309  // well as the fe_index fields
1310  switch (dim)
1311  {
1312  case 2:
1313  {
1314  const_cast<::Triangulation<dim, spacedim> &>(
1315  *dof_handler.tria)
1316  .clear_user_flags_line();
1317 
1318  break;
1319  }
1320 
1321  case 3:
1322  {
1323  const_cast<::Triangulation<dim, spacedim> &>(
1324  *dof_handler.tria)
1325  .clear_user_flags_quad();
1326 
1327  break;
1328  }
1329 
1330  default:
1331  Assert(false, ExcNotImplemented());
1332  }
1333 
1334  for (const auto &cell : dof_handler.active_cell_iterators())
1335  if (!cell->is_artificial())
1336  for (const auto face : cell->face_indices())
1337  if (!cell->face(face)->user_flag_set())
1338  {
1339  // Same decision tree as before
1340  if (cell->at_boundary(face) ||
1341  cell->face(face)->has_children() ||
1342  cell->neighbor_is_coarser(face) ||
1343  (!cell->at_boundary(face) &&
1344  cell->neighbor(face)->is_artificial()) ||
1345  (!cell->at_boundary(face) &&
1346  !cell->neighbor(face)->is_artificial() &&
1347  (cell->active_fe_index() ==
1348  cell->neighbor(face)->active_fe_index())))
1349  {
1350  const unsigned int fe = cell->active_fe_index();
1351  const unsigned int n_dofs =
1352  dof_handler.get_fe(fe)
1353  .template n_dofs_per_object<dim - 1>(face);
1354  const unsigned int offset =
1355  dof_handler
1356  .hp_object_fe_ptr[d][cell->face(face)->index()];
1357 
1358  dof_handler.hp_object_fe_indices[d][offset] = fe;
1359  dof_handler.object_dof_ptr[l][d][offset + 1] = n_dofs;
1360 
1361  for (unsigned int i = 0; i < n_dofs; ++i)
1362  dof_handler.object_dof_indices[l][d].push_back(
1364  }
1365  else
1366  {
1367  unsigned int fe_1 = cell->active_fe_index();
1368  unsigned int face_no_1 = face;
1369  unsigned int fe_2 =
1370  cell->neighbor(face)->active_fe_index();
1371  unsigned int face_no_2 = cell->neighbor_face_no(face);
1372 
1373  if (fe_2 < fe_1)
1374  {
1375  std::swap(fe_1, fe_2);
1376  std::swap(face_no_1, face_no_2);
1377  }
1378 
1379  const unsigned int n_dofs_1 =
1380  dof_handler.get_fe(fe_1)
1381  .template n_dofs_per_object<dim - 1>(face_no_1);
1382 
1383  const unsigned int n_dofs_2 =
1384  dof_handler.get_fe(fe_2)
1385  .template n_dofs_per_object<dim - 1>(face_no_2);
1386 
1387  const unsigned int offset =
1388  dof_handler
1389  .hp_object_fe_ptr[d][cell->face(face)->index()];
1390 
1391  dof_handler.hp_object_fe_indices[d].push_back(
1392  cell->active_fe_index());
1393  dof_handler.object_dof_ptr[l][d].push_back(
1394  dof_handler.object_dof_indices[l][d].size());
1395 
1396  dof_handler.hp_object_fe_indices[d][offset + 0] =
1397  fe_1;
1398  dof_handler.hp_object_fe_indices[d][offset + 1] =
1399  fe_2;
1400  dof_handler.object_dof_ptr[l][d][offset + 1] =
1401  n_dofs_1;
1402  dof_handler.object_dof_ptr[l][d][offset + 2] =
1403  n_dofs_2;
1404 
1405 
1406  for (unsigned int i = 0; i < n_dofs_1 + n_dofs_2; ++i)
1407  dof_handler.object_dof_indices[l][d].push_back(
1409  }
1410 
1411  // mark this face as visited
1412  cell->face(face)->set_user_flag();
1413  }
1414 
1415  for (unsigned int i = 1;
1416  i < dof_handler.object_dof_ptr[l][d].size();
1417  i++)
1418  dof_handler.object_dof_ptr[l][d][i] +=
1419  dof_handler.object_dof_ptr[l][d][i - 1];
1420 
1421  // at the end, restore the user flags for the faces
1422  switch (dim)
1423  {
1424  case 2:
1425  {
1426  const_cast<::Triangulation<dim, spacedim> &>(
1427  *dof_handler.tria)
1428  .load_user_flags_line(saved_face_user_flags);
1429 
1430  break;
1431  }
1432 
1433  case 3:
1434  {
1435  const_cast<::Triangulation<dim, spacedim> &>(
1436  *dof_handler.tria)
1437  .load_user_flags_quad(saved_face_user_flags);
1438 
1439  break;
1440  }
1441 
1442  default:
1443  Assert(false, ExcNotImplemented());
1444  }
1445  }
1446  }
1447 
1448 
1449 
1456  template <int spacedim>
1457  static void
1459  {
1460  Assert(dof_handler.fe_collection.size() > 0,
1462  Assert(dof_handler.tria->n_levels() > 0,
1463  ExcMessage("The current Triangulation must not be empty."));
1464  Assert(dof_handler.tria->n_levels() ==
1465  dof_handler.hp_cell_future_fe_indices.size(),
1466  ExcInternalError());
1467 
1469  reset_to_empty_objects(dof_handler);
1470 
1471  Threads::TaskGroup<> tasks;
1472  tasks +=
1473  Threads::new_task(&reserve_space_cells<1, spacedim>, dof_handler);
1474  tasks += Threads::new_task(&reserve_space_vertices<1, spacedim>,
1475  dof_handler);
1476  tasks.join_all();
1477  }
1478 
1479 
1480 
1481  template <int spacedim>
1482  static void
1484  {
1485  Assert(dof_handler.fe_collection.size() > 0,
1487  Assert(dof_handler.tria->n_levels() > 0,
1488  ExcMessage("The current Triangulation must not be empty."));
1489  Assert(dof_handler.tria->n_levels() ==
1490  dof_handler.hp_cell_future_fe_indices.size(),
1491  ExcInternalError());
1492 
1494  reset_to_empty_objects(dof_handler);
1495 
1496  Threads::TaskGroup<> tasks;
1497  tasks +=
1498  Threads::new_task(&reserve_space_cells<2, spacedim>, dof_handler);
1499  tasks +=
1500  Threads::new_task(&reserve_space_faces<2, spacedim>, dof_handler);
1501  tasks += Threads::new_task(&reserve_space_vertices<2, spacedim>,
1502  dof_handler);
1503  tasks.join_all();
1504  }
1505 
1506 
1507 
1508  template <int spacedim>
1509  static void
1511  {
1512  Assert(dof_handler.fe_collection.size() > 0,
1514  Assert(dof_handler.tria->n_levels() > 0,
1515  ExcMessage("The current Triangulation must not be empty."));
1516  Assert(dof_handler.tria->n_levels() ==
1517  dof_handler.hp_cell_future_fe_indices.size(),
1518  ExcInternalError());
1519 
1521  reset_to_empty_objects(dof_handler);
1522 
1523  Threads::TaskGroup<> tasks;
1524  tasks +=
1525  Threads::new_task(&reserve_space_cells<3, spacedim>, dof_handler);
1526  tasks +=
1527  Threads::new_task(&reserve_space_faces<3, spacedim>, dof_handler);
1528  tasks += Threads::new_task(&reserve_space_vertices<3, spacedim>,
1529  dof_handler);
1530 
1531  // While the tasks above are running, we can turn to line dofs
1532 
1533  // the situation here is pretty much like with vertices:
1534  // there can be an arbitrary number of finite elements
1535  // associated with each line.
1536  //
1537  // the algorithm we use is somewhat similar to what we do in
1538  // reserve_space_vertices()
1539  {
1540  // what we do first is to set up an array in which we
1541  // record whether a line is associated with any of the
1542  // given fe's, by setting a bit. in a later step, we
1543  // then actually allocate memory for the required dofs
1544  std::vector<std::vector<bool>> line_fe_association(
1545  dof_handler.fe_collection.size(),
1546  std::vector<bool>(dof_handler.tria->n_raw_lines(), false));
1547 
1548  for (const auto &cell : dof_handler.active_cell_iterators())
1549  if (!cell->is_artificial())
1550  for (const auto l : cell->line_indices())
1551  line_fe_association[cell->active_fe_index()]
1552  [cell->line_index(l)] = true;
1553 
1554  // first check which of the lines is used at all,
1555  // i.e. is associated with a finite element. we do this
1556  // since not all lines may actually be used, in which
1557  // case we do not have to allocate any memory at all
1558  std::vector<bool> line_is_used(dof_handler.tria->n_raw_lines(),
1559  false);
1560  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1561  ++line)
1562  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
1563  ++fe)
1564  if (line_fe_association[fe][line] == true)
1565  {
1566  line_is_used[line] = true;
1567  break;
1568  }
1569 
1570 
1571 
1572  const unsigned int d = 1;
1573  const unsigned int l = 0;
1574 
1575  dof_handler.hp_object_fe_ptr[d].clear();
1576  dof_handler.hp_object_fe_indices[d].clear();
1577  dof_handler.object_dof_ptr[l][d].clear();
1578  dof_handler.object_dof_indices[l][d].clear();
1579 
1580  dof_handler.hp_object_fe_ptr[d].reserve(
1581  dof_handler.tria->n_raw_lines() + 1);
1582 
1583  unsigned int line_slots_needed = 0;
1584  unsigned int fe_slots_needed = 0;
1585 
1586  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1587  ++line)
1588  {
1589  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1590 
1591  if (line_is_used[line] == true)
1592  {
1593  for (unsigned int fe = 0;
1594  fe < dof_handler.fe_collection.size();
1595  ++fe)
1596  if (line_fe_association[fe][line] == true)
1597  {
1598  fe_slots_needed++;
1599  line_slots_needed +=
1600  dof_handler.get_fe(fe).n_dofs_per_line();
1601  }
1602  }
1603  }
1604 
1605  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1606 
1607  // make sure that all entries have been set
1608  AssertDimension(dof_handler.hp_object_fe_ptr[d].size(),
1609  dof_handler.tria->n_raw_lines() + 1);
1610 
1611  dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
1612  dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
1613 
1614  dof_handler.object_dof_indices[l][d].reserve(line_slots_needed);
1615 
1616  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1617  ++line)
1618  if (line_is_used[line] == true)
1619  {
1620  for (unsigned int fe = 0;
1621  fe < dof_handler.fe_collection.size();
1622  ++fe)
1623  if (line_fe_association[fe][line] == true)
1624  {
1625  dof_handler.hp_object_fe_indices[d].push_back(fe);
1626  dof_handler.object_dof_ptr[l][d].push_back(
1627  dof_handler.object_dof_indices[l][d].size());
1628 
1629  for (unsigned int i = 0;
1630  i < dof_handler.get_fe(fe).n_dofs_per_line();
1631  i++)
1632  dof_handler.object_dof_indices[l][d].push_back(
1634  }
1635  }
1636 
1637  dof_handler.object_dof_ptr[l][d].push_back(
1638  dof_handler.object_dof_indices[l][d].size());
1639 
1640  // make sure that all entries have been set
1641  AssertDimension(dof_handler.hp_object_fe_indices[d].size(),
1642  fe_slots_needed);
1643  AssertDimension(dof_handler.object_dof_ptr[l][d].size(),
1644  fe_slots_needed + 1);
1645  AssertDimension(dof_handler.object_dof_indices[l][d].size(),
1646  line_slots_needed);
1647  }
1648 
1649  // Ensure that everything is done at this point.
1650  tasks.join_all();
1651  }
1652 
1653 
1654 
1666  template <int dim, int spacedim>
1667  static void
1669  {
1670  Assert(
1671  dof_handler.hp_capability_enabled == true,
1673 
1674  using active_fe_index_type =
1675  typename ::DoFHandler<dim, spacedim>::active_fe_index_type;
1676 
1677  if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1678  dynamic_cast<
1679  const ::parallel::shared::Triangulation<dim, spacedim>
1680  *>(&dof_handler.get_triangulation()))
1681  {
1682  // we have a shared triangulation. in this case, every processor
1683  // knows about all cells, but every processor only has knowledge
1684  // about the active FE index on the cells it owns.
1685  //
1686  // we can create a complete set of active FE indices by letting
1687  // every processor create a vector of indices for all cells,
1688  // filling only those on the cells it owns and setting the indices
1689  // on the other cells to zero. then we add all of these vectors
1690  // up, and because every vector entry has exactly one processor
1691  // that owns it, the sum is correct
1692  std::vector<active_fe_index_type> active_fe_indices(
1693  tr->n_active_cells(), 0u);
1694  for (const auto &cell : dof_handler.active_cell_iterators())
1695  if (cell->is_locally_owned())
1696  active_fe_indices[cell->active_cell_index()] =
1697  cell->active_fe_index();
1698 
1699  Utilities::MPI::sum(active_fe_indices,
1700  tr->get_communicator(),
1701  active_fe_indices);
1702 
1703  // now go back and fill the active FE index on all other
1704  // cells. we would like to call cell->set_active_fe_index(),
1705  // but that function does not allow setting these indices on
1706  // non-locally_owned cells. so we have to work around the
1707  // issue a little bit by accessing the underlying data
1708  // structures directly
1709  for (const auto &cell : dof_handler.active_cell_iterators())
1710  if (!cell->is_locally_owned())
1711  dof_handler
1712  .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1713  active_fe_indices[cell->active_cell_index()];
1714  }
1715  else if (const ::parallel::
1716  DistributedTriangulationBase<dim, spacedim> *tr =
1717  dynamic_cast<
1718  const ::parallel::
1719  DistributedTriangulationBase<dim, spacedim> *>(
1720  &dof_handler.get_triangulation()))
1721  {
1722  // For completely distributed meshes, use the function that is
1723  // able to move data from locally owned cells on one processor to
1724  // the corresponding ghost cells on others. To this end, we need
1725  // to have functions that can pack and unpack the data we want to
1726  // transport -- namely, the single unsigned int active_fe_index
1727  // objects
1728  auto pack =
1729  [](const typename ::DoFHandler<dim, spacedim>::
1730  active_cell_iterator &cell) -> active_fe_index_type {
1731  return cell->active_fe_index();
1732  };
1733 
1734  auto unpack =
1735  [&dof_handler](
1736  const typename ::DoFHandler<dim, spacedim>::
1737  active_cell_iterator & cell,
1738  const active_fe_index_type active_fe_index) -> void {
1739  // we would like to say
1740  // cell->set_active_fe_index(active_fe_index);
1741  // but this is not allowed on cells that are not
1742  // locally owned, and we are on a ghost cell
1743  dof_handler
1744  .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1745  active_fe_index;
1746  };
1747 
1749  active_fe_index_type,
1750  ::DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1751  }
1752  else
1753  {
1754  // a sequential triangulation. there is nothing we need to do here
1755  Assert(
1756  (dynamic_cast<
1757  const ::parallel::TriangulationBase<dim, spacedim> *>(
1758  &dof_handler.get_triangulation()) == nullptr),
1759  ExcInternalError());
1760  }
1761  }
1762 
1763 
1764 
1778  template <int dim, int spacedim>
1779  static void
1781  {
1782  Assert(
1783  dof_handler.hp_capability_enabled == true,
1785 
1786  using active_fe_index_type =
1787  typename ::DoFHandler<dim, spacedim>::active_fe_index_type;
1788 
1789  if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1790  dynamic_cast<
1791  const ::parallel::shared::Triangulation<dim, spacedim>
1792  *>(&dof_handler.get_triangulation()))
1793  {
1794  std::vector<active_fe_index_type> future_fe_indices(
1795  tr->n_active_cells(), 0u);
1796  for (const auto &cell : dof_handler.active_cell_iterators())
1797  if (cell->is_locally_owned())
1798  future_fe_indices[cell->active_cell_index()] =
1799  dof_handler
1800  .hp_cell_future_fe_indices[cell->level()][cell->index()];
1801 
1802  Utilities::MPI::sum(future_fe_indices,
1803  tr->get_communicator(),
1804  future_fe_indices);
1805 
1806  for (const auto &cell : dof_handler.active_cell_iterators())
1807  if (!cell->is_locally_owned())
1808  dof_handler
1809  .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1810  future_fe_indices[cell->active_cell_index()];
1811  }
1812  else if (const ::parallel::
1813  DistributedTriangulationBase<dim, spacedim> *tr =
1814  dynamic_cast<
1815  const ::parallel::
1816  DistributedTriangulationBase<dim, spacedim> *>(
1817  &dof_handler.get_triangulation()))
1818  {
1819  auto pack =
1820  [&dof_handler](
1821  const typename ::DoFHandler<dim, spacedim>::
1822  active_cell_iterator &cell) -> active_fe_index_type {
1823  return dof_handler
1824  .hp_cell_future_fe_indices[cell->level()][cell->index()];
1825  };
1826 
1827  auto unpack =
1828  [&dof_handler](
1829  const typename ::DoFHandler<dim, spacedim>::
1830  active_cell_iterator & cell,
1831  const active_fe_index_type future_fe_index) -> void {
1832  dof_handler
1833  .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1834  future_fe_index;
1835  };
1836 
1838  active_fe_index_type,
1839  ::DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1840  }
1841  else
1842  {
1843  Assert(
1844  (dynamic_cast<
1845  const ::parallel::TriangulationBase<dim, spacedim> *>(
1846  &dof_handler.get_triangulation()) == nullptr),
1847  ExcInternalError());
1848  }
1849  }
1850 
1851 
1852 
1873  template <int dim, int spacedim>
1874  static void
1876  DoFHandler<dim, spacedim> &dof_handler)
1877  {
1878  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1879 
1880  for (const auto &cell : dof_handler.active_cell_iterators())
1881  if (cell->is_locally_owned())
1882  {
1883  if (cell->refine_flag_set())
1884  {
1885  // Store the active FE index of each cell that will be
1886  // refined to and distribute it later on its children.
1887  // Pick their future index if flagged for p-refinement.
1888  fe_transfer->refined_cells_fe_index.insert(
1889  {cell, cell->future_fe_index()});
1890  }
1891  else if (cell->coarsen_flag_set())
1892  {
1893  // From all cells that will be coarsened, determine their
1894  // parent and calculate its proper active FE index, so that
1895  // it can be set after refinement. But first, check if that
1896  // particular cell has a parent at all.
1897  Assert(cell->level() > 0, ExcInternalError());
1898  const auto &parent = cell->parent();
1899 
1900  // Check if the active FE index for the current cell has
1901  // been determined already.
1902  if (fe_transfer->coarsened_cells_fe_index.find(parent) ==
1903  fe_transfer->coarsened_cells_fe_index.end())
1904  {
1905  // Find a suitable active FE index for the parent cell
1906  // based on the 'least dominant finite element' of its
1907  // children. Consider the childrens' hypothetical future
1908  // index when they have been flagged for p-refinement.
1909 #ifdef DEBUG
1910  for (const auto &child : parent->child_iterators())
1911  Assert(child->is_active() &&
1912  child->coarsen_flag_set(),
1913  typename ::Triangulation<
1914  dim>::ExcInconsistentCoarseningFlags());
1915 #endif
1916 
1917  const unsigned int fe_index = ::internal::hp::
1918  DoFHandlerImplementation::Implementation::
1919  dominated_future_fe_on_children<dim, spacedim>(
1920  parent);
1921 
1922  fe_transfer->coarsened_cells_fe_index.insert(
1923  {parent, fe_index});
1924  }
1925  }
1926  else
1927  {
1928  // No h-refinement is scheduled for this cell.
1929  // However, it may have p-refinement indicators, so we
1930  // choose a new active FE index based on its flags.
1931  if (cell->future_fe_index_set() == true)
1932  fe_transfer->persisting_cells_fe_index.insert(
1933  {cell, cell->future_fe_index()});
1934  }
1935  }
1936  }
1937 
1938 
1939 
1944  template <int dim, int spacedim>
1945  static void
1947  DoFHandler<dim, spacedim> &dof_handler)
1948  {
1949  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1950 
1951  // Set active FE indices on persisting cells.
1952  for (const auto &persist : fe_transfer->persisting_cells_fe_index)
1953  {
1954  const auto &cell = persist.first;
1955 
1956  if (cell->is_locally_owned())
1957  {
1958  Assert(cell->is_active(), ExcInternalError());
1959  cell->set_active_fe_index(persist.second);
1960  }
1961  }
1962 
1963  // Distribute active FE indices from all refined cells on their
1964  // respective children.
1965  for (const auto &refine : fe_transfer->refined_cells_fe_index)
1966  {
1967  const auto &parent = refine.first;
1968 
1969  for (const auto &child : parent->child_iterators())
1970  if (child->is_locally_owned())
1971  {
1972  Assert(child->is_active(), ExcInternalError());
1973  child->set_active_fe_index(refine.second);
1974  }
1975  }
1976 
1977  // Set active FE indices on coarsened cells that have been determined
1978  // before the actual coarsening happened.
1979  for (const auto &coarsen : fe_transfer->coarsened_cells_fe_index)
1980  {
1981  const auto &cell = coarsen.first;
1982 
1983  if (cell->is_locally_owned())
1984  {
1985  Assert(cell->is_active(), ExcInternalError());
1986  cell->set_active_fe_index(coarsen.second);
1987  }
1988  }
1989  }
1990 
1991 
2002  template <int dim, int spacedim>
2003  static unsigned int
2006  const std::vector<unsigned int> & children_fe_indices,
2007  const ::hp::FECollection<dim, spacedim> &fe_collection)
2008  {
2009  Assert(!children_fe_indices.empty(), ExcInternalError());
2010 
2011  // convert vector to set
2012  const std::set<unsigned int> children_fe_indices_set(
2013  children_fe_indices.begin(), children_fe_indices.end());
2014 
2015  const unsigned int dominated_fe_index =
2016  fe_collection.find_dominated_fe_extended(children_fe_indices_set,
2017  /*codim=*/0);
2018 
2019  Assert(dominated_fe_index != numbers::invalid_unsigned_int,
2021 
2022  return dominated_fe_index;
2023  }
2024 
2025 
2033  template <int dim, int spacedim>
2034  static unsigned int
2036  const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
2037  {
2038  Assert(
2039  !parent->is_active(),
2040  ExcMessage(
2041  "You ask for information on children of this cell which is only "
2042  "available for active cells. This cell has no children."));
2043 
2044  const auto &dof_handler = parent->get_dof_handler();
2045  Assert(
2046  dof_handler.has_hp_capabilities(),
2048 
2049  std::set<unsigned int> future_fe_indices_children;
2050  for (const auto &child : parent->child_iterators())
2051  {
2052  Assert(
2053  child->is_active(),
2054  ExcMessage(
2055  "You ask for information on children of this cell which is only "
2056  "available for active cells. One of its children is not active."));
2057 
2058  // Ghost siblings might occur on parallel::shared::Triangulation
2059  // objects. The public interface does not allow to access future
2060  // FE indices on ghost cells. However, we need this information
2061  // here and thus call the internal function that does not check
2062  // for cell ownership. This requires that future FE indices have
2063  // been communicated prior to calling this function.
2064  const unsigned int future_fe_index_child =
2065  ::internal::DoFCellAccessorImplementation::
2066  Implementation::future_fe_index<dim, spacedim, false>(*child);
2067 
2068  future_fe_indices_children.insert(future_fe_index_child);
2069  }
2070  Assert(!future_fe_indices_children.empty(), ExcInternalError());
2071 
2072  const unsigned int future_fe_index =
2073  dof_handler.fe_collection.find_dominated_fe_extended(
2074  future_fe_indices_children,
2075  /*codim=*/0);
2076 
2077  Assert(future_fe_index != numbers::invalid_unsigned_int,
2079 
2080  return future_fe_index;
2081  }
2082  };
2083 
2084 
2085 
2089  template <int dim, int spacedim>
2090  void
2092  {
2093  Implementation::communicate_future_fe_indices<dim, spacedim>(
2094  dof_handler);
2095  }
2096 
2097 
2098 
2102  template <int dim, int spacedim>
2103  unsigned int
2105  const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
2106  {
2107  return Implementation::dominated_future_fe_on_children<dim, spacedim>(
2108  parent);
2109  }
2110  } // namespace DoFHandlerImplementation
2111  } // namespace hp
2112 } // namespace internal
2113 
2114 
2115 
2116 template <int dim, int spacedim>
2118  : hp_capability_enabled(true)
2119  , tria(nullptr, typeid(*this).name())
2120  , mg_faces(nullptr)
2121 {}
2122 
2123 
2124 
2125 template <int dim, int spacedim>
2127  : DoFHandler()
2128 {
2129  reinit(tria);
2130 }
2131 
2132 
2133 
2134 template <int dim, int spacedim>
2136 {
2137  // unsubscribe all attachments to signals of the underlying triangulation
2138  for (auto &connection : this->tria_listeners)
2139  connection.disconnect();
2140  this->tria_listeners.clear();
2141 
2142  for (auto &connection : this->tria_listeners_for_transfer)
2143  connection.disconnect();
2144  this->tria_listeners_for_transfer.clear();
2145 
2146  // release allocated memory
2147  // virtual functions called in constructors and destructors never use the
2148  // override in a derived class
2149  // for clarity be explicit on which function is called
2151 
2152  // also release the policy. this needs to happen before the
2153  // current object disappears because the policy objects
2154  // store references to the DoFhandler object they work on
2155  this->policy.reset();
2156 }
2157 
2158 
2159 
2160 template <int dim, int spacedim>
2161 void
2163  const FiniteElement<dim, spacedim> &fe)
2164 {
2166 }
2167 
2168 
2169 
2170 template <int dim, int spacedim>
2171 void
2174 {
2175  this->reinit(tria);
2176  this->distribute_dofs(fe);
2177 }
2178 
2179 
2180 
2181 template <int dim, int spacedim>
2182 void
2184 {
2185  //
2186  // call destructor
2187  //
2188  // remove association with old triangulation
2189  for (auto &connection : this->tria_listeners)
2190  connection.disconnect();
2191  this->tria_listeners.clear();
2192 
2193  for (auto &connection : this->tria_listeners_for_transfer)
2194  connection.disconnect();
2195  this->tria_listeners_for_transfer.clear();
2196 
2197  // release allocated memory and policy
2199  this->policy.reset();
2200 
2201  // reset the finite element collection
2203 
2204  //
2205  // call constructor
2206  //
2207  // establish connection to new triangulation
2208  this->tria = &tria;
2209  this->setup_policy();
2210 
2211  // start in hp-mode and let distribute_dofs toggle it if necessary
2212  hp_capability_enabled = true;
2214  this->create_active_fe_table();
2215 }
2216 
2217 
2218 
2219 /*------------------------ Cell iterator functions ------------------------*/
2220 
2221 template <int dim, int spacedim>
2223 DoFHandler<dim, spacedim>::begin(const unsigned int level) const
2224 {
2226  this->get_triangulation().begin(level);
2227  if (cell == this->get_triangulation().end(level))
2228  return end(level);
2229  return cell_iterator(*cell, this);
2230 }
2231 
2232 
2233 
2234 template <int dim, int spacedim>
2237 {
2238  // level is checked in begin
2239  cell_iterator i = begin(level);
2240  if (i.state() != IteratorState::valid)
2241  return i;
2242  while (i->has_children())
2243  if ((++i).state() != IteratorState::valid)
2244  return i;
2245  return i;
2246 }
2247 
2248 
2249 
2250 template <int dim, int spacedim>
2253 {
2254  return cell_iterator(&this->get_triangulation(), -1, -1, this);
2255 }
2256 
2257 
2258 
2259 template <int dim, int spacedim>
2261 DoFHandler<dim, spacedim>::end(const unsigned int level) const
2262 {
2264  this->get_triangulation().end(level);
2265  if (cell.state() != IteratorState::valid)
2266  return end();
2267  return cell_iterator(*cell, this);
2268 }
2269 
2270 
2271 
2272 template <int dim, int spacedim>
2275 {
2277  this->get_triangulation().end_active(level);
2278  if (cell.state() != IteratorState::valid)
2279  return active_cell_iterator(end());
2280  return active_cell_iterator(*cell, this);
2281 }
2282 
2283 
2284 
2285 template <int dim, int spacedim>
2288 {
2289  Assert(this->has_level_dofs(),
2290  ExcMessage("You can only iterate over mg "
2291  "levels if mg dofs got distributed."));
2293  this->get_triangulation().begin(level);
2294  if (cell == this->get_triangulation().end(level))
2295  return end_mg(level);
2296  return level_cell_iterator(*cell, this);
2297 }
2298 
2299 
2300 
2301 template <int dim, int spacedim>
2303 DoFHandler<dim, spacedim>::end_mg(const unsigned int level) const
2304 {
2305  Assert(this->has_level_dofs(),
2306  ExcMessage("You can only iterate over mg "
2307  "levels if mg dofs got distributed."));
2309  this->get_triangulation().end(level);
2310  if (cell.state() != IteratorState::valid)
2311  return end();
2312  return level_cell_iterator(*cell, this);
2313 }
2314 
2315 
2316 
2317 template <int dim, int spacedim>
2320 {
2321  return level_cell_iterator(&this->get_triangulation(), -1, -1, this);
2322 }
2323 
2324 
2325 
2326 template <int dim, int spacedim>
2329 {
2331  begin(), end());
2332 }
2333 
2334 
2335 
2336 template <int dim, int spacedim>
2339 {
2340  return IteratorRange<
2342  end());
2343 }
2344 
2345 
2346 
2347 template <int dim, int spacedim>
2350 {
2352  begin_mg(), end_mg());
2353 }
2354 
2355 
2356 
2357 template <int dim, int spacedim>
2360  const unsigned int level) const
2361 {
2363  begin(level), end(level));
2364 }
2365 
2366 
2367 
2368 template <int dim, int spacedim>
2371  const unsigned int level) const
2372 {
2373  return IteratorRange<
2375  begin_active(level), end_active(level));
2376 }
2377 
2378 
2379 
2380 template <int dim, int spacedim>
2383  const unsigned int level) const
2384 {
2386  begin_mg(level), end_mg(level));
2387 }
2388 
2389 
2390 
2391 //---------------------------------------------------------------------------
2392 
2393 
2394 
2395 template <int dim, int spacedim>
2398 {
2399  Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2401 
2402  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2403 
2404  std::unordered_set<types::global_dof_index> boundary_dofs;
2405  std::vector<types::global_dof_index> dofs_on_face;
2406  dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2407 
2408  const IndexSet &owned_dofs = locally_owned_dofs();
2409 
2410  // loop over all faces to check whether they are at a
2411  // boundary. note that we need not take special care of single
2412  // lines in 3d (using @p{cell->has_boundary_lines}), since we do
2413  // not support boundaries of dimension dim-2, and so every
2414  // boundary line is also part of a boundary face.
2415  for (const auto &cell : this->active_cell_iterators())
2416  if (cell->is_locally_owned() && cell->at_boundary())
2417  {
2418  for (const auto iface : cell->face_indices())
2419  {
2420  const auto face = cell->face(iface);
2421  if (face->at_boundary())
2422  {
2423  const unsigned int dofs_per_face =
2424  cell->get_fe().n_dofs_per_face(iface);
2425  dofs_on_face.resize(dofs_per_face);
2426 
2427  face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2428  for (unsigned int i = 0; i < dofs_per_face; ++i)
2429  {
2430  const unsigned int global_idof_index = dofs_on_face[i];
2431  if (owned_dofs.is_element(global_idof_index))
2432  {
2433  boundary_dofs.insert(global_idof_index);
2434  }
2435  }
2436  }
2437  }
2438  }
2439  return boundary_dofs.size();
2440 }
2441 
2442 
2443 
2444 template <int dim, int spacedim>
2447  const std::set<types::boundary_id> &boundary_ids) const
2448 {
2449  Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2451 
2452  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2453  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2454  boundary_ids.end(),
2456 
2457  // same as above, but with additional checks for set of boundary
2458  // indicators
2459  std::unordered_set<types::global_dof_index> boundary_dofs;
2460  std::vector<types::global_dof_index> dofs_on_face;
2461  dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2462 
2463  const IndexSet &owned_dofs = locally_owned_dofs();
2464 
2465  for (const auto &cell : this->active_cell_iterators())
2466  if (cell->is_locally_owned() && cell->at_boundary())
2467  {
2468  for (const auto iface : cell->face_indices())
2469  {
2470  const auto face = cell->face(iface);
2471  const unsigned int boundary_id = face->boundary_id();
2472  if (face->at_boundary() &&
2473  (boundary_ids.find(boundary_id) != boundary_ids.end()))
2474  {
2475  const unsigned int dofs_per_face =
2476  cell->get_fe().n_dofs_per_face(iface);
2477  dofs_on_face.resize(dofs_per_face);
2478 
2479  face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2480  for (unsigned int i = 0; i < dofs_per_face; ++i)
2481  {
2482  const unsigned int global_idof_index = dofs_on_face[i];
2483  if (owned_dofs.is_element(global_idof_index))
2484  {
2485  boundary_dofs.insert(global_idof_index);
2486  }
2487  }
2488  }
2489  }
2490  }
2491  return boundary_dofs.size();
2492 }
2493 
2494 
2495 
2496 template <int dim, int spacedim>
2497 std::size_t
2499 {
2500  std::size_t mem = MemoryConsumption::memory_consumption(this->tria) +
2503 
2512 
2513 
2515  {
2516  // nothing to add
2517  }
2518  else
2519  {
2520  // collect size of multigrid data structures
2521 
2523 
2524  for (unsigned int level = 0; level < this->mg_levels.size(); ++level)
2525  mem += this->mg_levels[level]->memory_consumption();
2526 
2527  if (this->mg_faces != nullptr)
2529 
2530  for (unsigned int i = 0; i < this->mg_vertex_dofs.size(); ++i)
2531  mem += sizeof(MGVertexDoFs) +
2532  (1 + this->mg_vertex_dofs[i].get_finest_level() -
2533  this->mg_vertex_dofs[i].get_coarsest_level()) *
2534  sizeof(types::global_dof_index);
2535  }
2536 
2537  return mem;
2538 }
2539 
2540 
2541 
2542 template <int dim, int spacedim>
2543 void
2545 {
2547 }
2548 
2549 
2550 
2551 template <int dim, int spacedim>
2552 void
2554 {
2555  Assert(
2556  this->tria != nullptr,
2557  ExcMessage(
2558  "You need to set the Triangulation in the DoFHandler using reinit() or "
2559  "in the constructor before you can distribute DoFs."));
2560  Assert(this->tria->n_levels() > 0,
2561  ExcMessage("The Triangulation you are using is empty!"));
2562  Assert(ff.size() > 0, ExcMessage("The hp::FECollection given is empty!"));
2563 
2564  // don't create a new object if the one we have is already appropriate
2565  if (this->fe_collection != ff)
2566  {
2568 
2569  const bool contains_multiple_fes = (this->fe_collection.size() > 1);
2570 
2571  // disable hp-mode if only a single finite element has been registered
2572  if (hp_capability_enabled && !contains_multiple_fes)
2573  {
2574  hp_capability_enabled = false;
2575 
2576  // unsubscribe connections to signals that are only relevant for
2577  // hp-mode, since we only have a single element here
2578  for (auto &connection : this->tria_listeners_for_transfer)
2579  connection.disconnect();
2580  this->tria_listeners_for_transfer.clear();
2581 
2582  // release active and future finite element tables
2583  this->hp_cell_active_fe_indices.clear();
2584  this->hp_cell_active_fe_indices.shrink_to_fit();
2585  this->hp_cell_future_fe_indices.clear();
2586  this->hp_cell_future_fe_indices.shrink_to_fit();
2587  }
2588 
2589  // re-enabling hp-mode is not permitted since the active and future FE
2590  // tables are no longer available
2591  AssertThrow(
2592  hp_capability_enabled || !contains_multiple_fes,
2593  ExcMessage(
2594  "You cannot re-enable hp-capabilities after you registered a single "
2595  "finite element. Please create a new DoFHandler object instead."));
2596  }
2597 
2599  {
2600  // make sure every processor knows the active FE indices
2601  // on both its own cells and all ghost cells
2604 
2605  // make sure that the FE collection is large enough to
2606  // cover all FE indices presently in use on the mesh
2607  for (const auto &cell : this->active_cell_iterators())
2608  if (!cell->is_artificial())
2609  Assert(cell->active_fe_index() < this->fe_collection.size(),
2610  ExcInvalidFEIndex(cell->active_fe_index(),
2611  this->fe_collection.size()));
2612  }
2613 }
2614 
2615 
2616 
2617 template <int dim, int spacedim>
2618 void
2620  const FiniteElement<dim, spacedim> &fe)
2621 {
2623 }
2624 
2625 
2626 
2627 template <int dim, int spacedim>
2628 void
2631 {
2632  Assert(
2633  this->tria != nullptr,
2634  ExcMessage(
2635  "You need to set the Triangulation in the DoFHandler using reinit() or "
2636  "in the constructor before you can distribute DoFs."));
2637  Assert(this->tria->n_levels() > 0,
2638  ExcMessage("The Triangulation you are using is empty!"));
2639  Assert(ff.size() > 0, ExcMessage("The hp::FECollection given is empty!"));
2640 
2641  //
2642  // register the new finite element collection
2643  //
2644  // don't create a new object if the one we have is identical
2645  if (this->fe_collection != ff)
2646  {
2648 
2649  const bool contains_multiple_fes = (this->fe_collection.size() > 1);
2650 
2651  // disable hp-mode if only a single finite element has been registered
2652  if (hp_capability_enabled && !contains_multiple_fes)
2653  {
2654  hp_capability_enabled = false;
2655 
2656  // unsubscribe connections to signals that are only relevant for
2657  // hp-mode, since we only have a single element here
2658  for (auto &connection : this->tria_listeners_for_transfer)
2659  connection.disconnect();
2660  this->tria_listeners_for_transfer.clear();
2661 
2662  // release active and future finite element tables
2663  this->hp_cell_active_fe_indices.clear();
2664  this->hp_cell_active_fe_indices.shrink_to_fit();
2665  this->hp_cell_future_fe_indices.clear();
2666  this->hp_cell_future_fe_indices.shrink_to_fit();
2667  }
2668 
2669  // re-enabling hp-mode is not permitted since the active and future FE
2670  // tables are no longer available
2671  AssertThrow(
2672  hp_capability_enabled || !contains_multiple_fes,
2673  ExcMessage(
2674  "You cannot re-enable hp-capabilities after you registered a single "
2675  "finite element. Please call reinit() or create a new DoFHandler "
2676  "object instead."));
2677  }
2678 
2679  //
2680  // enumerate all degrees of freedom
2681  //
2683  {
2684  // make sure every processor knows the active FE indices
2685  // on both its own cells and all ghost cells
2688 
2689 #ifdef DEBUG
2690  // make sure that the FE collection is large enough to
2691  // cover all FE indices presently in use on the mesh
2692  for (const auto &cell : this->active_cell_iterators())
2693  {
2694  if (!cell->is_artificial())
2695  Assert(cell->active_fe_index() < this->fe_collection.size(),
2696  ExcInvalidFEIndex(cell->active_fe_index(),
2697  this->fe_collection.size()));
2698  if (cell->is_locally_owned())
2699  Assert(cell->future_fe_index() < this->fe_collection.size(),
2700  ExcInvalidFEIndex(cell->future_fe_index(),
2701  this->fe_collection.size()));
2702  }
2703 #endif
2704  }
2705 
2706  {
2707  // We would like to enumerate all dofs for shared::Triangulations. If an
2708  // underlying shared::Tria allows artificial cells, we need to restore the
2709  // true cell owners temporarily.
2710  // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
2711  // the current set of subdomain ids, set subdomain ids to the "true" owner
2712  // of each cell upon construction of the TemporarilyRestoreSubdomainIds
2713  // object, and later restore these flags when it is destroyed.
2715  spacedim>
2716  subdomain_modifier(this->get_triangulation());
2717 
2718  // Adjust size of levels to the triangulation. Note that we still have to
2719  // allocate space for all degrees of freedom on this mesh (including ghost
2720  // and cells that are entirely stored on different processors), though we
2721  // may not assign numbers to some of them (i.e. they will remain at
2722  // invalid_dof_index). We need to allocate the space because we will want
2723  // to be able to query the dof_indices on each cell, and simply be told
2724  // that we don't know them on some cell (i.e. get back invalid_dof_index)
2727  *this);
2728  else
2730  }
2731 
2732  // hand the actual work over to the policy
2733  this->number_cache = this->policy->distribute_dofs();
2734 
2735  // do some housekeeping: compress indices
2736  // if(hp_capability_enabled)
2737  // {
2738  // Threads::TaskGroup<> tg;
2739  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2740  // tg += Threads::new_task(
2741  // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2742  // *this->levels_hp[level],
2743  // this->fe_collection);
2744  // tg.join_all();
2745  // }
2746 
2747  // Initialize the block info object only if this is a sequential
2748  // triangulation. It doesn't work correctly yet if it is parallel and has not
2749  // yet been implemented for hp-mode.
2750  if (!hp_capability_enabled &&
2752  *>(&*this->tria) == nullptr)
2753  this->block_info_object.initialize(*this, false, true);
2754 }
2755 
2756 
2757 
2758 template <int dim, int spacedim>
2759 void
2761 {
2763 
2764  Assert(
2765  this->object_dof_indices.size() > 0,
2766  ExcMessage(
2767  "Distribute active DoFs using distribute_dofs() before calling distribute_mg_dofs()."));
2768 
2769  Assert(
2770  ((this->tria->get_mesh_smoothing() &
2773  ExcMessage(
2774  "The mesh smoothing requirement 'limit_level_difference_at_vertices' has to be set for using multigrid!"));
2775 
2776  this->clear_mg_space();
2777 
2779  this->mg_number_cache = this->policy->distribute_mg_dofs();
2780 
2781  // initialize the block info object only if this is a sequential
2782  // triangulation. it doesn't work correctly yet if it is parallel
2783  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2784  &*this->tria) == nullptr)
2785  this->block_info_object.initialize(*this, true, false);
2786 }
2787 
2788 
2789 
2790 template <int dim, int spacedim>
2791 void
2793 {
2795 
2796  this->block_info_object.initialize_local(*this);
2797 }
2798 
2799 
2800 
2801 template <int dim, int spacedim>
2802 void
2804 {
2805  // decide whether we need a sequential or a parallel distributed policy
2806  if (dynamic_cast<const ::parallel::shared::Triangulation<dim, spacedim>
2807  *>(&this->get_triangulation()) != nullptr)
2808  this->policy = std::make_unique<internal::DoFHandlerImplementation::Policy::
2809  ParallelShared<dim, spacedim>>(*this);
2810  else if (dynamic_cast<
2811  const ::parallel::DistributedTriangulationBase<dim, spacedim>
2812  *>(&this->get_triangulation()) == nullptr)
2813  this->policy = std::make_unique<
2815  *this);
2816  else
2817  this->policy =
2818  std::make_unique<internal::DoFHandlerImplementation::Policy::
2819  ParallelDistributed<dim, spacedim>>(*this);
2820 }
2821 
2822 
2823 
2824 template <int dim, int spacedim>
2825 void
2827 {
2828  // release memory
2829  this->clear_space();
2830  this->clear_mg_space();
2831 }
2832 
2833 
2834 
2835 template <int dim, int spacedim>
2836 void
2838 {
2839  cell_dof_cache_indices.clear();
2840 
2841  cell_dof_cache_ptr.clear();
2842 
2843  object_dof_indices.clear();
2844 
2845  object_dof_ptr.clear();
2846 
2847  this->number_cache.clear();
2848 
2849  this->hp_cell_active_fe_indices.clear();
2850  this->hp_cell_future_fe_indices.clear();
2851 }
2852 
2853 
2854 
2855 template <int dim, int spacedim>
2856 void
2858 {
2859  this->mg_levels.clear();
2860  this->mg_faces.reset();
2861 
2862  std::vector<MGVertexDoFs> tmp;
2863 
2864  std::swap(this->mg_vertex_dofs, tmp);
2865 
2866  this->mg_number_cache.clear();
2867 }
2868 
2869 
2870 
2871 template <int dim, int spacedim>
2872 void
2874  const std::vector<types::global_dof_index> &new_numbers)
2875 {
2877  {
2878  Assert(this->hp_cell_future_fe_indices.size() > 0,
2879  ExcMessage(
2880  "You need to distribute DoFs before you can renumber them."));
2881 
2882  AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2883 
2884 #ifdef DEBUG
2885  // assert that the new indices are consecutively numbered if we are
2886  // working on a single processor. this doesn't need to
2887  // hold in the case of a parallel mesh since we map the interval
2888  // [0...n_dofs()) into itself but only globally, not on each processor
2889  if (this->n_locally_owned_dofs() == this->n_dofs())
2890  {
2891  std::vector<types::global_dof_index> tmp(new_numbers);
2892  std::sort(tmp.begin(), tmp.end());
2893  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2895  for (; p != tmp.end(); ++p, ++i)
2896  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2897  }
2898  else
2899  for (const auto new_number : new_numbers)
2900  Assert(new_number < this->n_dofs(),
2901  ExcMessage(
2902  "New DoF index is not less than the total number of dofs."));
2903 #endif
2904 
2905  // uncompress the internal storage scheme of dofs on cells so that
2906  // we can access dofs in turns. uncompress in parallel, starting
2907  // with the most expensive levels (the highest ones)
2908  //{
2909  // Threads::TaskGroup<> tg;
2910  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2911  // tg += Threads::new_task(
2912  // &::internal::hp::DoFLevel::uncompress_data<dim, spacedim>,
2913  // *this->levels_hp[level],
2914  // this->fe_collection);
2915  // tg.join_all();
2916  //}
2917 
2918  // do the renumbering
2919  this->number_cache = this->policy->renumber_dofs(new_numbers);
2920 
2921  // now re-compress the dof indices
2922  //{
2923  // Threads::TaskGroup<> tg;
2924  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2925  // tg += Threads::new_task(
2926  // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2927  // *this->levels_hp[level],
2928  // this->fe_collection);
2929  // tg.join_all();
2930  //}
2931  }
2932  else
2933  {
2934  Assert(this->object_dof_indices.size() > 0,
2935  ExcMessage(
2936  "You need to distribute DoFs before you can renumber them."));
2937 
2938 #ifdef DEBUG
2939  if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
2940  &*this->tria) != nullptr)
2941  {
2942  Assert(new_numbers.size() == this->n_dofs() ||
2943  new_numbers.size() == this->n_locally_owned_dofs(),
2944  ExcMessage("Incorrect size of the input array."));
2945  }
2946  else if (dynamic_cast<
2948  &*this->tria) != nullptr)
2949  {
2950  AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2951  }
2952  else
2953  {
2954  AssertDimension(new_numbers.size(), this->n_dofs());
2955  }
2956 
2957  // assert that the new indices are consecutively numbered if we are
2958  // working on a single processor. this doesn't need to
2959  // hold in the case of a parallel mesh since we map the interval
2960  // [0...n_dofs()) into itself but only globally, not on each processor
2961  if (this->n_locally_owned_dofs() == this->n_dofs())
2962  {
2963  std::vector<types::global_dof_index> tmp(new_numbers);
2964  std::sort(tmp.begin(), tmp.end());
2965  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2967  for (; p != tmp.end(); ++p, ++i)
2968  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2969  }
2970  else
2971  for (const auto new_number : new_numbers)
2972  Assert(new_number < this->n_dofs(),
2973  ExcMessage(
2974  "New DoF index is not less than the total number of dofs."));
2975 #endif
2976 
2977  this->number_cache = this->policy->renumber_dofs(new_numbers);
2978  }
2979 }
2980 
2981 
2982 
2983 template <int dim, int spacedim>
2984 void
2986  const unsigned int level,
2987  const std::vector<types::global_dof_index> &new_numbers)
2988 {
2990 
2991  Assert(
2992  this->mg_levels.size() > 0 && this->object_dof_indices.size() > 0,
2993  ExcMessage(
2994  "You need to distribute active and level DoFs before you can renumber level DoFs."));
2995  AssertIndexRange(level, this->get_triangulation().n_global_levels());
2996  AssertDimension(new_numbers.size(),
2997  this->locally_owned_mg_dofs(level).n_elements());
2998 
2999 #ifdef DEBUG
3000  // assert that the new indices are consecutively numbered if we are working
3001  // on a single processor. this doesn't need to hold in the case of a
3002  // parallel mesh since we map the interval [0...n_dofs(level)) into itself
3003  // but only globally, not on each processor
3004  if (this->n_locally_owned_dofs() == this->n_dofs())
3005  {
3006  std::vector<types::global_dof_index> tmp(new_numbers);
3007  std::sort(tmp.begin(), tmp.end());
3008  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
3010  for (; p != tmp.end(); ++p, ++i)
3011  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
3012  }
3013  else
3014  for (const auto new_number : new_numbers)
3015  Assert(new_number < this->n_dofs(level),
3016  ExcMessage(
3017  "New DoF index is not less than the total number of dofs."));
3018 #endif
3019 
3020  this->mg_number_cache[level] =
3021  this->policy->renumber_mg_dofs(level, new_numbers);
3022 }
3023 
3024 
3025 
3026 template <int dim, int spacedim>
3027 unsigned int
3029 {
3030  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
3031 
3032  switch (dim)
3033  {
3034  case 1:
3035  return this->fe_collection.max_dofs_per_vertex();
3036  case 2:
3037  return (3 * this->fe_collection.max_dofs_per_vertex() +
3038  2 * this->fe_collection.max_dofs_per_line());
3039  case 3:
3040  // we need to take refinement of one boundary face into
3041  // consideration here; in fact, this function returns what
3042  // #max_coupling_between_dofs<2> returns
3043  //
3044  // we assume here, that only four faces meet at the boundary;
3045  // this assumption is not justified and needs to be fixed some
3046  // time. fortunately, omitting it for now does no harm since
3047  // the matrix will cry foul if its requirements are not
3048  // satisfied
3049  return (19 * this->fe_collection.max_dofs_per_vertex() +
3050  28 * this->fe_collection.max_dofs_per_line() +
3051  8 * this->fe_collection.max_dofs_per_quad());
3052  default:
3053  Assert(false, ExcNotImplemented());
3054  return 0;
3055  }
3056 }
3057 
3058 
3059 
3060 template <int dim, int spacedim>
3061 unsigned int
3063 {
3064  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
3067 }
3068 
3069 
3070 
3071 template <int dim, int spacedim>
3072 template <int structdim>
3074 DoFHandler<dim, spacedim>::get_dof_index(const unsigned int obj_level,
3075  const unsigned int obj_index,
3076  const unsigned int fe_index,
3077  const unsigned int local_index) const
3078 {
3080  {
3081  Assert(false, ExcNotImplemented());
3083  }
3084  else
3085  {
3087  *this,
3088  this->mg_levels[obj_level],
3089  this->mg_faces,
3090  obj_index,
3091  fe_index,
3092  local_index,
3093  std::integral_constant<int, structdim>());
3094  }
3095 }
3096 
3097 
3098 
3099 template <int dim, int spacedim>
3100 template <int structdim>
3101 void
3103  const unsigned int obj_level,
3104  const unsigned int obj_index,
3105  const unsigned int fe_index,
3106  const unsigned int local_index,
3108 {
3110  {
3111  Assert(false, ExcNotImplemented());
3112  return;
3113  }
3114  else
3115  {
3117  *this,
3118  this->mg_levels[obj_level],
3119  this->mg_faces,
3120  obj_index,
3121  fe_index,
3122  local_index,
3123  global_index,
3124  std::integral_constant<int, structdim>());
3125  }
3126 }
3127 
3128 
3129 
3130 template <int dim, int spacedim>
3131 void
3133  const std::vector<unsigned int> &active_fe_indices)
3134 {
3135  Assert(active_fe_indices.size() == this->get_triangulation().n_active_cells(),
3136  ExcDimensionMismatch(active_fe_indices.size(),
3137  this->get_triangulation().n_active_cells()));
3138 
3139  this->create_active_fe_table();
3140  // we could set the values directly, since they are stored as
3141  // protected data of this object, but for simplicity we use the
3142  // cell-wise access. this way we also have to pass some debug-mode
3143  // tests which we would have to duplicate ourselves otherwise
3144  for (const auto &cell : this->active_cell_iterators())
3145  if (cell->is_locally_owned())
3146  cell->set_active_fe_index(active_fe_indices[cell->active_cell_index()]);
3147 }
3148 
3149 
3150 
3151 template <int dim, int spacedim>
3152 void
3154  std::vector<unsigned int> &active_fe_indices) const
3155 {
3156  active_fe_indices.resize(this->get_triangulation().n_active_cells());
3157 
3158  // we could try to extract the values directly, since they are
3159  // stored as protected data of this object, but for simplicity we
3160  // use the cell-wise access.
3161  for (const auto &cell : this->active_cell_iterators())
3162  if (!cell->is_artificial())
3163  active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
3164 }
3165 
3166 
3167 
3168 template <int dim, int spacedim>
3169 void
3171 {
3172  // make sure this is called during initialization in hp-mode
3174 
3175  // connect functions to signals of the underlying triangulation
3176  this->tria_listeners.push_back(this->tria->signals.create.connect(
3177  [this]() { this->reinit(*(this->tria)); }));
3178  this->tria_listeners.push_back(
3179  this->tria->signals.clear.connect([this]() { this->clear(); }));
3180 
3181  // attach corresponding callback functions dealing with the transfer of
3182  // active FE indices depending on the type of triangulation
3183  if (dynamic_cast<
3184  const ::parallel::fullydistributed::Triangulation<dim, spacedim>
3185  *>(&this->get_triangulation()))
3186  {
3187  // no transfer of active FE indices for this class
3188  }
3189  else if (dynamic_cast<
3190  const ::parallel::distributed::Triangulation<dim, spacedim>
3191  *>(&this->get_triangulation()))
3192  {
3193  // repartitioning signals
3194  this->tria_listeners_for_transfer.push_back(
3195  this->tria->signals.pre_distributed_repartition.connect([this]() {
3196  internal::hp::DoFHandlerImplementation::Implementation::
3197  ensure_absence_of_future_fe_indices<dim, spacedim>(*this);
3198  }));
3199  this->tria_listeners_for_transfer.push_back(
3200  this->tria->signals.pre_distributed_repartition.connect(
3201  [this]() { this->pre_distributed_transfer_action(); }));
3202  this->tria_listeners_for_transfer.push_back(
3203  this->tria->signals.post_distributed_repartition.connect(
3204  [this] { this->post_distributed_transfer_action(); }));
3205 
3206  // refinement signals
3207  this->tria_listeners_for_transfer.push_back(
3208  this->tria->signals.post_p4est_refinement.connect(
3209  [this]() { this->pre_distributed_transfer_action(); }));
3210  this->tria_listeners_for_transfer.push_back(
3211  this->tria->signals.post_distributed_refinement.connect(
3212  [this]() { this->post_distributed_transfer_action(); }));
3213 
3214  // serialization signals
3215  this->tria_listeners_for_transfer.push_back(
3216  this->tria->signals.post_distributed_save.connect(
3217  [this]() { this->active_fe_index_transfer.reset(); }));
3218  this->tria_listeners_for_transfer.push_back(
3219  this->tria->signals.post_distributed_load.connect(
3220  [this]() { this->update_active_fe_table(); }));
3221  }
3222  else if (dynamic_cast<
3223  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3224  &this->get_triangulation()) != nullptr)
3225  {
3226  // partitioning signals
3227  this->tria_listeners_for_transfer.push_back(
3228  this->tria->signals.pre_partition.connect([this]() {
3229  internal::hp::DoFHandlerImplementation::Implementation::
3230  ensure_absence_of_future_fe_indices(*this);
3231  }));
3232 
3233  // refinement signals
3234  this->tria_listeners_for_transfer.push_back(
3235  this->tria->signals.pre_refinement.connect([this]() {
3236  internal::hp::DoFHandlerImplementation::Implementation::
3237  communicate_future_fe_indices(*this);
3238  }));
3239  this->tria_listeners_for_transfer.push_back(
3240  this->tria->signals.pre_refinement.connect(
3241  [this] { this->pre_transfer_action(); }));
3242  this->tria_listeners_for_transfer.push_back(
3243  this->tria->signals.post_refinement.connect(
3244  [this] { this->post_transfer_action(); }));
3245  }
3246  else
3247  {
3248  // refinement signals
3249  this->tria_listeners_for_transfer.push_back(
3250  this->tria->signals.pre_refinement.connect(
3251  [this] { this->pre_transfer_action(); }));
3252  this->tria_listeners_for_transfer.push_back(
3253  this->tria->signals.post_refinement.connect(
3254  [this] { this->post_transfer_action(); }));
3255  }
3256 }
3257 
3258 
3259 
3260 template <int dim, int spacedim>
3261 void
3263 {
3265 
3266 
3267  // Create sufficiently many hp::DoFLevels.
3268  // while (this->levels_hp.size() < this->tria->n_levels())
3269  // this->levels_hp.emplace_back(new ::internal::hp::DoFLevel);
3270 
3271  this->hp_cell_active_fe_indices.resize(this->tria->n_levels());
3272  this->hp_cell_future_fe_indices.resize(this->tria->n_levels());
3273 
3274  // then make sure that on each level we have the appropriate size
3275  // of active FE indices; preset them to zero, i.e. the default FE
3276  for (unsigned int level = 0; level < this->hp_cell_future_fe_indices.size();
3277  ++level)
3278  {
3279  if (this->hp_cell_active_fe_indices[level].size() == 0 &&
3280  this->hp_cell_future_fe_indices[level].size() == 0)
3281  {
3282  this->hp_cell_active_fe_indices[level].resize(
3283  this->tria->n_raw_cells(level), 0);
3284  this->hp_cell_future_fe_indices[level].resize(
3285  this->tria->n_raw_cells(level), invalid_active_fe_index);
3286  }
3287  else
3288  {
3289  // Either the active FE indices have size zero because
3290  // they were just created, or the correct size. Other
3291  // sizes indicate that something went wrong.
3292  Assert(this->hp_cell_active_fe_indices[level].size() ==
3293  this->tria->n_raw_cells(level) &&
3294  this->hp_cell_future_fe_indices[level].size() ==
3295  this->tria->n_raw_cells(level),
3296  ExcInternalError());
3297  }
3298 
3299  // it may be that the previous table was compressed; in that
3300  // case, restore the correct active FE index. the fact that
3301  // this no longer matches the indices in the table is of no
3302  // importance because the current function is called at a
3303  // point where we have to recreate the dof_indices tables in
3304  // the levels anyway
3305  // this->levels_hp[level]->normalize_active_fe_indices();
3306  }
3307 }
3308 
3309 
3310 
3311 template <int dim, int spacedim>
3312 void
3314 {
3315  // // Normally only one level is added, but if this Triangulation
3316  // // is created by copy_triangulation, it can be more than one level.
3317  // while (this->levels_hp.size() < this->tria->n_levels())
3318  // this->levels_hp.emplace_back(new ::internal::hp::DoFLevel);
3319  //
3320  // // Coarsening can lead to the loss of levels. Hence remove them.
3321  // while (this->levels_hp.size() > this->tria->n_levels())
3322  // {
3323  // // drop the last element. that also releases the memory pointed to
3324  // this->levels_hp.pop_back();
3325  // }
3326 
3327  this->hp_cell_active_fe_indices.resize(this->tria->n_levels());
3328  this->hp_cell_active_fe_indices.shrink_to_fit();
3329 
3330  this->hp_cell_future_fe_indices.resize(this->tria->n_levels());
3331  this->hp_cell_future_fe_indices.shrink_to_fit();
3332 
3333  for (unsigned int i = 0; i < this->hp_cell_future_fe_indices.size(); ++i)
3334  {
3335  // Resize active FE indices vectors. Use zero indicator to extend.
3336  this->hp_cell_active_fe_indices[i].resize(this->tria->n_raw_cells(i), 0);
3337 
3338  // Resize future FE indices vectors. Make sure that all
3339  // future FE indices have been cleared after refinement happened.
3340  //
3341  // We have used future FE indices to update all active FE indices
3342  // before refinement happened, thus we are safe to clear them now.
3343  this->hp_cell_future_fe_indices[i].assign(this->tria->n_raw_cells(i),
3345  }
3346 }
3347 
3348 
3349 template <int dim, int spacedim>
3350 void
3352 {
3353  Assert(this->active_fe_index_transfer == nullptr, ExcInternalError());
3354 
3355  this->active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3356 
3359 }
3360 
3361 
3362 
3363 template <int dim, int spacedim>
3364 void
3366 {
3367 #ifndef DEAL_II_WITH_P4EST
3368  Assert(false,
3369  ExcMessage(
3370  "You are attempting to use a functionality that is only available "
3371  "if deal.II was configured to use p4est, but cmake did not find a "
3372  "valid p4est library."));
3373 #else
3374  // the implementation below requires a p:d:T currently
3375  Assert(
3377  &this->get_triangulation()) != nullptr),
3378  ExcNotImplemented());
3379 
3381 
3382  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3383 
3384  // If we work on a p::d::Triangulation, we have to transfer all
3385  // active FE indices since ownership of cells may change. We will
3386  // use our p::d::CellDataTransfer member to achieve this. Further,
3387  // we prepare the values in such a way that they will correspond to
3388  // the active FE indices on the new mesh.
3389 
3390  // Gather all current future FE indices.
3391  active_fe_index_transfer->active_fe_indices.resize(
3393 
3394  for (const auto &cell : active_cell_iterators())
3395  if (cell->is_locally_owned())
3396  active_fe_index_transfer->active_fe_indices[cell->active_cell_index()] =
3397  cell->future_fe_index();
3398 
3399  // Create transfer object and attach to it.
3400  const auto *distributed_tria =
3402  &this->get_triangulation());
3403 
3404  active_fe_index_transfer->cell_data_transfer = std::make_unique<
3405  parallel::distributed::
3406  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3407  *distributed_tria,
3408  /*transfer_variable_size_data=*/false,
3409  /*refinement_strategy=*/
3410  &::AdaptationStrategies::Refinement::
3411  preserve<dim, spacedim, unsigned int>,
3412  /*coarsening_strategy=*/
3413  [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3414  const std::vector<unsigned int> &children_fe_indices)
3415  -> unsigned int {
3416  return ::internal::hp::DoFHandlerImplementation::Implementation::
3417  determine_fe_from_children<dim, spacedim>(parent,
3418  children_fe_indices,
3419  fe_collection);
3420  });
3421 
3422  active_fe_index_transfer->cell_data_transfer
3423  ->prepare_for_coarsening_and_refinement(
3424  active_fe_index_transfer->active_fe_indices);
3425 #endif
3426 }
3427 
3428 
3429 
3430 template <int dim, int spacedim>
3431 void
3433 {
3435 
3436  Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
3437 
3440 
3441  // We have to distribute the information about active FE indices
3442  // of all cells (including the artificial ones) on all processors,
3443  // if a parallel::shared::Triangulation has been used.
3446 
3447  // Free memory.
3448  this->active_fe_index_transfer.reset();
3449 }
3450 
3451 
3452 
3453 template <int dim, int spacedim>
3454 void
3456 {
3457 #ifndef DEAL_II_WITH_P4EST
3458  Assert(false, ExcInternalError());
3459 #else
3461 
3462  Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
3463 
3464  // Unpack active FE indices.
3465  this->active_fe_index_transfer->active_fe_indices.resize(
3467  this->active_fe_index_transfer->cell_data_transfer->unpack(
3468  this->active_fe_index_transfer->active_fe_indices);
3469 
3470  // Update all locally owned active FE indices.
3471  this->set_active_fe_indices(
3472  this->active_fe_index_transfer->active_fe_indices);
3473 
3474  // Update active FE indices on ghost cells.
3477 
3478  // Free memory.
3479  this->active_fe_index_transfer.reset();
3480 #endif
3481 }
3482 
3483 
3484 
3485 template <int dim, int spacedim>
3486 void
3488 {
3489 #ifndef DEAL_II_WITH_P4EST
3490  Assert(false,
3491  ExcMessage(
3492  "You are attempting to use a functionality that is only available "
3493  "if deal.II was configured to use p4est, but cmake did not find a "
3494  "valid p4est library."));
3495 #else
3496  // the implementation below requires a p:d:T currently
3497  Assert(
3499  &this->get_triangulation()) != nullptr),
3500  ExcNotImplemented());
3501 
3503 
3504  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3505 
3506  // Create transfer object and attach to it.
3507  const auto *distributed_tria =
3509  &this->get_triangulation());
3510 
3511  active_fe_index_transfer->cell_data_transfer = std::make_unique<
3512  parallel::distributed::
3513  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3514  *distributed_tria,
3515  /*transfer_variable_size_data=*/false,
3516  /*refinement_strategy=*/
3517  &::AdaptationStrategies::Refinement::
3518  preserve<dim, spacedim, unsigned int>,
3519  /*coarsening_strategy=*/
3520  [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3521  const std::vector<unsigned int> &children_fe_indices)
3522  -> unsigned int {
3523  return ::internal::hp::DoFHandlerImplementation::Implementation::
3524  determine_fe_from_children<dim, spacedim>(parent,
3525  children_fe_indices,
3526  fe_collection);
3527  });
3528 
3529  // If we work on a p::d::Triangulation, we have to transfer all
3530  // active FE indices since ownership of cells may change.
3531 
3532  // Gather all current active FE indices
3533  get_active_fe_indices(active_fe_index_transfer->active_fe_indices);
3534 
3535  // Attach to transfer object
3536  active_fe_index_transfer->cell_data_transfer->prepare_for_serialization(
3537  active_fe_index_transfer->active_fe_indices);
3538 #endif
3539 }
3540 
3541 
3542 
3543 template <int dim, int spacedim>
3544 void
3546 {
3547 #ifndef DEAL_II_WITH_P4EST
3548  Assert(false,
3549  ExcMessage(
3550  "You are attempting to use a functionality that is only available "
3551  "if deal.II was configured to use p4est, but cmake did not find a "
3552  "valid p4est library."));
3553 #else
3554  // the implementation below requires a p:d:T currently
3555  Assert(
3557  &this->get_triangulation()) != nullptr),
3558  ExcNotImplemented());
3559 
3561 
3562  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3563 
3564  // Create transfer object and attach to it.
3565  const auto *distributed_tria =
3567  &this->get_triangulation());
3568 
3569  active_fe_index_transfer->cell_data_transfer = std::make_unique<
3570  parallel::distributed::
3571  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3572  *distributed_tria,
3573  /*transfer_variable_size_data=*/false,
3574  /*refinement_strategy=*/
3575  &::AdaptationStrategies::Refinement::
3576  preserve<dim, spacedim, unsigned int>,
3577  /*coarsening_strategy=*/
3578  [this](const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3579  const std::vector<unsigned int> &children_fe_indices)
3580  -> unsigned int {
3581  return ::internal::hp::DoFHandlerImplementation::Implementation::
3582  determine_fe_from_children<dim, spacedim>(parent,
3583  children_fe_indices,
3584  fe_collection);
3585  });
3586 
3587  // Unpack active FE indices.
3588  active_fe_index_transfer->active_fe_indices.resize(
3590  active_fe_index_transfer->cell_data_transfer->deserialize(
3591  active_fe_index_transfer->active_fe_indices);
3592 
3593  // Update all locally owned active FE indices.
3594  set_active_fe_indices(active_fe_index_transfer->active_fe_indices);
3595 
3596  // Update active FE indices on ghost cells.
3599 
3600  // Free memory.
3601  active_fe_index_transfer.reset();
3602 #endif
3603 }
3604 
3605 
3606 
3607 template <int dim, int spacedim>
3609  : coarsest_level(numbers::invalid_unsigned_int)
3610  , finest_level(0)
3611 {}
3612 
3613 
3614 
3615 template <int dim, int spacedim>
3616 void
3618  const unsigned int cl,
3619  const unsigned int fl,
3620  const unsigned int dofs_per_vertex)
3621 {
3622  coarsest_level = cl;
3623  finest_level = fl;
3624 
3626  {
3627  const unsigned int n_levels = finest_level - coarsest_level + 1;
3628  const unsigned int n_indices = n_levels * dofs_per_vertex;
3629 
3630  indices = std::make_unique<types::global_dof_index[]>(n_indices);
3631  std::fill(indices.get(),
3632  indices.get() + n_indices,
3634  }
3635  else
3636  indices.reset();
3637 }
3638 
3639 
3640 
3641 template <int dim, int spacedim>
3642 unsigned int
3644 {
3645  return coarsest_level;
3646 }
3647 
3648 
3649 
3650 template <int dim, int spacedim>
3651 unsigned int
3653 {
3654  return finest_level;
3655 }
3656 
3657 /*-------------- Explicit Instantiations -------------------------------*/
3658 #include "dof_handler.inst"
3659 
3660 
3661 
IteratorRange< active_cell_iterator > active_cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcInvalidFEIndex(int arg1, int arg2)
unsigned int max_couplings_between_dofs() const
void connect_to_triangulation_signals()
unsigned int get_coarsest_level() const
static void reserve_space_mg(DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:518
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1587
const Triangulation< dim, spacedim > & get_triangulation() const
level_cell_iterator end_mg() const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
std::array< std::vector< offset_type >, dim+1 > hp_object_fe_ptr
Definition: dof_handler.h:1567
static void reserve_cells(DoFHandler< dim, spacedim > &dof_handler, const unsigned int n_inner_dofs_per_cell)
Definition: dof_handler.cc:290
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
static void reserve_space_mg(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:448
cell_iterator begin(const unsigned int level=0) const
static ::ExceptionBase & ExcNoDominatedFiniteElementOnChildren()
void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
static void set_dof_index(const DoFHandler< 3, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 3 >> &, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 3 >> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const std::integral_constant< int, 1 >)
Definition: dof_handler.cc:881
unsigned short int active_fe_index_type
Definition: dof_handler.h:529
Task< RT > new_task(const std::function< RT()> &function)
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})
static types::global_dof_index get_dof_index(const DoFHandler< 3, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 3 >> &, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 3 >> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const std::integral_constant< int, 1 >)
Definition: dof_handler.cc:743
static ::ExceptionBase & ExcOnlyAvailableWithHP()
void pre_transfer_action()
const ::Triangulation< dim, spacedim > & tria
cell_iterator end() const
::internal::DoFHandlerImplementation::NumberCache number_cache
Definition: dof_handler.h:1515
static void reserve_space(::DoFHandler< 3, spacedim > &dof_handler)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
void clear()
std::vector< std::array< std::vector< offset_type >, dim+1 > > object_dof_ptr
Definition: dof_handler.h:1554
static unsigned int max_couplings_between_dofs(const DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:101
static unsigned int max_couplings_between_dofs(const DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:111
IteratorRange< level_cell_iterator > mg_cell_iterators_on_level(const unsigned int level) const
void deserialize_active_fe_indices()
void update_active_fe_table()
void clear_mg_space()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
virtual std::size_t memory_consumption() const
static void reserve_space_cells(DoFHandler< dim, spacedim > &dof_handler)
void initialize(const Triangulation< dim, spacedim > &tria, const FiniteElement< dim, spacedim > &fe)
void setup_policy()
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition: dof_handler.h:1594
static ::ExceptionBase & ExcNoFESelected()
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
static void distribute_fe_indices_on_refined_cells(DoFHandler< dim, spacedim > &dof_handler)
std::vector< std::vector< active_fe_index_type > > hp_cell_future_fe_indices
Definition: dof_handler.h:1581
static void set_dof_index(const DoFHandler< 2, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 2 >> &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 2 >> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const std::integral_constant< int, 2 >)
Definition: dof_handler.cc:857
static void communicate_active_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
std::vector< std::vector< types::global_dof_index > > cell_dof_cache_indices
Definition: dof_handler.h:1529
static void set_dof_index(const DoFHandler< 3, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 3 >> &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 3 >> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const std::integral_constant< int, 3 >)
Definition: dof_handler.cc:929
void set_dof_index(const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const
void post_distributed_transfer_action()
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
active_cell_iterator begin_active(const unsigned int level=0) const
static void communicate_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
types::global_dof_index n_locally_owned_dofs() const
BlockInfo block_info_object
Definition: dof_handler.h:1479
static ::ExceptionBase & ExcMessage(std::string arg1)
static void reserve_space(::DoFHandler< 1, spacedim > &dof_handler)
static unsigned int dominated_future_fe_on_children(const typename DoFHandler< dim, spacedim >::cell_iterator &parent)
static void set_dof_index(const DoFHandler< 2, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 2 >> &, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 2 >> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const std::integral_constant< int, 1 >)
Definition: dof_handler.cc:833
static types::global_dof_index get_dof_index(const DoFHandler< 2, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 2 >> &, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 2 >> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const std::integral_constant< int, 1 >)
Definition: dof_handler.cc:701
unsigned int dominated_future_fe_on_children(const typename DoFHandler< dim, spacedim >::cell_iterator &parent)
static const char T
T sum(const T &t, const MPI_Comm &mpi_communicator)
std::vector<::internal::DoFHandlerImplementation::NumberCache > mg_number_cache
Definition: dof_handler.h:1521
void distribute_mg_dofs()
void post_transfer_action()
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:339
#define Assert(cond, exc)
Definition: exceptions.h:1461
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1498
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::unique_ptr< ActiveFEIndexTransfer > active_fe_index_transfer
Definition: dof_handler.h:1606
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition: dof_handler.h:1543
active_cell_iterator end_active(const unsigned int level) const
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
Definition: dof_handler.h:1600
static unsigned int determine_fe_from_children(const typename Triangulation< dim, spacedim >::cell_iterator &, const std::vector< unsigned int > &children_fe_indices, const ::hp::FECollection< dim, spacedim > &fe_collection)
types::global_dof_index n_boundary_dofs() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:1491
types::global_dof_index n_dofs() const
void pre_distributed_transfer_action()
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:502
level_cell_iterator begin_mg(const unsigned int level=0) const
std::unique_ptr<::internal::DoFHandlerImplementation::Policy::PolicyBase< dim, spacedim > > policy
Definition: dof_handler.h:1506
static unsigned int max_couplings_between_dofs(const DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:230
static types::global_dof_index get_dof_index(const DoFHandler< 3, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 3 >> &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 3 >> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const std::integral_constant< int, 3 >)
Definition: dof_handler.cc:787
const hp::FECollection< dim, spacedim > & get_fe_collection() const
void set_fe(const FiniteElement< dim, spacedim > &fe)
void coarsen(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold)
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
static void reserve_space_mg(DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:595
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13704
unsigned int get_finest_level() const
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1219
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Definition: hp.h:117
void set_active_fe_indices(const std::vector< unsigned int > &active_fe_indices)
std::vector< std::vector< active_fe_index_type > > hp_cell_active_fe_indices
Definition: dof_handler.h:1574
std::array< std::vector< active_fe_index_type >, dim+1 > hp_object_fe_indices
Definition: dof_handler.h:1562
void initialize(const DoFHandler< dim, spacedim > &, bool levels_only=false, bool active_only=false)
Fill the object with values describing block structure of the DoFHandler.
Definition: block_info.cc:30
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
std::unique_ptr< types::global_dof_index[]> indices
Definition: dof_handler.h:1429
void initialize_local_block_info()
std::string policy_to_string(const ::internal::DoFHandlerImplementation::Policy::PolicyBase< dim, spacedim > &policy)
Definition: dof_handler.cc:56
static const active_fe_index_type invalid_active_fe_index
Definition: dof_handler.h:540
void reinit(const Triangulation< dim, spacedim > &tria)
static types::global_dof_index get_dof_index(const DoFHandler< 3, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 3 >> &, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 3 >> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const std::integral_constant< int, 2 >)
Definition: dof_handler.cc:765
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
void initialize_local(const DoFHandler< dim, spacedim > &)
Initialize block structure on cells and compute renumbering between cell dofs and block cell dofs...
Definition: block_info.cc:62
IteratorRange< level_cell_iterator > mg_cell_iterators() const
unsigned int coarsest_level
Definition: dof_handler.h:1413
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
static void set_dof_index(const DoFHandler< 3, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 3 >> &, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 3 >> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const std::integral_constant< int, 2 >)
Definition: dof_handler.cc:905
void clear_space()
static ::ExceptionBase & ExcNotImplementedWithHP()
static void ensure_absence_of_future_fe_indices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:972
void init(const unsigned int coarsest_level, const unsigned int finest_level, const unsigned int dofs_per_vertex)
static void reset_to_empty_objects(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:266
static void collect_fe_indices_on_cells_to_be_refined(DoFHandler< dim, spacedim > &dof_handler)
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
static void reserve_space(::DoFHandler< 2, spacedim > &dof_handler)
types::global_dof_index get_dof_index(const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition: utilities.h:1322
static void reserve_space_vertices(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:992
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcNotImplemented()
Iterator points to a valid object.
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
bool is_element(const size_type index) const
Definition: index_set.h:1770
static ::ExceptionBase & ExcNewNumbersNotConsecutive(types::global_dof_index arg1)
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
static void set_dof_index(const DoFHandler< 1, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 1 >> &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 1 >> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const std::integral_constant< int, 1 >)
Definition: dof_handler.cc:809
const types::boundary_id internal_face_boundary_id
Definition: types.h:255
unsigned int max_couplings_between_boundary_dofs() const
void get_active_fe_indices(std::vector< unsigned int > &active_fe_indices) const
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
const types::global_dof_index invalid_dof_index
Definition: types.h:211
IteratorState::IteratorStates state() const
virtual ~DoFHandler() override
std::vector< boost::signals2::connection > tria_listeners
Definition: dof_handler.h:1612
std::vector< boost::signals2::connection > tria_listeners_for_transfer
Definition: dof_handler.h:1619
const IndexSet & locally_owned_dofs() const
static types::global_dof_index get_dof_index(const DoFHandler< 1, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 1 >> &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 1 >> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const std::integral_constant< int, 1 >)
Definition: dof_handler.cc:678
bool has_level_dofs() const
static types::global_dof_index get_dof_index(const DoFHandler< 2, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< 2 >> &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< 2 >> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const std::integral_constant< int, 2 >)
Definition: dof_handler.cc:721
unsigned int boundary_id
Definition: types.h:129
size_type n_elements() const
Definition: index_set.h:1837
static ::ExceptionBase & ExcNoFESelected()
IteratorRange< cell_iterator > cell_iterators() const
static void reserve_space_faces(DoFHandler< dim, spacedim > &dof_handler)
bool hp_capability_enabled
Definition: dof_handler.h:1485
void prepare_for_serialization_of_active_fe_indices()
static void reserve_space(DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_handler.cc:398
std::vector< std::vector< offset_type > > cell_dof_cache_ptr
Definition: dof_handler.h:1535
void create_active_fe_table()
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()