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