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