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