Reference documentation for deal.II version Git d2d482dcec 2020-10-26 10:29:03 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_handler.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 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 
259  template <int spacedim>
260  static void reserve_space(DoFHandler<1, spacedim> &dof_handler)
261  {
262  dof_handler.object_dof_indices[0][0].resize(
263  dof_handler.tria->n_vertices() *
264  dof_handler.get_fe().n_dofs_per_vertex(),
266 
267  for (unsigned int i = 0; i < dof_handler.tria->n_levels(); ++i)
268  {
269  dof_handler.object_dof_indices[i][1].resize(
270  dof_handler.tria->n_raw_cells(i) *
271  dof_handler.get_fe().n_dofs_per_line(),
273 
274  dof_handler.object_dof_ptr[i][1].reserve(
275  dof_handler.tria->n_raw_cells(i) + 1);
276  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1;
277  j++)
278  dof_handler.object_dof_ptr[i][1].push_back(
279  j * dof_handler.get_fe().n_dofs_per_line());
280 
281  dof_handler.cell_dof_cache_indices[i].resize(
282  dof_handler.tria->n_raw_cells(i) *
283  dof_handler.get_fe().n_dofs_per_cell(),
285 
286  dof_handler.cell_dof_cache_ptr[i].reserve(
287  dof_handler.tria->n_raw_cells(i) + 1);
288  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1;
289  j++)
290  dof_handler.cell_dof_cache_ptr[i].push_back(
291  j * dof_handler.get_fe().n_dofs_per_cell());
292  }
293 
294  dof_handler.object_dof_indices[0][0].resize(
295  dof_handler.tria->n_vertices() *
296  dof_handler.get_fe().n_dofs_per_vertex(),
298  }
299 
300  template <int spacedim>
301  static void reserve_space(DoFHandler<2, spacedim> &dof_handler)
302  {
303  dof_handler.object_dof_indices[0][0].resize(
304  dof_handler.tria->n_vertices() *
305  dof_handler.get_fe().n_dofs_per_vertex(),
307 
308  for (unsigned int i = 0; i < dof_handler.tria->n_levels(); ++i)
309  {
310  dof_handler.object_dof_indices[i][2].resize(
311  dof_handler.tria->n_raw_cells(i) *
312  dof_handler.get_fe().n_dofs_per_quad(
313  0 /*note: in 2D there is only one quad*/),
315 
316  dof_handler.object_dof_ptr[i][2].reserve(
317  dof_handler.tria->n_raw_cells(i) + 1);
318  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1;
319  j++)
320  dof_handler.object_dof_ptr[i][2].push_back(
321  j * dof_handler.get_fe().n_dofs_per_quad(
322  0 /*note: in 2D there is only one quad*/));
323 
324  dof_handler.cell_dof_cache_indices[i].resize(
325  dof_handler.tria->n_raw_cells(i) *
326  dof_handler.get_fe().n_dofs_per_cell(),
328 
329  dof_handler.cell_dof_cache_ptr[i].reserve(
330  dof_handler.tria->n_raw_cells(i) + 1);
331  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1;
332  j++)
333  dof_handler.cell_dof_cache_ptr[i].push_back(
334  j * dof_handler.get_fe().n_dofs_per_cell());
335  }
336 
337  dof_handler.object_dof_indices[0][0].resize(
338  dof_handler.tria->n_vertices() *
339  dof_handler.get_fe().n_dofs_per_vertex(),
341 
342  if (dof_handler.tria->n_cells() > 0)
343  {
344  // line
345  dof_handler.object_dof_ptr[0][1].reserve(
346  dof_handler.tria->n_raw_lines() + 1);
347  for (unsigned int i = 0; i < dof_handler.tria->n_raw_lines() + 1;
348  i++)
349  dof_handler.object_dof_ptr[0][1].push_back(
350  i * dof_handler.get_fe().n_dofs_per_line());
351 
352  dof_handler.object_dof_indices[0][1].resize(
353  dof_handler.tria->n_raw_lines() *
354  dof_handler.get_fe().n_dofs_per_line(),
356  }
357  }
358 
359  template <int spacedim>
360  static void reserve_space(DoFHandler<3, spacedim> &dof_handler)
361  {
362  const unsigned int dim = 3;
363 
364  dof_handler.object_dof_indices[0][0].resize(
365  dof_handler.tria->n_vertices() *
366  dof_handler.get_fe().n_dofs_per_vertex(),
368 
369  for (unsigned int i = 0; i < dof_handler.tria->n_levels(); ++i)
370  {
371  dof_handler.object_dof_indices[i][3].resize(
372  dof_handler.tria->n_raw_cells(i) *
373  dof_handler.get_fe().n_dofs_per_hex(),
375 
376  dof_handler.object_dof_ptr[i][3].reserve(
377  dof_handler.tria->n_raw_cells(i) + 1);
378  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1;
379  j++)
380  dof_handler.object_dof_ptr[i][3].push_back(
381  j * dof_handler.get_fe().n_dofs_per_hex());
382 
383  dof_handler.cell_dof_cache_indices[i].resize(
384  dof_handler.tria->n_raw_cells(i) *
385  dof_handler.get_fe().n_dofs_per_cell(),
387 
388  dof_handler.cell_dof_cache_ptr[i].reserve(
389  dof_handler.tria->n_raw_cells(i) + 1);
390  for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1;
391  j++)
392  dof_handler.cell_dof_cache_ptr[i].push_back(
393  j * dof_handler.get_fe().n_dofs_per_cell());
394  }
395 
396  dof_handler.object_dof_indices[0][0].resize(
397  dof_handler.tria->n_vertices() *
398  dof_handler.get_fe().n_dofs_per_vertex(),
400 
401  if (dof_handler.tria->n_cells() > 0)
402  {
403  // lines
404  dof_handler.object_dof_ptr[0][1].reserve(
405  dof_handler.tria->n_raw_lines() + 1);
406  for (unsigned int i = 0; i < dof_handler.tria->n_raw_lines() + 1;
407  i++)
408  dof_handler.object_dof_ptr[0][1].push_back(
409  i * dof_handler.get_fe().n_dofs_per_line());
410 
411  dof_handler.object_dof_indices[0][1].resize(
412  dof_handler.tria->n_raw_lines() *
413  dof_handler.get_fe().n_dofs_per_line(),
415 
416  // faces
417  {
418  dof_handler.object_dof_ptr[0][2].assign(
419  dof_handler.tria->n_raw_quads() + 1, -1);
420  // determine for each face the number of dofs
421  for (const auto &cell : dof_handler.tria->cell_iterators())
422  for (const auto face_index : cell->face_indices())
423  {
424  const auto &face = cell->face(face_index);
425  const auto n_dofs_per_quad =
426  dof_handler.get_fe().n_dofs_per_quad(face_index);
427 
428  auto &n_dofs_per_quad_target =
429  dof_handler.object_dof_ptr[0][2][face->index() + 1];
430 
431  // make sure that either the face has not been visited or
432  // the face has the same number of dofs assigned
433  Assert(
434  (n_dofs_per_quad_target ==
435  static_cast<
437  -1) ||
438  n_dofs_per_quad_target == n_dofs_per_quad),
440 
441  n_dofs_per_quad_target = n_dofs_per_quad;
442  }
443 
444  // convert the absolute numbers to CRS
445  dof_handler.object_dof_ptr[0][2][0] = 0;
446  for (unsigned int i = 1; i < dof_handler.tria->n_raw_quads() + 1;
447  i++)
448  {
449  if (dof_handler.object_dof_ptr[0][2][i] ==
450  static_cast<
452  dof_handler.object_dof_ptr[0][2][i] =
453  dof_handler.object_dof_ptr[0][2][i - 1];
454  else
455  dof_handler.object_dof_ptr[0][2][i] +=
456  dof_handler.object_dof_ptr[0][2][i - 1];
457  }
458 
459  // allocate memory for indices
460  dof_handler.object_dof_indices[0][2].resize(
461  dof_handler.object_dof_ptr[0][2].back(),
463  }
464  }
465  }
466 
467  template <int spacedim>
468  static void reserve_space_mg(DoFHandler<1, spacedim> &dof_handler)
469  {
470  Assert(dof_handler.get_triangulation().n_levels() > 0,
471  ExcMessage("Invalid triangulation"));
472  dof_handler.clear_mg_space();
473 
474  const ::Triangulation<1, spacedim> &tria =
475  dof_handler.get_triangulation();
476  const unsigned int dofs_per_line =
477  dof_handler.get_fe().n_dofs_per_line();
478  const unsigned int n_levels = tria.n_levels();
479 
480  for (unsigned int i = 0; i < n_levels; ++i)
481  {
482  dof_handler.mg_levels.emplace_back(
484  dof_handler.mg_levels.back()->dof_object.dofs =
485  std::vector<types::global_dof_index>(tria.n_raw_lines(i) *
486  dofs_per_line,
488  }
489 
490  const unsigned int n_vertices = tria.n_vertices();
491 
492  dof_handler.mg_vertex_dofs.resize(n_vertices);
493 
494  std::vector<unsigned int> max_level(n_vertices, 0);
495  std::vector<unsigned int> min_level(n_vertices, n_levels);
496 
497  for (typename ::Triangulation<1, spacedim>::cell_iterator cell =
498  tria.begin();
499  cell != tria.end();
500  ++cell)
501  {
502  const unsigned int level = cell->level();
503 
504  for (const auto vertex : cell->vertex_indices())
505  {
506  const unsigned int vertex_index = cell->vertex_index(vertex);
507 
508  if (min_level[vertex_index] > level)
509  min_level[vertex_index] = level;
510 
511  if (max_level[vertex_index] < level)
512  max_level[vertex_index] = level;
513  }
514  }
515 
516  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
517  if (tria.vertex_used(vertex))
518  {
519  Assert(min_level[vertex] < n_levels, ExcInternalError());
520  Assert(max_level[vertex] >= min_level[vertex],
521  ExcInternalError());
522  dof_handler.mg_vertex_dofs[vertex].init(
523  min_level[vertex],
524  max_level[vertex],
525  dof_handler.get_fe().n_dofs_per_vertex());
526  }
527 
528  else
529  {
530  Assert(min_level[vertex] == n_levels, ExcInternalError());
531  Assert(max_level[vertex] == 0, ExcInternalError());
532  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
533  }
534  }
535 
536  template <int spacedim>
537  static void reserve_space_mg(DoFHandler<2, spacedim> &dof_handler)
538  {
539  Assert(dof_handler.get_triangulation().n_levels() > 0,
540  ExcMessage("Invalid triangulation"));
541  dof_handler.clear_mg_space();
542 
543  const ::FiniteElement<2, spacedim> &fe = dof_handler.get_fe();
544  const ::Triangulation<2, spacedim> &tria =
545  dof_handler.get_triangulation();
546  const unsigned int n_levels = tria.n_levels();
547 
548  for (unsigned int i = 0; i < n_levels; ++i)
549  {
550  dof_handler.mg_levels.emplace_back(
551  std::make_unique<
553  dof_handler.mg_levels.back()->dof_object.dofs =
554  std::vector<types::global_dof_index>(
555  tria.n_raw_quads(i) *
556  fe.n_dofs_per_quad(0 /*note: in 2D there is only one quad*/),
558  }
559 
560  dof_handler.mg_faces =
561  std::make_unique<internal::DoFHandlerImplementation::DoFFaces<2>>();
562  dof_handler.mg_faces->lines.dofs =
563  std::vector<types::global_dof_index>(tria.n_raw_lines() *
564  fe.n_dofs_per_line(),
566 
567  const unsigned int n_vertices = tria.n_vertices();
568 
569  dof_handler.mg_vertex_dofs.resize(n_vertices);
570 
571  std::vector<unsigned int> max_level(n_vertices, 0);
572  std::vector<unsigned int> min_level(n_vertices, n_levels);
573 
574  for (typename ::Triangulation<2, spacedim>::cell_iterator cell =
575  tria.begin();
576  cell != tria.end();
577  ++cell)
578  {
579  const unsigned int level = cell->level();
580 
581  for (const auto vertex : cell->vertex_indices())
582  {
583  const unsigned int vertex_index = cell->vertex_index(vertex);
584 
585  if (min_level[vertex_index] > level)
586  min_level[vertex_index] = level;
587 
588  if (max_level[vertex_index] < level)
589  max_level[vertex_index] = level;
590  }
591  }
592 
593  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
594  if (tria.vertex_used(vertex))
595  {
596  Assert(min_level[vertex] < n_levels, ExcInternalError());
597  Assert(max_level[vertex] >= min_level[vertex],
598  ExcInternalError());
599  dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
600  max_level[vertex],
601  fe.n_dofs_per_vertex());
602  }
603 
604  else
605  {
606  Assert(min_level[vertex] == n_levels, ExcInternalError());
607  Assert(max_level[vertex] == 0, ExcInternalError());
608  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
609  }
610  }
611 
612  template <int spacedim>
613  static void reserve_space_mg(DoFHandler<3, spacedim> &dof_handler)
614  {
615  Assert(dof_handler.get_triangulation().n_levels() > 0,
616  ExcMessage("Invalid triangulation"));
617  dof_handler.clear_mg_space();
618 
619  const ::FiniteElement<3, spacedim> &fe = dof_handler.get_fe();
620  const ::Triangulation<3, spacedim> &tria =
621  dof_handler.get_triangulation();
622  const unsigned int n_levels = tria.n_levels();
623 
624  for (unsigned int i = 0; i < n_levels; ++i)
625  {
626  dof_handler.mg_levels.emplace_back(
627  std::make_unique<
629  dof_handler.mg_levels.back()->dof_object.dofs =
630  std::vector<types::global_dof_index>(tria.n_raw_hexs(i) *
631  fe.n_dofs_per_hex(),
633  }
634 
635  dof_handler.mg_faces =
636  std::make_unique<internal::DoFHandlerImplementation::DoFFaces<3>>();
637  dof_handler.mg_faces->lines.dofs =
638  std::vector<types::global_dof_index>(tria.n_raw_lines() *
639  fe.n_dofs_per_line(),
641 
642  // TODO: the implementation makes the assumption that all faces have the
643  // same number of dofs
644  AssertDimension(fe.n_unique_faces(), 1);
645  dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index>(
646  tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/),
648 
649  const unsigned int n_vertices = tria.n_vertices();
650 
651  dof_handler.mg_vertex_dofs.resize(n_vertices);
652 
653  std::vector<unsigned int> max_level(n_vertices, 0);
654  std::vector<unsigned int> min_level(n_vertices, n_levels);
655 
656  for (typename ::Triangulation<3, spacedim>::cell_iterator cell =
657  tria.begin();
658  cell != tria.end();
659  ++cell)
660  {
661  const unsigned int level = cell->level();
662 
663  for (const auto vertex : cell->vertex_indices())
664  {
665  const unsigned int vertex_index = cell->vertex_index(vertex);
666 
667  if (min_level[vertex_index] > level)
668  min_level[vertex_index] = level;
669 
670  if (max_level[vertex_index] < level)
671  max_level[vertex_index] = level;
672  }
673  }
674 
675  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
676  if (tria.vertex_used(vertex))
677  {
678  Assert(min_level[vertex] < n_levels, ExcInternalError());
679  Assert(max_level[vertex] >= min_level[vertex],
680  ExcInternalError());
681  dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
682  max_level[vertex],
683  fe.n_dofs_per_vertex());
684  }
685 
686  else
687  {
688  Assert(min_level[vertex] == n_levels, ExcInternalError());
689  Assert(max_level[vertex] == 0, ExcInternalError());
690  dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
691  }
692  }
693 
694  template <int spacedim>
697  const DoFHandler<1, spacedim> &dof_handler,
699  &mg_level,
701  &,
702  const unsigned int obj_index,
703  const unsigned int fe_index,
704  const unsigned int local_index,
705  const std::integral_constant<int, 1>)
706  {
707  Assert(dof_handler.hp_capability_enabled == false,
709 
710  return mg_level->dof_object.get_dof_index(
711  static_cast<const DoFHandler<1, spacedim> &>(dof_handler),
712  obj_index,
713  fe_index,
714  local_index);
715  }
716 
717  template <int spacedim>
720  const DoFHandler<2, spacedim> &dof_handler,
722  &,
724  & mg_faces,
725  const unsigned int obj_index,
726  const unsigned int fe_index,
727  const unsigned int local_index,
728  const std::integral_constant<int, 1>)
729  {
730  return mg_faces->lines.get_dof_index(
731  static_cast<const DoFHandler<2, spacedim> &>(dof_handler),
732  obj_index,
733  fe_index,
734  local_index);
735  }
736 
737  template <int spacedim>
740  const DoFHandler<2, spacedim> &dof_handler,
742  &mg_level,
744  &,
745  const unsigned int obj_index,
746  const unsigned int fe_index,
747  const unsigned int local_index,
748  const std::integral_constant<int, 2>)
749  {
750  Assert(dof_handler.hp_capability_enabled == false,
752  return mg_level->dof_object.get_dof_index(
753  static_cast<const DoFHandler<2, spacedim> &>(dof_handler),
754  obj_index,
755  fe_index,
756  local_index);
757  }
758 
759  template <int spacedim>
762  const DoFHandler<3, spacedim> &dof_handler,
764  &,
766  & mg_faces,
767  const unsigned int obj_index,
768  const unsigned int fe_index,
769  const unsigned int local_index,
770  const std::integral_constant<int, 1>)
771  {
772  Assert(dof_handler.hp_capability_enabled == false,
774  return mg_faces->lines.get_dof_index(
775  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
776  obj_index,
777  fe_index,
778  local_index);
779  }
780 
781  template <int spacedim>
784  const DoFHandler<3, spacedim> &dof_handler,
786  &,
788  & mg_faces,
789  const unsigned int obj_index,
790  const unsigned int fe_index,
791  const unsigned int local_index,
792  const std::integral_constant<int, 2>)
793  {
794  Assert(dof_handler.hp_capability_enabled == false,
796  return mg_faces->quads.get_dof_index(
797  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
798  obj_index,
799  fe_index,
800  local_index);
801  }
802 
803  template <int spacedim>
806  const DoFHandler<3, spacedim> &dof_handler,
808  &mg_level,
810  &,
811  const unsigned int obj_index,
812  const unsigned int fe_index,
813  const unsigned int local_index,
814  const std::integral_constant<int, 3>)
815  {
816  Assert(dof_handler.hp_capability_enabled == false,
818  return mg_level->dof_object.get_dof_index(
819  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
820  obj_index,
821  fe_index,
822  local_index);
823  }
824 
825  template <int spacedim>
826  static void
828  const DoFHandler<1, spacedim> &dof_handler,
830  &mg_level,
832  &,
833  const unsigned int obj_index,
834  const unsigned int fe_index,
835  const unsigned int local_index,
837  const std::integral_constant<int, 1>)
838  {
839  Assert(dof_handler.hp_capability_enabled == false,
841  mg_level->dof_object.set_dof_index(
842  static_cast<const DoFHandler<1, spacedim> &>(dof_handler),
843  obj_index,
844  fe_index,
845  local_index,
846  global_index);
847  }
848 
849  template <int spacedim>
850  static void
852  const DoFHandler<2, spacedim> &dof_handler,
854  &,
856  & mg_faces,
857  const unsigned int obj_index,
858  const unsigned int fe_index,
859  const unsigned int local_index,
861  const std::integral_constant<int, 1>)
862  {
863  Assert(dof_handler.hp_capability_enabled == false,
865  mg_faces->lines.set_dof_index(
866  static_cast<const DoFHandler<2, spacedim> &>(dof_handler),
867  obj_index,
868  fe_index,
869  local_index,
870  global_index);
871  }
872 
873  template <int spacedim>
874  static void
876  const DoFHandler<2, spacedim> &dof_handler,
878  &mg_level,
880  &,
881  const unsigned int obj_index,
882  const unsigned int fe_index,
883  const unsigned int local_index,
885  const std::integral_constant<int, 2>)
886  {
887  Assert(dof_handler.hp_capability_enabled == false,
889  mg_level->dof_object.set_dof_index(
890  static_cast<const DoFHandler<2, spacedim> &>(dof_handler),
891  obj_index,
892  fe_index,
893  local_index,
894  global_index);
895  }
896 
897  template <int spacedim>
898  static void
900  const DoFHandler<3, spacedim> &dof_handler,
902  &,
904  & mg_faces,
905  const unsigned int obj_index,
906  const unsigned int fe_index,
907  const unsigned int local_index,
909  const std::integral_constant<int, 1>)
910  {
911  Assert(dof_handler.hp_capability_enabled == false,
913  mg_faces->lines.set_dof_index(
914  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
915  obj_index,
916  fe_index,
917  local_index,
918  global_index);
919  }
920 
921  template <int spacedim>
922  static void
924  const DoFHandler<3, spacedim> &dof_handler,
926  &,
928  & mg_faces,
929  const unsigned int obj_index,
930  const unsigned int fe_index,
931  const unsigned int local_index,
933  const std::integral_constant<int, 2>)
934  {
935  Assert(dof_handler.hp_capability_enabled == false,
937  mg_faces->quads.set_dof_index(
938  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
939  obj_index,
940  fe_index,
941  local_index,
942  global_index);
943  }
944 
945  template <int spacedim>
946  static void
948  const DoFHandler<3, spacedim> &dof_handler,
950  &mg_level,
952  &,
953  const unsigned int obj_index,
954  const unsigned int fe_index,
955  const unsigned int local_index,
957  const std::integral_constant<int, 3>)
958  {
959  Assert(dof_handler.hp_capability_enabled == false,
961  mg_level->dof_object.set_dof_index(
962  static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
963  obj_index,
964  fe_index,
965  local_index,
966  global_index);
967  }
968  };
969  } // namespace DoFHandlerImplementation
970 
971  namespace hp
972  {
973  namespace DoFHandlerImplementation
974  {
980  {
986  template <int dim, int spacedim>
987  static void
989  DoFHandler<dim, spacedim> &dof_handler)
990  {
991  (void)dof_handler;
992  for (const auto &cell : dof_handler.active_cell_iterators())
993  if (cell->is_locally_owned())
994  Assert(
995  !cell->future_fe_index_set(),
996  ExcMessage(
997  "There shouldn't be any cells flagged for p-adaptation when partitioning."));
998  }
999 
1000 
1001 
1006  template <int dim, int spacedim>
1007  static void
1009  {
1010  // Release all space except the fields for active_fe_indices and
1011  // refinement flags which we have to back up before
1012  {
1013  std::vector<std::vector<
1015  active_fe_backup(dof_handler.hp_cell_active_fe_indices.size()),
1016  future_fe_backup(dof_handler.hp_cell_future_fe_indices.size());
1017  for (unsigned int level = 0;
1018  level < dof_handler.hp_cell_future_fe_indices.size();
1019  ++level)
1020  {
1021  active_fe_backup[level] =
1022  std::move(dof_handler.hp_cell_active_fe_indices[level]);
1023  future_fe_backup[level] =
1024  std::move(dof_handler.hp_cell_future_fe_indices[level]);
1025  }
1026 
1027  // delete all levels and set them up newly, since vectors
1028  // are troublesome if you want to change their size
1029  dof_handler.clear_space();
1030 
1031  dof_handler.object_dof_indices.resize(dof_handler.tria->n_levels());
1032  dof_handler.object_dof_ptr.resize(dof_handler.tria->n_levels());
1033  dof_handler.cell_dof_cache_indices.resize(
1034  dof_handler.tria->n_levels());
1035  dof_handler.cell_dof_cache_ptr.resize(dof_handler.tria->n_levels());
1036  dof_handler.hp_cell_active_fe_indices.resize(
1037  dof_handler.tria->n_levels());
1038  dof_handler.hp_cell_future_fe_indices.resize(
1039  dof_handler.tria->n_levels());
1040 
1041  for (unsigned int level = 0; level < dof_handler.tria->n_levels();
1042  ++level)
1043  {
1044  // recover backups
1045  dof_handler.hp_cell_active_fe_indices[level] =
1046  std::move(active_fe_backup[level]);
1047  dof_handler.hp_cell_future_fe_indices[level] =
1048  std::move(future_fe_backup[level]);
1049  }
1050  }
1051  }
1052 
1053 
1054 
1059  template <int dim, int spacedim>
1060  static void
1062  {
1063  // The final step in all of the reserve_space() functions is to set
1064  // up vertex dof information. since vertices are sequentially
1065  // numbered, what we do first is to set up an array in which
1066  // we record whether a vertex is associated with any of the
1067  // given fe's, by setting a bit. in a later step, we then
1068  // actually allocate memory for the required dofs
1069  //
1070  // in the following, we only need to consider vertices that are
1071  // adjacent to either a locally owned or a ghost cell; we never
1072  // store anything on vertices that are only surrounded by
1073  // artificial cells. so figure out that subset of vertices
1074  // first
1075  std::vector<bool> locally_used_vertices(
1076  dof_handler.tria->n_vertices(), false);
1077  for (const auto &cell : dof_handler.active_cell_iterators())
1078  if (!cell->is_artificial())
1079  for (const auto v : cell->vertex_indices())
1080  locally_used_vertices[cell->vertex_index(v)] = true;
1081 
1082  std::vector<std::vector<bool>> vertex_fe_association(
1083  dof_handler.fe_collection.size(),
1084  std::vector<bool>(dof_handler.tria->n_vertices(), false));
1085 
1086  for (const auto &cell : dof_handler.active_cell_iterators())
1087  if (!cell->is_artificial())
1088  for (const auto v : cell->vertex_indices())
1089  vertex_fe_association[cell->active_fe_index()]
1090  [cell->vertex_index(v)] = true;
1091 
1092  // in debug mode, make sure that each vertex is associated
1093  // with at least one fe (note that except for unused
1094  // vertices, all vertices are actually active). this is of
1095  // course only true for vertices that are part of either
1096  // ghost or locally owned cells
1097 #ifdef DEBUG
1098  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
1099  if (locally_used_vertices[v] == true)
1100  if (dof_handler.tria->vertex_used(v) == true)
1101  {
1102  unsigned int fe = 0;
1103  for (; fe < dof_handler.fe_collection.size(); ++fe)
1104  if (vertex_fe_association[fe][v] == true)
1105  break;
1106  Assert(fe != dof_handler.fe_collection.size(),
1107  ExcInternalError());
1108  }
1109 #endif
1110 
1111  const unsigned int d = 0;
1112  const unsigned int l = 0;
1113 
1114  dof_handler.hp_object_fe_ptr[d].clear();
1115  dof_handler.hp_object_fe_indices[d].clear();
1116  dof_handler.object_dof_ptr[l][d].clear();
1117  dof_handler.object_dof_indices[l][d].clear();
1118 
1119  dof_handler.hp_object_fe_ptr[d].reserve(
1120  dof_handler.tria->n_vertices() + 1);
1121 
1122  unsigned int vertex_slots_needed = 0;
1123  unsigned int fe_slots_needed = 0;
1124 
1125  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
1126  {
1127  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1128 
1129  if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
1130  {
1131  for (unsigned int fe = 0;
1132  fe < dof_handler.fe_collection.size();
1133  ++fe)
1134  if (vertex_fe_association[fe][v] == true)
1135  {
1136  fe_slots_needed++;
1137  vertex_slots_needed +=
1138  dof_handler.get_fe(fe).n_dofs_per_vertex();
1139  }
1140  }
1141  }
1142 
1143  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1144 
1145  dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
1146  dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
1147 
1148  dof_handler.object_dof_indices[l][d].reserve(vertex_slots_needed);
1149 
1150  for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
1151  if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
1152  {
1153  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
1154  ++fe)
1155  if (vertex_fe_association[fe][v] == true)
1156  {
1157  dof_handler.hp_object_fe_indices[d].push_back(fe);
1158  dof_handler.object_dof_ptr[l][d].push_back(
1159  dof_handler.object_dof_indices[l][d].size());
1160 
1161  for (unsigned int i = 0;
1162  i < dof_handler.get_fe(fe).n_dofs_per_vertex();
1163  i++)
1164  dof_handler.object_dof_indices[l][d].push_back(
1166  }
1167  }
1168 
1169 
1170  dof_handler.object_dof_ptr[l][d].push_back(
1171  dof_handler.object_dof_indices[l][d].size());
1172 
1173  AssertDimension(vertex_slots_needed,
1174  dof_handler.object_dof_indices[l][d].size());
1175  AssertDimension(fe_slots_needed,
1176  dof_handler.hp_object_fe_indices[d].size());
1177  AssertDimension(fe_slots_needed + 1,
1178  dof_handler.object_dof_ptr[l][d].size());
1179  AssertDimension(dof_handler.tria->n_vertices() + 1,
1180  dof_handler.hp_object_fe_ptr[d].size());
1181 
1182  dof_handler.object_dof_indices[l][d].assign(
1183  vertex_slots_needed, numbers::invalid_dof_index);
1184  }
1185 
1186 
1187 
1192  template <int dim, int spacedim>
1193  static void
1195  {
1196  (void)dof_handler;
1197  // count how much space we need on each level for the cell
1198  // dofs and set the dof_*_offsets data. initially set the
1199  // latter to an invalid index, and only later set it to
1200  // something reasonable for active dof_handler.cells
1201  //
1202  // note that for dof_handler.cells, the situation is simpler
1203  // than for other (lower dimensional) objects since exactly
1204  // one finite element is used for it
1205  for (unsigned int level = 0; level < dof_handler.tria->n_levels();
1206  ++level)
1207  {
1208  dof_handler.object_dof_ptr[level][dim] =
1209  std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
1210  dof_handler.tria->n_raw_cells(level),
1211  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
1212  -1));
1213  dof_handler.cell_dof_cache_ptr[level] =
1214  std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
1215  dof_handler.tria->n_raw_cells(level),
1216  static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
1217  -1));
1218 
1219  types::global_dof_index next_free_dof = 0;
1220  types::global_dof_index cache_size = 0;
1221 
1222  for (auto cell :
1224  if (cell->is_active() && !cell->is_artificial())
1225  {
1226  dof_handler.object_dof_ptr[level][dim][cell->index()] =
1227  next_free_dof;
1228  next_free_dof +=
1229  cell->get_fe().template n_dofs_per_object<dim>();
1230 
1231  dof_handler.cell_dof_cache_ptr[level][cell->index()] =
1232  cache_size;
1233  cache_size += cell->get_fe().n_dofs_per_cell();
1234  }
1235 
1236  dof_handler.object_dof_indices[level][dim] =
1237  std::vector<types::global_dof_index>(
1238  next_free_dof, numbers::invalid_dof_index);
1239  dof_handler.cell_dof_cache_indices[level] =
1240  std::vector<types::global_dof_index>(
1241  cache_size, numbers::invalid_dof_index);
1242  }
1243  }
1244 
1245 
1246 
1251  template <int dim, int spacedim>
1252  static void
1254  {
1255  // FACE DOFS
1256  //
1257  // Count face dofs, then allocate as much space
1258  // as we need and prime the linked list for faces (see the
1259  // description in hp::DoFLevel) with the indices we will
1260  // need. Note that our task is more complicated than for the
1261  // cell case above since two adjacent cells may have different
1262  // active_fe_indices, in which case we need to allocate
1263  // *two* sets of face dofs for the same face. But they don't
1264  // *have* to be different, and so we need to prepare for this
1265  // as well.
1266  //
1267  // The way we do things is that we loop over all active
1268  // cells (these are the only ones that have DoFs
1269  // anyway) and all their faces. We note in the
1270  // user flags whether we have previously visited a face and
1271  // if so skip it (consequently, we have to save and later
1272  // restore the face flags)
1273  {
1274  std::vector<bool> saved_face_user_flags;
1275  switch (dim)
1276  {
1277  case 2:
1278  {
1279  const_cast<::Triangulation<dim, spacedim> &>(
1280  *dof_handler.tria)
1281  .save_user_flags_line(saved_face_user_flags);
1282  const_cast<::Triangulation<dim, spacedim> &>(
1283  *dof_handler.tria)
1284  .clear_user_flags_line();
1285 
1286  break;
1287  }
1288 
1289  case 3:
1290  {
1291  const_cast<::Triangulation<dim, spacedim> &>(
1292  *dof_handler.tria)
1293  .save_user_flags_quad(saved_face_user_flags);
1294  const_cast<::Triangulation<dim, spacedim> &>(
1295  *dof_handler.tria)
1296  .clear_user_flags_quad();
1297 
1298  break;
1299  }
1300 
1301  default:
1302  Assert(false, ExcNotImplemented());
1303  }
1304 
1305  const unsigned int d = dim - 1;
1306  const unsigned int l = 0;
1307 
1308  dof_handler.hp_object_fe_ptr[d].clear();
1309  dof_handler.hp_object_fe_indices[d].clear();
1310  dof_handler.object_dof_ptr[l][d].clear();
1311  dof_handler.object_dof_indices[l][d].clear();
1312 
1313  dof_handler.hp_object_fe_ptr[d].resize(
1314  dof_handler.tria->n_raw_faces() + 1);
1315 
1316  // An array to hold how many slots (see the hp::DoFLevel
1317  // class) we will have to store on each level
1318  unsigned int n_face_slots = 0;
1319 
1320  for (const auto &cell : dof_handler.active_cell_iterators())
1321  if (!cell->is_artificial())
1322  for (const auto face : cell->face_indices())
1323  if (cell->face(face)->user_flag_set() == false)
1324  {
1325  unsigned int fe_slots_needed = 0;
1326 
1327  if (cell->at_boundary(face) ||
1328  cell->face(face)->has_children() ||
1329  cell->neighbor_is_coarser(face) ||
1330  (!cell->at_boundary(face) &&
1331  cell->neighbor(face)->is_artificial()) ||
1332  (!cell->at_boundary(face) &&
1333  !cell->neighbor(face)->is_artificial() &&
1334  (cell->active_fe_index() ==
1335  cell->neighbor(face)->active_fe_index())))
1336  {
1337  fe_slots_needed = 1;
1338  n_face_slots +=
1339  dof_handler.get_fe(cell->active_fe_index())
1340  .template n_dofs_per_object<dim - 1>(face);
1341  }
1342  else
1343  {
1344  fe_slots_needed = 2;
1345  n_face_slots +=
1346  dof_handler.get_fe(cell->active_fe_index())
1347  .template n_dofs_per_object<dim - 1>(face) +
1348  dof_handler
1349  .get_fe(cell->neighbor(face)->active_fe_index())
1350  .template n_dofs_per_object<dim - 1>(
1351  cell->neighbor_face_no(face));
1352  }
1353 
1354  // mark this face as visited
1355  cell->face(face)->set_user_flag();
1356 
1357  dof_handler
1358  .hp_object_fe_ptr[d][cell->face(face)->index() + 1] =
1359  fe_slots_needed;
1360  }
1361 
1362  for (unsigned int i = 1; i < dof_handler.hp_object_fe_ptr[d].size();
1363  i++)
1364  dof_handler.hp_object_fe_ptr[d][i] +=
1365  dof_handler.hp_object_fe_ptr[d][i - 1];
1366 
1367 
1368  dof_handler.hp_object_fe_indices[d].resize(
1369  dof_handler.hp_object_fe_ptr[d].back());
1370  dof_handler.object_dof_ptr[l][d].resize(
1371  dof_handler.hp_object_fe_ptr[d].back() + 1);
1372 
1373  dof_handler.object_dof_indices[l][d].reserve(n_face_slots);
1374 
1375 
1376  // With the memory now allocated, loop over the
1377  // dof_handler cells again and prime the _offset values as
1378  // well as the fe_index fields
1379  switch (dim)
1380  {
1381  case 2:
1382  {
1383  const_cast<::Triangulation<dim, spacedim> &>(
1384  *dof_handler.tria)
1385  .clear_user_flags_line();
1386 
1387  break;
1388  }
1389 
1390  case 3:
1391  {
1392  const_cast<::Triangulation<dim, spacedim> &>(
1393  *dof_handler.tria)
1394  .clear_user_flags_quad();
1395 
1396  break;
1397  }
1398 
1399  default:
1400  Assert(false, ExcNotImplemented());
1401  }
1402 
1403  for (const auto &cell : dof_handler.active_cell_iterators())
1404  if (!cell->is_artificial())
1405  for (const auto face : cell->face_indices())
1406  if (!cell->face(face)->user_flag_set())
1407  {
1408  // Same decision tree as before
1409  if (cell->at_boundary(face) ||
1410  cell->face(face)->has_children() ||
1411  cell->neighbor_is_coarser(face) ||
1412  (!cell->at_boundary(face) &&
1413  cell->neighbor(face)->is_artificial()) ||
1414  (!cell->at_boundary(face) &&
1415  !cell->neighbor(face)->is_artificial() &&
1416  (cell->active_fe_index() ==
1417  cell->neighbor(face)->active_fe_index())))
1418  {
1419  const unsigned int fe = cell->active_fe_index();
1420  const unsigned int n_dofs =
1421  dof_handler.get_fe(fe)
1422  .template n_dofs_per_object<dim - 1>(face);
1423  const unsigned int offset =
1424  dof_handler
1425  .hp_object_fe_ptr[d][cell->face(face)->index()];
1426 
1427  dof_handler.hp_object_fe_indices[d][offset] = fe;
1428  dof_handler.object_dof_ptr[l][d][offset + 1] = n_dofs;
1429 
1430  for (unsigned int i = 0; i < n_dofs; i++)
1431  dof_handler.object_dof_indices[l][d].push_back(
1433  }
1434  else
1435  {
1436  unsigned int fe_1 = cell->active_fe_index();
1437  unsigned int face_no_1 = face;
1438  unsigned int fe_2 =
1439  cell->neighbor(face)->active_fe_index();
1440  unsigned int face_no_2 = cell->neighbor_face_no(face);
1441 
1442  if (fe_2 < fe_1)
1443  {
1444  std::swap(fe_1, fe_2);
1445  std::swap(face_no_1, face_no_2);
1446  }
1447 
1448  const unsigned int n_dofs_1 =
1449  dof_handler.get_fe(fe_1)
1450  .template n_dofs_per_object<dim - 1>(face_no_1);
1451 
1452  const unsigned int n_dofs_2 =
1453  dof_handler.get_fe(fe_2)
1454  .template n_dofs_per_object<dim - 1>(face_no_2);
1455 
1456  const unsigned int offset =
1457  dof_handler
1458  .hp_object_fe_ptr[d][cell->face(face)->index()];
1459 
1460  dof_handler.hp_object_fe_indices[d].push_back(
1461  cell->active_fe_index());
1462  dof_handler.object_dof_ptr[l][d].push_back(
1463  dof_handler.object_dof_indices[l][d].size());
1464 
1465  dof_handler.hp_object_fe_indices[d][offset + 0] =
1466  fe_1;
1467  dof_handler.hp_object_fe_indices[d][offset + 1] =
1468  fe_2;
1469  dof_handler.object_dof_ptr[l][d][offset + 1] =
1470  n_dofs_1;
1471  dof_handler.object_dof_ptr[l][d][offset + 2] =
1472  n_dofs_2;
1473 
1474 
1475  for (unsigned int i = 0; i < n_dofs_1 + n_dofs_2; i++)
1476  dof_handler.object_dof_indices[l][d].push_back(
1478  }
1479 
1480  // mark this face as visited
1481  cell->face(face)->set_user_flag();
1482  }
1483 
1484  for (unsigned int i = 1;
1485  i < dof_handler.object_dof_ptr[l][d].size();
1486  i++)
1487  dof_handler.object_dof_ptr[l][d][i] +=
1488  dof_handler.object_dof_ptr[l][d][i - 1];
1489 
1490  // at the end, restore the user flags for the faces
1491  switch (dim)
1492  {
1493  case 2:
1494  {
1495  const_cast<::Triangulation<dim, spacedim> &>(
1496  *dof_handler.tria)
1497  .load_user_flags_line(saved_face_user_flags);
1498 
1499  break;
1500  }
1501 
1502  case 3:
1503  {
1504  const_cast<::Triangulation<dim, spacedim> &>(
1505  *dof_handler.tria)
1506  .load_user_flags_quad(saved_face_user_flags);
1507 
1508  break;
1509  }
1510 
1511  default:
1512  Assert(false, ExcNotImplemented());
1513  }
1514  }
1515  }
1516 
1517 
1518 
1525  template <int spacedim>
1526  static void reserve_space(::DoFHandler<1, spacedim> &dof_handler)
1527  {
1528  Assert(dof_handler.fe_collection.size() > 0,
1530  Assert(dof_handler.tria->n_levels() > 0,
1531  ExcMessage("The current Triangulation must not be empty."));
1532  Assert(dof_handler.tria->n_levels() ==
1533  dof_handler.hp_cell_future_fe_indices.size(),
1534  ExcInternalError());
1535 
1536  reserve_space_release_space(dof_handler);
1537 
1538  Threads::TaskGroup<> tasks;
1539  tasks +=
1540  Threads::new_task(&reserve_space_cells<1, spacedim>, dof_handler);
1541  tasks += Threads::new_task(&reserve_space_vertices<1, spacedim>,
1542  dof_handler);
1543  tasks.join_all();
1544  }
1545 
1546 
1547 
1548  template <int spacedim>
1549  static void reserve_space(::DoFHandler<2, spacedim> &dof_handler)
1550  {
1551  Assert(dof_handler.fe_collection.size() > 0,
1553  Assert(dof_handler.tria->n_levels() > 0,
1554  ExcMessage("The current Triangulation must not be empty."));
1555  Assert(dof_handler.tria->n_levels() ==
1556  dof_handler.hp_cell_future_fe_indices.size(),
1557  ExcInternalError());
1558 
1559  reserve_space_release_space(dof_handler);
1560 
1561  Threads::TaskGroup<> tasks;
1562  tasks +=
1563  Threads::new_task(&reserve_space_cells<2, spacedim>, dof_handler);
1564  tasks +=
1565  Threads::new_task(&reserve_space_faces<2, spacedim>, dof_handler);
1566  tasks += Threads::new_task(&reserve_space_vertices<2, spacedim>,
1567  dof_handler);
1568  tasks.join_all();
1569  }
1570 
1571 
1572 
1573  template <int spacedim>
1574  static void reserve_space(::DoFHandler<3, spacedim> &dof_handler)
1575  {
1576  Assert(dof_handler.fe_collection.size() > 0,
1578  Assert(dof_handler.tria->n_levels() > 0,
1579  ExcMessage("The current Triangulation must not be empty."));
1580  Assert(dof_handler.tria->n_levels() ==
1581  dof_handler.hp_cell_future_fe_indices.size(),
1582  ExcInternalError());
1583 
1584  reserve_space_release_space(dof_handler);
1585 
1586  Threads::TaskGroup<> tasks;
1587  tasks +=
1588  Threads::new_task(&reserve_space_cells<3, spacedim>, dof_handler);
1589  tasks +=
1590  Threads::new_task(&reserve_space_faces<3, spacedim>, dof_handler);
1591  tasks += Threads::new_task(&reserve_space_vertices<3, spacedim>,
1592  dof_handler);
1593 
1594  // While the tasks above are running, we can turn to line dofs
1595 
1596  // the situation here is pretty much like with vertices:
1597  // there can be an arbitrary number of finite elements
1598  // associated with each line.
1599  //
1600  // the algorithm we use is somewhat similar to what we do in
1601  // reserve_space_vertices()
1602  {
1603  // what we do first is to set up an array in which we
1604  // record whether a line is associated with any of the
1605  // given fe's, by setting a bit. in a later step, we
1606  // then actually allocate memory for the required dofs
1607  std::vector<std::vector<bool>> line_fe_association(
1608  dof_handler.fe_collection.size(),
1609  std::vector<bool>(dof_handler.tria->n_raw_lines(), false));
1610 
1611  for (const auto &cell : dof_handler.active_cell_iterators())
1612  if (!cell->is_artificial())
1613  for (const auto l : cell->line_indices())
1614  line_fe_association[cell->active_fe_index()]
1615  [cell->line_index(l)] = true;
1616 
1617  // first check which of the lines is used at all,
1618  // i.e. is associated with a finite element. we do this
1619  // since not all lines may actually be used, in which
1620  // case we do not have to allocate any memory at all
1621  std::vector<bool> line_is_used(dof_handler.tria->n_raw_lines(),
1622  false);
1623  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1624  ++line)
1625  for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
1626  ++fe)
1627  if (line_fe_association[fe][line] == true)
1628  {
1629  line_is_used[line] = true;
1630  break;
1631  }
1632 
1633 
1634 
1635  const unsigned int d = 1;
1636  const unsigned int l = 0;
1637 
1638  dof_handler.hp_object_fe_ptr[d].clear();
1639  dof_handler.hp_object_fe_indices[d].clear();
1640  dof_handler.object_dof_ptr[l][d].clear();
1641  dof_handler.object_dof_indices[l][d].clear();
1642 
1643  dof_handler.hp_object_fe_ptr[d].reserve(
1644  dof_handler.tria->n_raw_lines() + 1);
1645 
1646  unsigned int line_slots_needed = 0;
1647  unsigned int fe_slots_needed = 0;
1648 
1649  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1650  ++line)
1651  {
1652  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1653 
1654  if (line_is_used[line] == true)
1655  {
1656  for (unsigned int fe = 0;
1657  fe < dof_handler.fe_collection.size();
1658  ++fe)
1659  if (line_fe_association[fe][line] == true)
1660  {
1661  fe_slots_needed++;
1662  line_slots_needed +=
1663  dof_handler.get_fe(fe).n_dofs_per_line();
1664  }
1665  }
1666  }
1667 
1668  dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1669 
1670  // make sure that all entries have been set
1671  AssertDimension(dof_handler.hp_object_fe_ptr[d].size(),
1672  dof_handler.tria->n_raw_lines() + 1);
1673 
1674  dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
1675  dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
1676 
1677  dof_handler.object_dof_indices[l][d].reserve(line_slots_needed);
1678 
1679  for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1680  ++line)
1681  if (line_is_used[line] == true)
1682  {
1683  for (unsigned int fe = 0;
1684  fe < dof_handler.fe_collection.size();
1685  ++fe)
1686  if (line_fe_association[fe][line] == true)
1687  {
1688  dof_handler.hp_object_fe_indices[d].push_back(fe);
1689  dof_handler.object_dof_ptr[l][d].push_back(
1690  dof_handler.object_dof_indices[l][d].size());
1691 
1692  for (unsigned int i = 0;
1693  i < dof_handler.get_fe(fe).n_dofs_per_line();
1694  i++)
1695  dof_handler.object_dof_indices[l][d].push_back(
1697  }
1698  }
1699 
1700  dof_handler.object_dof_ptr[l][d].push_back(
1701  dof_handler.object_dof_indices[l][d].size());
1702 
1703  // make sure that all entries have been set
1704  AssertDimension(dof_handler.hp_object_fe_indices[d].size(),
1705  fe_slots_needed);
1706  AssertDimension(dof_handler.object_dof_ptr[l][d].size(),
1707  fe_slots_needed + 1);
1708  AssertDimension(dof_handler.object_dof_indices[l][d].size(),
1709  line_slots_needed);
1710  }
1711 
1712  // Ensure that everything is done at this point.
1713  tasks.join_all();
1714  }
1715 
1716 
1717 
1729  template <int dim, int spacedim>
1730  static void
1732  {
1733  Assert(
1734  dof_handler.hp_capability_enabled == true,
1736 
1737  if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1738  dynamic_cast<
1739  const ::parallel::shared::Triangulation<dim, spacedim>
1740  *>(&dof_handler.get_triangulation()))
1741  {
1742  // we have a shared triangulation. in this case, every processor
1743  // knows about all cells, but every processor only has knowledge
1744  // about the active_fe_index on the cells it owns.
1745  //
1746  // we can create a complete set of active_fe_indices by letting
1747  // every processor create a vector of indices for all cells,
1748  // filling only those on the cells it owns and setting the indices
1749  // on the other cells to zero. then we add all of these vectors
1750  // up, and because every vector entry has exactly one processor
1751  // that owns it, the sum is correct
1752  std::vector<unsigned int> active_fe_indices(tr->n_active_cells(),
1753  0u);
1754  for (const auto &cell : dof_handler.active_cell_iterators())
1755  if (cell->is_locally_owned())
1756  active_fe_indices[cell->active_cell_index()] =
1757  cell->active_fe_index();
1758 
1759  Utilities::MPI::sum(active_fe_indices,
1760  tr->get_communicator(),
1761  active_fe_indices);
1762 
1763  // now go back and fill the active_fe_index on all other
1764  // cells. we would like to call cell->set_active_fe_index(),
1765  // but that function does not allow setting these indices on
1766  // non-locally_owned cells. so we have to work around the
1767  // issue a little bit by accessing the underlying data
1768  // structures directly
1769  for (const auto &cell : dof_handler.active_cell_iterators())
1770  if (!cell->is_locally_owned())
1771  dof_handler
1772  .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1773  active_fe_indices[cell->active_cell_index()];
1774  }
1775  else if (const ::parallel::
1776  DistributedTriangulationBase<dim, spacedim> *tr =
1777  dynamic_cast<
1778  const ::parallel::
1779  DistributedTriangulationBase<dim, spacedim> *>(
1780  &dof_handler.get_triangulation()))
1781  {
1782  // For completely distributed meshes, use the function that is
1783  // able to move data from locally owned cells on one processor to
1784  // the corresponding ghost cells on others. To this end, we need
1785  // to have functions that can pack and unpack the data we want to
1786  // transport -- namely, the single unsigned int active_fe_index
1787  // objects
1788  auto pack = [](const typename ::DoFHandler<dim, spacedim>::
1789  active_cell_iterator &cell) -> unsigned int {
1790  return cell->active_fe_index();
1791  };
1792 
1793  auto unpack = [&dof_handler](
1794  const typename ::DoFHandler<dim, spacedim>::
1795  active_cell_iterator &cell,
1796  const unsigned int active_fe_index) -> void {
1797  // we would like to say
1798  // cell->set_active_fe_index(active_fe_index);
1799  // but this is not allowed on cells that are not
1800  // locally owned, and we are on a ghost cell
1801  dof_handler
1802  .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1803  active_fe_index;
1804  };
1805 
1807  unsigned int,
1809  static_cast<::DoFHandler<dim, spacedim> &>(dof_handler),
1810  pack,
1811  unpack);
1812  }
1813  else
1814  {
1815  // a sequential triangulation. there is nothing we need to do here
1816  Assert(
1817  (dynamic_cast<
1818  const ::parallel::TriangulationBase<dim, spacedim> *>(
1819  &dof_handler.get_triangulation()) == nullptr),
1820  ExcInternalError());
1821  }
1822  }
1823 
1824 
1825 
1846  template <int dim, int spacedim>
1847  static void
1849  DoFHandler<dim, spacedim> &dof_handler)
1850  {
1851  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1852 
1853  for (const auto &cell : dof_handler.active_cell_iterators())
1854  if (cell->is_locally_owned())
1855  {
1856  if (cell->refine_flag_set())
1857  {
1858  // Store the active_fe_index of each cell that will be
1859  // refined to and distribute it later on its children.
1860  // Pick their future index if flagged for p-refinement.
1861  fe_transfer->refined_cells_fe_index.insert(
1862  {cell, cell->future_fe_index()});
1863  }
1864  else if (cell->coarsen_flag_set())
1865  {
1866  // From all cells that will be coarsened, determine their
1867  // parent and calculate its proper active_fe_index, so that
1868  // it can be set after refinement. But first, check if that
1869  // particular cell has a parent at all.
1870  Assert(cell->level() > 0, ExcInternalError());
1871  const auto &parent = cell->parent();
1872 
1873  // Check if the active_fe_index for the current cell has
1874  // been determined already.
1875  if (fe_transfer->coarsened_cells_fe_index.find(parent) ==
1876  fe_transfer->coarsened_cells_fe_index.end())
1877  {
1878  // Find a suitable active_fe_index for the parent cell
1879  // based on the 'least dominant finite element' of its
1880  // children. Consider the childrens' hypothetical future
1881  // index when they have been flagged for p-refinement.
1882 #ifdef DEBUG
1883  for (const auto &child : parent->child_iterators())
1884  Assert(child->is_active() &&
1885  child->coarsen_flag_set(),
1887  dim>::ExcInconsistentCoarseningFlags());
1888 #endif
1889 
1890  const unsigned int fe_index =
1891  parent->dominated_future_fe_on_children();
1892 
1893  fe_transfer->coarsened_cells_fe_index.insert(
1894  {parent, fe_index});
1895  }
1896  }
1897  else
1898  {
1899  // No h-refinement is scheduled for this cell.
1900  // However, it may have p-refinement indicators, so we
1901  // choose a new active_fe_index based on its flags.
1902  if (cell->future_fe_index_set() == true)
1903  fe_transfer->persisting_cells_fe_index.insert(
1904  {cell, cell->future_fe_index()});
1905  }
1906  }
1907  }
1908 
1909 
1910 
1915  template <int dim, int spacedim>
1916  static void
1918  DoFHandler<dim, spacedim> &dof_handler)
1919  {
1920  const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1921 
1922  // Set active_fe_indices on persisting cells.
1923  for (const auto &persist : fe_transfer->persisting_cells_fe_index)
1924  {
1925  const auto &cell = persist.first;
1926 
1927  if (cell->is_locally_owned())
1928  {
1929  Assert(cell->is_active(), ExcInternalError());
1930  cell->set_active_fe_index(persist.second);
1931  }
1932  }
1933 
1934  // Distribute active_fe_indices from all refined cells on their
1935  // respective children.
1936  for (const auto &refine : fe_transfer->refined_cells_fe_index)
1937  {
1938  const auto &parent = refine.first;
1939 
1940  for (const auto &child : parent->child_iterators())
1941  {
1942  Assert(child->is_locally_owned() && child->is_active(),
1943  ExcInternalError());
1944  child->set_active_fe_index(refine.second);
1945  }
1946  }
1947 
1948  // Set active_fe_indices on coarsened cells that have been determined
1949  // before the actual coarsening happened.
1950  for (const auto &coarsen : fe_transfer->coarsened_cells_fe_index)
1951  {
1952  const auto &cell = coarsen.first;
1953  Assert(cell->is_locally_owned() && cell->is_active(),
1954  ExcInternalError());
1955  cell->set_active_fe_index(coarsen.second);
1956  }
1957  }
1958 
1959 
1970  template <int dim, int spacedim>
1971  static unsigned int
1974  const std::vector<unsigned int> & children_fe_indices,
1975  ::hp::FECollection<dim, spacedim> &fe_collection)
1976  {
1977  Assert(!children_fe_indices.empty(), ExcInternalError());
1978 
1979  // convert vector to set
1980  const std::set<unsigned int> children_fe_indices_set(
1981  children_fe_indices.begin(), children_fe_indices.end());
1982 
1983  const unsigned int dominated_fe_index =
1984  fe_collection.find_dominated_fe_extended(children_fe_indices_set,
1985  /*codim=*/0);
1986 
1987  Assert(dominated_fe_index != numbers::invalid_unsigned_int,
1988  (typename ::DoFCellAccessor<dim, spacedim, false>::
1989  ExcNoDominatedFiniteElementOnChildren()));
1990 
1991  return dominated_fe_index;
1992  }
1993  };
1994  } // namespace DoFHandlerImplementation
1995  } // namespace hp
1996 } // namespace internal
1997 
1998 
1999 
2000 template <int dim, int spacedim>
2001 DoFHandler<dim, spacedim>::DoFHandler(const bool hp_capability_enabled)
2002  : hp_capability_enabled(hp_capability_enabled)
2003  , tria(nullptr, typeid(*this).name())
2004  , mg_faces(nullptr)
2005 {}
2006 
2007 
2008 
2009 template <int dim, int spacedim>
2011  const bool hp_capability_enabled)
2012  : hp_capability_enabled(hp_capability_enabled)
2013  , tria(&tria, typeid(*this).name())
2014  , mg_faces(nullptr)
2015 {
2016  if (hp_capability_enabled)
2017  {
2019  this->create_active_fe_table();
2020  }
2021  else
2022  {
2023  this->setup_policy();
2024  }
2025 }
2026 
2027 template <int dim, int spacedim>
2029 {
2031  {
2032  // unsubscribe as a listener to refinement of the underlying
2033  // triangulation
2034  for (auto &connection : this->tria_listeners)
2035  connection.disconnect();
2036  this->tria_listeners.clear();
2037 
2038  // ...and release allocated memory
2039  // virtual functions called in constructors and destructors never use the
2040  // override in a derived class
2041  // for clarity be explicit on which function is called
2043  }
2044  else
2045  {
2046  // release allocated memory
2047  // virtual functions called in constructors and destructors never use the
2048  // override in a derived class
2049  // for clarity be explicit on which function is called
2051 
2052  // also release the policy. this needs to happen before the
2053  // current object disappears because the policy objects
2054  // store references to the DoFhandler object they work on
2055  this->policy.reset();
2056  }
2057 }
2058 
2059 
2060 
2061 template <int dim, int spacedim>
2062 void
2064  const FiniteElement<dim, spacedim> &fe)
2065 {
2067 }
2068 
2069 
2070 
2071 template <int dim, int spacedim>
2072 void
2075 {
2077  {
2078  this->clear();
2079 
2080  if (this->tria != &tria)
2081  {
2082  for (auto &connection : this->tria_listeners)
2083  connection.disconnect();
2084  this->tria_listeners.clear();
2085 
2086  this->tria = &tria;
2087 
2089  }
2090 
2091  this->create_active_fe_table();
2092 
2093  this->distribute_dofs(fe);
2094  }
2095  else
2096  {
2097  this->tria = &tria;
2098  // this->faces = nullptr;
2099  this->number_cache.n_global_dofs = 0;
2100 
2101  this->setup_policy();
2102 
2103  this->distribute_dofs(fe);
2104  }
2105 }
2106 
2107 
2108 
2109 /*------------------------ Cell iterator functions ------------------------*/
2110 
2111 template <int dim, int spacedim>
2113 DoFHandler<dim, spacedim>::begin(const unsigned int level) const
2114 {
2116  this->get_triangulation().begin(level);
2117  if (cell == this->get_triangulation().end(level))
2118  return end(level);
2119  return cell_iterator(*cell, this);
2120 }
2121 
2122 
2123 
2124 template <int dim, int spacedim>
2127 {
2128  // level is checked in begin
2129  cell_iterator i = begin(level);
2130  if (i.state() != IteratorState::valid)
2131  return i;
2132  while (i->has_children())
2133  if ((++i).state() != IteratorState::valid)
2134  return i;
2135  return i;
2136 }
2137 
2138 
2139 
2140 template <int dim, int spacedim>
2143 {
2144  return cell_iterator(&this->get_triangulation(), -1, -1, this);
2145 }
2146 
2147 
2148 
2149 template <int dim, int spacedim>
2151 DoFHandler<dim, spacedim>::end(const unsigned int level) const
2152 {
2154  this->get_triangulation().end(level);
2155  if (cell.state() != IteratorState::valid)
2156  return end();
2157  return cell_iterator(*cell, this);
2158 }
2159 
2160 
2161 
2162 template <int dim, int spacedim>
2165 {
2167  this->get_triangulation().end_active(level);
2168  if (cell.state() != IteratorState::valid)
2169  return active_cell_iterator(end());
2170  return active_cell_iterator(*cell, this);
2171 }
2172 
2173 
2174 
2175 template <int dim, int spacedim>
2178 {
2179  Assert(this->has_level_dofs(),
2180  ExcMessage("You can only iterate over mg "
2181  "levels if mg dofs got distributed."));
2183  this->get_triangulation().begin(level);
2184  if (cell == this->get_triangulation().end(level))
2185  return end_mg(level);
2186  return level_cell_iterator(*cell, this);
2187 }
2188 
2189 
2190 
2191 template <int dim, int spacedim>
2193 DoFHandler<dim, spacedim>::end_mg(const unsigned int level) const
2194 {
2195  Assert(this->has_level_dofs(),
2196  ExcMessage("You can only iterate over mg "
2197  "levels if mg dofs got distributed."));
2199  this->get_triangulation().end(level);
2200  if (cell.state() != IteratorState::valid)
2201  return end();
2202  return level_cell_iterator(*cell, this);
2203 }
2204 
2205 
2206 
2207 template <int dim, int spacedim>
2210 {
2211  return level_cell_iterator(&this->get_triangulation(), -1, -1, this);
2212 }
2213 
2214 
2215 
2216 template <int dim, int spacedim>
2219 {
2221  begin(), end());
2222 }
2223 
2224 
2225 
2226 template <int dim, int spacedim>
2229 {
2230  return IteratorRange<
2232  end());
2233 }
2234 
2235 
2236 
2237 template <int dim, int spacedim>
2240 {
2242  begin_mg(), end_mg());
2243 }
2244 
2245 
2246 
2247 template <int dim, int spacedim>
2250  const unsigned int level) const
2251 {
2253  begin(level), end(level));
2254 }
2255 
2256 
2257 
2258 template <int dim, int spacedim>
2261  const unsigned int level) const
2262 {
2263  return IteratorRange<
2265  begin_active(level), end_active(level));
2266 }
2267 
2268 
2269 
2270 template <int dim, int spacedim>
2273  const unsigned int level) const
2274 {
2276  begin_mg(level), end_mg(level));
2277 }
2278 
2279 
2280 
2281 //---------------------------------------------------------------------------
2282 
2283 
2284 
2285 template <int dim, int spacedim>
2288 {
2289  Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2291 
2292  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2293 
2294  std::unordered_set<types::global_dof_index> boundary_dofs;
2295  std::vector<types::global_dof_index> dofs_on_face;
2296  dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2297 
2298  const IndexSet &owned_dofs = locally_owned_dofs();
2299 
2300  // loop over all faces to check whether they are at a
2301  // boundary. note that we need not take special care of single
2302  // lines in 3d (using @p{cell->has_boundary_lines}), since we do
2303  // not support boundaries of dimension dim-2, and so every
2304  // boundary line is also part of a boundary face.
2305  for (const auto &cell : this->active_cell_iterators())
2306  if (cell->is_locally_owned() && cell->at_boundary())
2307  {
2308  for (const auto iface : cell->face_indices())
2309  {
2310  const auto face = cell->face(iface);
2311  if (face->at_boundary())
2312  {
2313  const unsigned int dofs_per_face =
2314  cell->get_fe().n_dofs_per_face(iface);
2315  dofs_on_face.resize(dofs_per_face);
2316 
2317  face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2318  for (unsigned int i = 0; i < dofs_per_face; ++i)
2319  {
2320  const unsigned int global_idof_index = dofs_on_face[i];
2321  if (owned_dofs.is_element(global_idof_index))
2322  {
2323  boundary_dofs.insert(global_idof_index);
2324  }
2325  }
2326  }
2327  }
2328  }
2329  return boundary_dofs.size();
2330 }
2331 
2332 
2333 
2334 template <int dim, int spacedim>
2337  const std::set<types::boundary_id> &boundary_ids) const
2338 {
2339  Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2341 
2342  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2343  Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2344  boundary_ids.end(),
2346 
2347  // same as above, but with additional checks for set of boundary
2348  // indicators
2349  std::unordered_set<types::global_dof_index> boundary_dofs;
2350  std::vector<types::global_dof_index> dofs_on_face;
2351  dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2352 
2353  const IndexSet &owned_dofs = locally_owned_dofs();
2354 
2355  for (const auto &cell : this->active_cell_iterators())
2356  if (cell->is_locally_owned() && cell->at_boundary())
2357  {
2358  for (const auto iface : cell->face_indices())
2359  {
2360  const auto face = cell->face(iface);
2361  const unsigned int boundary_id = face->boundary_id();
2362  if (face->at_boundary() &&
2363  (boundary_ids.find(boundary_id) != boundary_ids.end()))
2364  {
2365  const unsigned int dofs_per_face =
2366  cell->get_fe().n_dofs_per_face(iface);
2367  dofs_on_face.resize(dofs_per_face);
2368 
2369  face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2370  for (unsigned int i = 0; i < dofs_per_face; ++i)
2371  {
2372  const unsigned int global_idof_index = dofs_on_face[i];
2373  if (owned_dofs.is_element(global_idof_index))
2374  {
2375  boundary_dofs.insert(global_idof_index);
2376  }
2377  }
2378  }
2379  }
2380  }
2381  return boundary_dofs.size();
2382 }
2383 
2384 
2385 
2386 template <int dim, int spacedim>
2387 std::size_t
2389 {
2390  std::size_t mem = MemoryConsumption::memory_consumption(this->tria) +
2393 
2402 
2403 
2405  {
2406  // nothing to add
2407  }
2408  else
2409  {
2410  // collect size of multigrid data structures
2411 
2413 
2414  for (unsigned int level = 0; level < this->mg_levels.size(); ++level)
2415  mem += this->mg_levels[level]->memory_consumption();
2416 
2417  if (this->mg_faces != nullptr)
2419 
2420  for (unsigned int i = 0; i < this->mg_vertex_dofs.size(); ++i)
2421  mem += sizeof(MGVertexDoFs) +
2422  (1 + this->mg_vertex_dofs[i].get_finest_level() -
2423  this->mg_vertex_dofs[i].get_coarsest_level()) *
2424  sizeof(types::global_dof_index);
2425  }
2426 
2427  return mem;
2428 }
2429 
2430 
2431 
2432 template <int dim, int spacedim>
2433 void
2435 {
2437 }
2438 
2439 
2440 
2441 template <int dim, int spacedim>
2442 void
2444 {
2445  Assert(
2446  this->tria != nullptr,
2447  ExcMessage(
2448  "You need to set the Triangulation in the DoFHandler using initialize() or "
2449  "in the constructor before you can distribute DoFs."));
2450  Assert(this->tria->n_levels() > 0,
2451  ExcMessage("The Triangulation you are using is empty!"));
2452  Assert(ff.size() > 0, ExcMessage("The hp::FECollection given is empty!"));
2453 
2454  // don't create a new object if the one we have is already appropriate
2455  if (this->fe_collection != ff)
2457 
2459  {
2460  // ensure that the active_fe_indices vectors are initialized correctly
2461  this->create_active_fe_table();
2462 
2463  // make sure every processor knows the active_fe_indices
2464  // on both its own cells and all ghost cells
2467 
2468  // make sure that the fe collection is large enough to
2469  // cover all fe indices presently in use on the mesh
2470  for (const auto &cell : this->active_cell_iterators())
2471  if (!cell->is_artificial())
2472  Assert(cell->active_fe_index() < this->fe_collection.size(),
2473  ExcInvalidFEIndex(cell->active_fe_index(),
2474  this->fe_collection.size()));
2475  }
2476 }
2477 
2478 
2479 
2480 template <int dim, int spacedim>
2481 void
2483  const FiniteElement<dim, spacedim> &fe)
2484 {
2486 }
2487 
2488 
2489 
2490 template <int dim, int spacedim>
2491 void
2494 {
2496  {
2497  object_dof_indices.resize(this->tria->n_levels());
2498  object_dof_ptr.resize(this->tria->n_levels());
2499  cell_dof_cache_indices.resize(this->tria->n_levels());
2500  cell_dof_cache_ptr.resize(this->tria->n_levels());
2501  hp_cell_active_fe_indices.resize(this->tria->n_levels());
2502  hp_cell_future_fe_indices.resize(this->tria->n_levels());
2503  // assign the fe_collection and initialize all active_fe_indices
2504  this->set_fe(ff);
2505 
2506  // If an underlying shared::Tria allows artificial cells,
2507  // then save the current set of subdomain ids, and set
2508  // subdomain ids to the "true" owner of each cell. we later
2509  // restore these flags
2510  std::vector<types::subdomain_id> saved_subdomain_ids;
2512  (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
2513  &this->get_triangulation()));
2514  if (shared_tria != nullptr && shared_tria->with_artificial_cells())
2515  {
2516  saved_subdomain_ids.resize(shared_tria->n_active_cells());
2517 
2518  const std::vector<types::subdomain_id> &true_subdomain_ids =
2519  shared_tria->get_true_subdomain_ids_of_cells();
2520 
2521  for (const auto &cell : shared_tria->active_cell_iterators())
2522  {
2523  const unsigned int index = cell->active_cell_index();
2524  saved_subdomain_ids[index] = cell->subdomain_id();
2525  cell->set_subdomain_id(true_subdomain_ids[index]);
2526  }
2527  }
2528 
2529  // then allocate space for all the other tables
2531  reserve_space(*this);
2532 
2533  // now undo the subdomain modification
2534  if (shared_tria != nullptr && shared_tria->with_artificial_cells())
2535  for (const auto &cell : shared_tria->active_cell_iterators())
2536  cell->set_subdomain_id(
2537  saved_subdomain_ids[cell->active_cell_index()]);
2538 
2539 
2540  // Clear user flags because we will need them. But first we save
2541  // them and make sure that we restore them later such that at the
2542  // end of this function the Triangulation will be in the same
2543  // state as it was at the beginning of this function.
2544  std::vector<bool> user_flags;
2545  this->tria->save_user_flags(user_flags);
2546  const_cast<Triangulation<dim, spacedim> &>(*this->tria)
2547  .clear_user_flags();
2548 
2549 
2551 
2552  // Now for the real work:
2553  this->number_cache = this->policy->distribute_dofs();
2554 
2556 
2557  // do some housekeeping: compress indices
2558  //{
2559  // Threads::TaskGroup<> tg;
2560  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2561  // tg += Threads::new_task(
2562  // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2563  // *this->levels_hp[level],
2564  // this->fe_collection);
2565  // tg.join_all();
2566  //}
2567 
2568  // finally restore the user flags
2569  const_cast<Triangulation<dim, spacedim> &>(*this->tria)
2570  .load_user_flags(user_flags);
2571  }
2572  else
2573  {
2574  // first, assign the finite_element
2575  this->set_fe(ff);
2576 
2577  // delete all levels and set them up newly. note that we still have to
2578  // allocate space for all degrees of freedom on this mesh (including ghost
2579  // and cells that are entirely stored on different processors), though we
2580  // may not assign numbers to some of them (i.e. they will remain at
2581  // invalid_dof_index). We need to allocate the space because we will want
2582  // to be able to query the dof_indices on each cell, and simply be told
2583  // that we don't know them on some cell (i.e. get back invalid_dof_index)
2584  this->clear_space();
2585  object_dof_indices.resize(this->tria->n_levels());
2586  object_dof_ptr.resize(this->tria->n_levels());
2587  cell_dof_cache_indices.resize(this->tria->n_levels());
2588  cell_dof_cache_ptr.resize(this->tria->n_levels());
2590 
2591  // hand things off to the policy
2592  this->number_cache = this->policy->distribute_dofs();
2593 
2594  // initialize the block info object only if this is a sequential
2595  // triangulation. it doesn't work correctly yet if it is parallel
2596  if (dynamic_cast<
2598  &*this->tria) == nullptr)
2599  this->block_info_object.initialize(*this, false, true);
2600  }
2601 }
2602 
2603 
2604 
2605 template <int dim, int spacedim>
2606 void
2608 {
2610 
2611  Assert(
2612  this->object_dof_indices.size() > 0,
2613  ExcMessage(
2614  "Distribute active DoFs using distribute_dofs() before calling distribute_mg_dofs()."));
2615 
2616  Assert(
2617  ((this->tria->get_mesh_smoothing() &
2620  ExcMessage(
2621  "The mesh smoothing requirement 'limit_level_difference_at_vertices' has to be set for using multigrid!"));
2622 
2623  this->clear_mg_space();
2624 
2626  this->mg_number_cache = this->policy->distribute_mg_dofs();
2627 
2628  // initialize the block info object only if this is a sequential
2629  // triangulation. it doesn't work correctly yet if it is parallel
2630  if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2631  &*this->tria) == nullptr)
2632  this->block_info_object.initialize(*this, true, false);
2633 }
2634 
2635 
2636 
2637 template <int dim, int spacedim>
2638 void
2640 {
2642 
2643  this->block_info_object.initialize_local(*this);
2644 }
2645 
2646 
2647 
2648 template <int dim, int spacedim>
2649 void
2651 {
2652  // decide whether we need a sequential or a parallel distributed policy
2653  if (dynamic_cast<const ::parallel::shared::Triangulation<dim, spacedim>
2654  *>(&this->get_triangulation()) != nullptr)
2655  this->policy = std::make_unique<internal::DoFHandlerImplementation::Policy::
2656  ParallelShared<dim, spacedim>>(*this);
2657  else if (dynamic_cast<
2658  const ::parallel::DistributedTriangulationBase<dim, spacedim>
2659  *>(&this->get_triangulation()) == nullptr)
2660  this->policy = std::make_unique<
2662  *this);
2663  else
2664  this->policy =
2665  std::make_unique<internal::DoFHandlerImplementation::Policy::
2666  ParallelDistributed<dim, spacedim>>(*this);
2667 }
2668 
2669 
2670 
2671 template <int dim, int spacedim>
2672 void
2674 {
2676  {
2677  // release memory
2678  this->clear_space();
2679  }
2680  else
2681  {
2682  // release memory
2683  this->clear_space();
2684  this->clear_mg_space();
2685  }
2686 }
2687 
2688 
2689 
2690 template <int dim, int spacedim>
2691 void
2693 {
2694  cell_dof_cache_indices.clear();
2695 
2696  cell_dof_cache_ptr.clear();
2697 
2698  object_dof_indices.clear();
2699 
2700  object_dof_ptr.clear();
2701 
2703  {
2704  this->hp_cell_active_fe_indices.clear();
2705  this->hp_cell_future_fe_indices.clear();
2706 
2707  object_dof_indices.clear();
2708  }
2709  else
2710  {
2711  this->number_cache.clear();
2712  }
2713 }
2714 
2715 
2716 
2717 template <int dim, int spacedim>
2718 void
2720 {
2721  this->mg_levels.clear();
2722  this->mg_faces.reset();
2723 
2724  std::vector<MGVertexDoFs> tmp;
2725 
2726  std::swap(this->mg_vertex_dofs, tmp);
2727 
2728  this->mg_number_cache.clear();
2729 }
2730 
2731 
2732 
2733 template <int dim, int spacedim>
2734 void
2736  const std::vector<types::global_dof_index> &new_numbers)
2737 {
2739  {
2740  Assert(this->hp_cell_future_fe_indices.size() > 0,
2741  ExcMessage(
2742  "You need to distribute DoFs before you can renumber them."));
2743 
2744  AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2745 
2746 #ifdef DEBUG
2747  // assert that the new indices are consecutively numbered if we are
2748  // working on a single processor. this doesn't need to
2749  // hold in the case of a parallel mesh since we map the interval
2750  // [0...n_dofs()) into itself but only globally, not on each processor
2751  if (this->n_locally_owned_dofs() == this->n_dofs())
2752  {
2753  std::vector<types::global_dof_index> tmp(new_numbers);
2754  std::sort(tmp.begin(), tmp.end());
2755  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2757  for (; p != tmp.end(); ++p, ++i)
2758  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2759  }
2760  else
2761  for (const auto new_number : new_numbers)
2762  Assert(new_number < this->n_dofs(),
2763  ExcMessage(
2764  "New DoF index is not less than the total number of dofs."));
2765 #endif
2766 
2767  // uncompress the internal storage scheme of dofs on cells so that
2768  // we can access dofs in turns. uncompress in parallel, starting
2769  // with the most expensive levels (the highest ones)
2770  //{
2771  // Threads::TaskGroup<> tg;
2772  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2773  // tg += Threads::new_task(
2774  // &::internal::hp::DoFLevel::uncompress_data<dim, spacedim>,
2775  // *this->levels_hp[level],
2776  // this->fe_collection);
2777  // tg.join_all();
2778  //}
2779 
2780  // do the renumbering
2781  this->number_cache = this->policy->renumber_dofs(new_numbers);
2782 
2783  // now re-compress the dof indices
2784  //{
2785  // Threads::TaskGroup<> tg;
2786  // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2787  // tg += Threads::new_task(
2788  // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2789  // *this->levels_hp[level],
2790  // this->fe_collection);
2791  // tg.join_all();
2792  //}
2793  }
2794  else
2795  {
2796  Assert(this->object_dof_indices.size() > 0,
2797  ExcMessage(
2798  "You need to distribute DoFs before you can renumber them."));
2799 
2800 #ifdef DEBUG
2801  if (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
2802  &*this->tria) != nullptr)
2803  {
2804  Assert(new_numbers.size() == this->n_dofs() ||
2805  new_numbers.size() == this->n_locally_owned_dofs(),
2806  ExcMessage("Incorrect size of the input array."));
2807  }
2808  else if (dynamic_cast<
2810  &*this->tria) != nullptr)
2811  {
2812  AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2813  }
2814  else
2815  {
2816  AssertDimension(new_numbers.size(), this->n_dofs());
2817  }
2818 
2819  // assert that the new indices are consecutively numbered if we are
2820  // working on a single processor. this doesn't need to
2821  // hold in the case of a parallel mesh since we map the interval
2822  // [0...n_dofs()) into itself but only globally, not on each processor
2823  if (this->n_locally_owned_dofs() == this->n_dofs())
2824  {
2825  std::vector<types::global_dof_index> tmp(new_numbers);
2826  std::sort(tmp.begin(), tmp.end());
2827  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2829  for (; p != tmp.end(); ++p, ++i)
2830  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2831  }
2832  else
2833  for (const auto new_number : new_numbers)
2834  Assert(new_number < this->n_dofs(),
2835  ExcMessage(
2836  "New DoF index is not less than the total number of dofs."));
2837 #endif
2838 
2839  this->number_cache = this->policy->renumber_dofs(new_numbers);
2840  }
2841 }
2842 
2843 
2844 
2845 template <int dim, int spacedim>
2846 void
2848  const unsigned int level,
2849  const std::vector<types::global_dof_index> &new_numbers)
2850 {
2852 
2853  Assert(
2854  this->mg_levels.size() > 0 && this->object_dof_indices.size() > 0,
2855  ExcMessage(
2856  "You need to distribute active and level DoFs before you can renumber level DoFs."));
2857  AssertIndexRange(level, this->get_triangulation().n_global_levels());
2858  AssertDimension(new_numbers.size(),
2859  this->locally_owned_mg_dofs(level).n_elements());
2860 
2861 #ifdef DEBUG
2862  // assert that the new indices are consecutively numbered if we are working
2863  // on a single processor. this doesn't need to hold in the case of a
2864  // parallel mesh since we map the interval [0...n_dofs(level)) into itself
2865  // but only globally, not on each processor
2866  if (this->n_locally_owned_dofs() == this->n_dofs())
2867  {
2868  std::vector<types::global_dof_index> tmp(new_numbers);
2869  std::sort(tmp.begin(), tmp.end());
2870  std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2872  for (; p != tmp.end(); ++p, ++i)
2873  Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2874  }
2875  else
2876  for (const auto new_number : new_numbers)
2877  Assert(new_number < this->n_dofs(level),
2878  ExcMessage(
2879  "New DoF index is not less than the total number of dofs."));
2880 #endif
2881 
2882  this->mg_number_cache[level] =
2883  this->policy->renumber_mg_dofs(level, new_numbers);
2884 }
2885 
2886 
2887 
2888 template <int dim, int spacedim>
2889 unsigned int
2891 {
2892  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2893 
2894  switch (dim)
2895  {
2896  case 1:
2897  return this->fe_collection.max_dofs_per_vertex();
2898  case 2:
2899  return (3 * this->fe_collection.max_dofs_per_vertex() +
2900  2 * this->fe_collection.max_dofs_per_line());
2901  case 3:
2902  // we need to take refinement of one boundary face into
2903  // consideration here; in fact, this function returns what
2904  // #max_coupling_between_dofs<2> returns
2905  //
2906  // we assume here, that only four faces meet at the boundary;
2907  // this assumption is not justified and needs to be fixed some
2908  // time. fortunately, omitting it for now does no harm since
2909  // the matrix will cry foul if its requirements are not
2910  // satisfied
2911  return (19 * this->fe_collection.max_dofs_per_vertex() +
2912  28 * this->fe_collection.max_dofs_per_line() +
2913  8 * this->fe_collection.max_dofs_per_quad());
2914  default:
2915  Assert(false, ExcNotImplemented());
2916  return 0;
2917  }
2918 }
2919 
2920 
2921 
2922 template <int dim, int spacedim>
2923 unsigned int
2925 {
2926  Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2929 }
2930 
2931 
2932 
2933 template <int dim, int spacedim>
2934 template <int structdim>
2936 DoFHandler<dim, spacedim>::get_dof_index(const unsigned int obj_level,
2937  const unsigned int obj_index,
2938  const unsigned int fe_index,
2939  const unsigned int local_index) const
2940 {
2942  {
2943  Assert(false, ExcNotImplemented());
2945  }
2946  else
2947  {
2949  *this,
2950  this->mg_levels[obj_level],
2951  this->mg_faces,
2952  obj_index,
2953  fe_index,
2954  local_index,
2955  std::integral_constant<int, structdim>());
2956  }
2957 }
2958 
2959 
2960 
2961 template <int dim, int spacedim>
2962 template <int structdim>
2963 void
2965  const unsigned int obj_level,
2966  const unsigned int obj_index,
2967  const unsigned int fe_index,
2968  const unsigned int local_index,
2970 {
2972  {
2973  Assert(false, ExcNotImplemented());
2974  return;
2975  }
2976  else
2977  {
2979  *this,
2980  this->mg_levels[obj_level],
2981  this->mg_faces,
2982  obj_index,
2983  fe_index,
2984  local_index,
2985  global_index,
2986  std::integral_constant<int, structdim>());
2987  }
2988 }
2989 
2990 
2991 
2992 template <int dim, int spacedim>
2993 void
2995  const std::vector<unsigned int> &active_fe_indices)
2996 {
2997  Assert(active_fe_indices.size() == this->get_triangulation().n_active_cells(),
2998  ExcDimensionMismatch(active_fe_indices.size(),
2999  this->get_triangulation().n_active_cells()));
3000 
3001  this->create_active_fe_table();
3002  // we could set the values directly, since they are stored as
3003  // protected data of this object, but for simplicity we use the
3004  // cell-wise access. this way we also have to pass some debug-mode
3005  // tests which we would have to duplicate ourselves otherwise
3006  for (const auto &cell : this->active_cell_iterators())
3007  if (cell->is_locally_owned())
3008  cell->set_active_fe_index(active_fe_indices[cell->active_cell_index()]);
3009 }
3010 
3011 
3012 
3013 template <int dim, int spacedim>
3014 void
3016  std::vector<unsigned int> &active_fe_indices) const
3017 {
3018  active_fe_indices.resize(this->get_triangulation().n_active_cells());
3019 
3020  // we could try to extract the values directly, since they are
3021  // stored as protected data of this object, but for simplicity we
3022  // use the cell-wise access.
3023  for (const auto &cell : this->active_cell_iterators())
3024  if (!cell->is_artificial())
3025  active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
3026 }
3027 
3028 
3029 
3030 template <int dim, int spacedim>
3031 void
3033 {
3034  // connect functions to signals of the underlying triangulation
3035  this->tria_listeners.push_back(this->tria->signals.pre_refinement.connect(
3036  [this]() { this->pre_refinement_action(); }));
3037  this->tria_listeners.push_back(this->tria->signals.post_refinement.connect(
3038  [this]() { this->post_refinement_action(); }));
3039  this->tria_listeners.push_back(this->tria->signals.create.connect(
3040  [this]() { this->post_refinement_action(); }));
3041 
3042  // decide whether we need a sequential or a parallel shared/distributed
3043  // policy and attach corresponding callback functions dealing with the
3044  // transfer of active_fe_indices
3045  if (dynamic_cast<
3046  const ::parallel::DistributedTriangulationBase<dim, spacedim> *>(
3047  &this->get_triangulation()))
3048  {
3049  this->policy =
3050  std::make_unique<internal::DoFHandlerImplementation::Policy::
3051  ParallelDistributed<dim, spacedim>>(*this);
3052 
3053  // repartitioning signals
3054  this->tria_listeners.push_back(
3055  this->tria->signals.pre_distributed_repartition.connect([this]() {
3056  internal::hp::DoFHandlerImplementation::Implementation::
3057  ensure_absence_of_future_fe_indices<dim, spacedim>(*this);
3058  }));
3059  this->tria_listeners.push_back(
3060  this->tria->signals.pre_distributed_repartition.connect(
3061  [this]() { this->pre_distributed_active_fe_index_transfer(); }));
3062  this->tria_listeners.push_back(
3063  this->tria->signals.post_distributed_repartition.connect(
3064  [this] { this->post_distributed_active_fe_index_transfer(); }));
3065 
3066  // refinement signals
3067  this->tria_listeners.push_back(
3068  this->tria->signals.pre_distributed_refinement.connect(
3069  [this]() { this->pre_distributed_active_fe_index_transfer(); }));
3070  this->tria_listeners.push_back(
3071  this->tria->signals.post_distributed_refinement.connect(
3072  [this]() { this->post_distributed_active_fe_index_transfer(); }));
3073 
3074  // serialization signals
3075  this->tria_listeners.push_back(
3076  this->tria->signals.post_distributed_save.connect([this]() {
3077  this->post_distributed_serialization_of_active_fe_indices();
3078  }));
3079  }
3080  else if (dynamic_cast<
3081  const ::parallel::shared::Triangulation<dim, spacedim> *>(
3082  &this->get_triangulation()) != nullptr)
3083  {
3084  this->policy =
3085  std::make_unique<internal::DoFHandlerImplementation::Policy::
3086  ParallelShared<dim, spacedim>>(*this);
3087 
3088  // partitioning signals
3089  this->tria_listeners.push_back(
3090  this->tria->signals.pre_partition.connect([this]() {
3091  internal::hp::DoFHandlerImplementation::Implementation::
3092  ensure_absence_of_future_fe_indices(*this);
3093  }));
3094 
3095  // refinement signals
3096  this->tria_listeners.push_back(this->tria->signals.pre_refinement.connect(
3097  [this] { this->pre_active_fe_index_transfer(); }));
3098  this->tria_listeners.push_back(
3099  this->tria->signals.post_refinement.connect(
3100  [this] { this->post_active_fe_index_transfer(); }));
3101  }
3102  else
3103  {
3104  this->policy = std::make_unique<
3106  *this);
3107 
3108  // refinement signals
3109  this->tria_listeners.push_back(this->tria->signals.pre_refinement.connect(
3110  [this] { this->pre_active_fe_index_transfer(); }));
3111  this->tria_listeners.push_back(
3112  this->tria->signals.post_refinement.connect(
3113  [this] { this->post_active_fe_index_transfer(); }));
3114  }
3115 }
3116 
3117 
3118 
3119 template <int dim, int spacedim>
3120 void
3122 {
3124 
3125 
3126  // Create sufficiently many hp::DoFLevels.
3127  // while (this->levels_hp.size() < this->tria->n_levels())
3128  // this->levels_hp.emplace_back(new ::internal::hp::DoFLevel);
3129 
3130  while (this->hp_cell_active_fe_indices.size() < this->tria->n_levels())
3131  this->hp_cell_active_fe_indices.push_back({});
3132 
3133  while (this->hp_cell_future_fe_indices.size() < this->tria->n_levels())
3134  this->hp_cell_future_fe_indices.push_back({});
3135 
3136  // then make sure that on each level we have the appropriate size
3137  // of active_fe_indices; preset them to zero, i.e. the default FE
3138  for (unsigned int level = 0; level < this->hp_cell_future_fe_indices.size();
3139  ++level)
3140  {
3141  if (this->hp_cell_active_fe_indices[level].size() == 0 &&
3142  this->hp_cell_future_fe_indices[level].size() == 0)
3143  {
3144  this->hp_cell_active_fe_indices[level].resize(
3145  this->tria->n_raw_cells(level), 0);
3146  this->hp_cell_future_fe_indices[level].resize(
3147  this->tria->n_raw_cells(level), invalid_active_fe_index);
3148  }
3149  else
3150  {
3151  // Either the active_fe_indices have size zero because
3152  // they were just created, or the correct size. Other
3153  // sizes indicate that something went wrong.
3154  Assert(this->hp_cell_active_fe_indices[level].size() ==
3155  this->tria->n_raw_cells(level) &&
3156  this->hp_cell_future_fe_indices[level].size() ==
3157  this->tria->n_raw_cells(level),
3158  ExcInternalError());
3159  }
3160 
3161  // it may be that the previous table was compressed; in that
3162  // case, restore the correct active_fe_index. the fact that
3163  // this no longer matches the indices in the table is of no
3164  // importance because the current function is called at a
3165  // point where we have to recreate the dof_indices tables in
3166  // the levels anyway
3167  // this->levels_hp[level]->normalize_active_fe_indices();
3168  }
3169 }
3170 
3171 
3172 
3173 template <int dim, int spacedim>
3174 void
3176 {
3178 }
3179 
3180 
3181 
3182 template <int dim, int spacedim>
3183 void
3185 {
3186  // // Normally only one level is added, but if this Triangulation
3187  // // is created by copy_triangulation, it can be more than one level.
3188  // while (this->levels_hp.size() < this->tria->n_levels())
3189  // this->levels_hp.emplace_back(new ::internal::hp::DoFLevel);
3190  //
3191  // // Coarsening can lead to the loss of levels. Hence remove them.
3192  // while (this->levels_hp.size() > this->tria->n_levels())
3193  // {
3194  // // drop the last element. that also releases the memory pointed to
3195  // this->levels_hp.pop_back();
3196  // }
3197 
3198  while (this->hp_cell_active_fe_indices.size() < this->tria->n_levels())
3199  this->hp_cell_active_fe_indices.push_back({});
3200 
3201  while (this->hp_cell_active_fe_indices.size() > this->tria->n_levels())
3202  this->hp_cell_active_fe_indices.pop_back();
3203 
3204  while (this->hp_cell_future_fe_indices.size() < this->tria->n_levels())
3205  this->hp_cell_future_fe_indices.push_back({});
3206 
3207  while (this->hp_cell_future_fe_indices.size() > this->tria->n_levels())
3208  this->hp_cell_future_fe_indices.pop_back();
3209 
3210 
3211 
3212  Assert(this->hp_cell_future_fe_indices.size() == this->tria->n_levels(),
3213  ExcInternalError());
3214  for (unsigned int i = 0; i < this->hp_cell_future_fe_indices.size(); ++i)
3215  {
3216  // Resize active_fe_indices vectors. Use zero indicator to extend.
3217  this->hp_cell_active_fe_indices[i].resize(this->tria->n_raw_cells(i), 0);
3218 
3219  // Resize future_fe_indices vectors. Make sure that all
3220  // future_fe_indices have been cleared after refinement happened.
3221  //
3222  // We have used future_fe_indices to update all active_fe_indices
3223  // before refinement happened, thus we are safe to clear them now.
3224  this->hp_cell_future_fe_indices[i].assign(this->tria->n_raw_cells(i),
3226  }
3227 }
3228 
3229 
3230 template <int dim, int spacedim>
3231 void
3233 {
3234  // Finite elements need to be assigned to each cell by calling
3235  // distribute_dofs() first to make this functionality available.
3236  if (this->fe_collection.size() > 0)
3237  {
3238  Assert(this->active_fe_index_transfer == nullptr, ExcInternalError());
3239 
3240  this->active_fe_index_transfer =
3241  std::make_unique<ActiveFEIndexTransfer>();
3242 
3245  }
3246 }
3247 
3248 
3249 
3250 template <int dim, int spacedim>
3251 void
3253 {
3254 #ifndef DEAL_II_WITH_P4EST
3255  Assert(false,
3256  ExcMessage(
3257  "You are attempting to use a functionality that is only available "
3258  "if deal.II was configured to use p4est, but cmake did not find a "
3259  "valid p4est library."));
3260 #else
3261  // the implementation below requires a p:d:T currently
3262  Assert(
3264  &this->get_triangulation()) != nullptr),
3265  ExcNotImplemented());
3266 
3267  // Finite elements need to be assigned to each cell by calling
3268  // distribute_dofs() first to make this functionality available.
3269  if (fe_collection.size() > 0)
3270  {
3272 
3273  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3274 
3275  // If we work on a p::d::Triangulation, we have to transfer all
3276  // active_fe_indices since ownership of cells may change. We will
3277  // use our p::d::CellDataTransfer member to achieve this. Further,
3278  // we prepare the values in such a way that they will correspond to
3279  // the active_fe_indices on the new mesh.
3280 
3281  // Gather all current future_fe_indices.
3282  active_fe_index_transfer->active_fe_indices.resize(
3284 
3285  for (const auto &cell : active_cell_iterators())
3286  if (cell->is_locally_owned())
3288  ->active_fe_indices[cell->active_cell_index()] =
3289  cell->future_fe_index();
3290 
3291  // Create transfer object and attach to it.
3292  const auto *distributed_tria = dynamic_cast<
3294  &this->get_triangulation());
3295 
3296  active_fe_index_transfer->cell_data_transfer = std::make_unique<
3297  parallel::distributed::
3298  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3299  *distributed_tria,
3300  /*transfer_variable_size_data=*/false,
3301  /*refinement_strategy=*/
3302  &::AdaptationStrategies::Refinement::
3303  preserve<dim, spacedim, unsigned int>,
3304  /*coarsening_strategy=*/
3305  [this](
3306  const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3307  const std::vector<unsigned int> &children_fe_indices)
3308  -> unsigned int {
3309  return ::internal::hp::DoFHandlerImplementation::
3310  Implementation::determine_fe_from_children<dim, spacedim>(
3311  parent, children_fe_indices, fe_collection);
3312  });
3313 
3314  active_fe_index_transfer->cell_data_transfer
3315  ->prepare_for_coarsening_and_refinement(
3316  active_fe_index_transfer->active_fe_indices);
3317  }
3318 #endif
3319 }
3320 
3321 
3322 
3323 template <int dim, int spacedim>
3324 void
3326 {
3327  // Finite elements need to be assigned to each cell by calling
3328  // distribute_dofs() first to make this functionality available.
3329  if (this->fe_collection.size() > 0)
3330  {
3331  Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
3332 
3335 
3336  // We have to distribute the information about active_fe_indices
3337  // of all cells (including the artificial ones) on all processors,
3338  // if a parallel::shared::Triangulation has been used.
3341 
3342  // Free memory.
3343  this->active_fe_index_transfer.reset();
3344  }
3345 }
3346 
3347 
3348 
3349 template <int dim, int spacedim>
3350 void
3352 {
3353 #ifndef DEAL_II_WITH_P4EST
3354  Assert(false, ExcInternalError());
3355 #else
3356  // Finite elements need to be assigned to each cell by calling
3357  // distribute_dofs() first to make this functionality available.
3358  if (this->fe_collection.size() > 0)
3359  {
3360  Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
3361 
3362  // Unpack active_fe_indices.
3363  this->active_fe_index_transfer->active_fe_indices.resize(
3366  this->active_fe_index_transfer->cell_data_transfer->unpack(
3367  this->active_fe_index_transfer->active_fe_indices);
3368 
3369  // Update all locally owned active_fe_indices.
3370  this->set_active_fe_indices(
3371  this->active_fe_index_transfer->active_fe_indices);
3372 
3373  // Update active_fe_indices on ghost cells.
3376 
3377  // Free memory.
3378  this->active_fe_index_transfer.reset();
3379  }
3380 #endif
3381 }
3382 
3383 
3384 
3385 template <int dim, int spacedim>
3386 void
3388 {
3389 #ifndef DEAL_II_WITH_P4EST
3390  Assert(false,
3391  ExcMessage(
3392  "You are attempting to use a functionality that is only available "
3393  "if deal.II was configured to use p4est, but cmake did not find a "
3394  "valid p4est library."));
3395 #else
3396  // the implementation below requires a p:d:T currently
3397  Assert(
3399  &this->get_triangulation()) != nullptr),
3400  ExcNotImplemented());
3401 
3402  // Finite elements need to be assigned to each cell by calling
3403  // distribute_dofs() first to make this functionality available.
3404  if (fe_collection.size() > 0)
3405  {
3407 
3408  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3409 
3410  // Create transfer object and attach to it.
3411  const auto *distributed_tria = dynamic_cast<
3413  &this->get_triangulation());
3414 
3415  active_fe_index_transfer->cell_data_transfer = std::make_unique<
3416  parallel::distributed::
3417  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3418  *distributed_tria,
3419  /*transfer_variable_size_data=*/false,
3420  /*refinement_strategy=*/
3421  &::AdaptationStrategies::Refinement::
3422  preserve<dim, spacedim, unsigned int>,
3423  /*coarsening_strategy=*/
3424  [this](
3425  const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3426  const std::vector<unsigned int> &children_fe_indices)
3427  -> unsigned int {
3428  return ::internal::hp::DoFHandlerImplementation::
3429  Implementation::determine_fe_from_children<dim, spacedim>(
3430  parent, children_fe_indices, fe_collection);
3431  });
3432 
3433  // If we work on a p::d::Triangulation, we have to transfer all
3434  // active fe indices since ownership of cells may change.
3435 
3436  // Gather all current active_fe_indices
3437  get_active_fe_indices(active_fe_index_transfer->active_fe_indices);
3438 
3439  // Attach to transfer object
3440  active_fe_index_transfer->cell_data_transfer->prepare_for_serialization(
3441  active_fe_index_transfer->active_fe_indices);
3442  }
3443 #endif
3444 }
3445 
3446 
3447 
3448 template <int dim, int spacedim>
3449 void
3451 {
3452 #ifndef DEAL_II_WITH_P4EST
3453  Assert(false,
3454  ExcMessage(
3455  "You are attempting to use a functionality that is only available "
3456  "if deal.II was configured to use p4est, but cmake did not find a "
3457  "valid p4est library."));
3458 #else
3459  if (this->fe_collection.size() > 0)
3460  {
3461  Assert(this->active_fe_index_transfer != nullptr, ExcInternalError());
3462 
3463  // Free memory.
3464  this->active_fe_index_transfer.reset();
3465  }
3466 #endif
3467 }
3468 
3469 
3470 
3471 template <int dim, int spacedim>
3472 void
3474 {
3475 #ifndef DEAL_II_WITH_P4EST
3476  Assert(false,
3477  ExcMessage(
3478  "You are attempting to use a functionality that is only available "
3479  "if deal.II was configured to use p4est, but cmake did not find a "
3480  "valid p4est library."));
3481 #else
3482  // the implementation below requires a p:d:T currently
3483  Assert(
3485  &this->get_triangulation()) != nullptr),
3486  ExcNotImplemented());
3487 
3488  // Finite elements need to be assigned to each cell by calling
3489  // distribute_dofs() first to make this functionality available.
3490  if (fe_collection.size() > 0)
3491  {
3493 
3494  active_fe_index_transfer = std::make_unique<ActiveFEIndexTransfer>();
3495 
3496  // Create transfer object and attach to it.
3497  const auto *distributed_tria = dynamic_cast<
3499  &this->get_triangulation());
3500 
3501  active_fe_index_transfer->cell_data_transfer = std::make_unique<
3502  parallel::distributed::
3503  CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>(
3504  *distributed_tria,
3505  /*transfer_variable_size_data=*/false,
3506  /*refinement_strategy=*/
3507  &::AdaptationStrategies::Refinement::
3508  preserve<dim, spacedim, unsigned int>,
3509  /*coarsening_strategy=*/
3510  [this](
3511  const typename Triangulation<dim, spacedim>::cell_iterator &parent,
3512  const std::vector<unsigned int> &children_fe_indices)
3513  -> unsigned int {
3514  return ::internal::hp::DoFHandlerImplementation::
3515  Implementation::determine_fe_from_children<dim, spacedim>(
3516  parent, children_fe_indices, fe_collection);
3517  });
3518 
3519  // Unpack active_fe_indices.
3520  active_fe_index_transfer->active_fe_indices.resize(
3522  active_fe_index_transfer->cell_data_transfer->deserialize(
3523  active_fe_index_transfer->active_fe_indices);
3524 
3525  // Update all locally owned active_fe_indices.
3526  set_active_fe_indices(active_fe_index_transfer->active_fe_indices);
3527 
3528  // Update active_fe_indices on ghost cells.
3531 
3532  // Free memory.
3533  active_fe_index_transfer.reset();
3534  }
3535 #endif
3536 }
3537 
3538 
3539 
3540 template <int dim, int spacedim>
3542  : coarsest_level(numbers::invalid_unsigned_int)
3543  , finest_level(0)
3544 {}
3545 
3546 
3547 
3548 template <int dim, int spacedim>
3549 void
3551  const unsigned int cl,
3552  const unsigned int fl,
3553  const unsigned int dofs_per_vertex)
3554 {
3555  coarsest_level = cl;
3556  finest_level = fl;
3557 
3559  {
3560  const unsigned int n_levels = finest_level - coarsest_level + 1;
3561  const unsigned int n_indices = n_levels * dofs_per_vertex;
3562 
3563  indices = std::make_unique<types::global_dof_index[]>(n_indices);
3564  std::fill(indices.get(),
3565  indices.get() + n_indices,
3567  }
3568  else
3569  indices.reset();
3570 }
3571 
3572 
3573 
3574 template <int dim, int spacedim>
3575 unsigned int
3577 {
3578  return coarsest_level;
3579 }
3580 
3581 
3582 
3583 template <int dim, int spacedim>
3584 unsigned int
3586 {
3587  return finest_level;
3588 }
3589 
3590 /*-------------- Explicit Instantiations -------------------------------*/
3591 #include "dof_handler.inst"
3592 
3593 
3594 
static ::ExceptionBase & ExcNotAvailableWithoutHP()
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
const bool hp_capability_enabled
Definition: dof_handler.h:425
unsigned int max_couplings_between_dofs() const
unsigned int n_active_cells() const
Definition: tria.cc:11816
unsigned int get_coarsest_level() const
static void reserve_space_mg(DoFHandler< 2, spacedim > &dof_handler)
Definition: dof_handler.cc:537
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition: dof_handler.h:1511
static void reserve_space(DoFHandler< 3, spacedim > &dof_handler)
Definition: dof_handler.cc:360
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:1491
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1580
static void reserve_space_mg(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:468
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:899
unsigned short int active_fe_index_type
Definition: dof_handler.h:440
static ::ExceptionBase & ExcInvalidFEIndex(int arg1, int arg2)
Task< RT > new_task(const std::function< RT()> &function)
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:761
const std::vector< types::subdomain_id > & get_true_subdomain_ids_of_cells() const
Definition: shared_tria.cc:333
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:11326
cell_iterator end() const
::internal::DoFHandlerImplementation::NumberCache number_cache
Definition: dof_handler.h:1439
static void reserve_space(::DoFHandler< 3, spacedim > &dof_handler)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1648
void clear()
std::vector< std::array< std::vector< offset_type >, dim+1 > > object_dof_ptr
Definition: dof_handler.h:1478
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
void post_distributed_active_fe_index_transfer()
IteratorRange< level_cell_iterator > mg_cell_iterators_on_level(const unsigned int level) const
void clear_mg_space()
unsigned int size() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1533
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:1518
void pre_active_fe_index_transfer()
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:1505
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:875
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:1453
void setup_policy_and_listeners()
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:947
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
active_cell_iterator begin_active(const unsigned int level=0) const
BlockInfo block_info_object
Definition: dof_handler.h:1409
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNoFESelected()
static void reserve_space(::DoFHandler< 1, spacedim > &dof_handler)
void pre_refinement_action()
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:851
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:719
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:1445
void distribute_mg_dofs()
#define Assert(cond, exc)
Definition: exceptions.h:1423
hp::FECollection< dim, spacedim > fe_collection
Definition: dof_handler.h:1422
void post_active_fe_index_transfer()
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void pre_distributed_active_fe_index_transfer()
types::global_dof_index n_dofs() const
std::unique_ptr< ActiveFEIndexTransfer > active_fe_index_transfer
Definition: dof_handler.h:1530
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition: dof_handler.h:1467
active_cell_iterator end_active(const unsigned int level) const
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
Definition: dof_handler.h:1524
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:369
void post_refinement_action()
static void reserve_space_release_space(DoFHandler< dim, spacedim > &dof_handler)
unsigned int level
Definition: grid_out.cc:4337
SmartPointer< const Triangulation< dim, spacedim >, DoFHandler< dim, spacedim > > tria
Definition: dof_handler.h:1415
typename LevelSelector::cell_iterator level_cell_iterator
Definition: dof_handler.h:407
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:1430
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:805
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:613
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:11766
TriangulationBase< dim, spacedim > Triangulation
Definition: tria_base.h:302
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:474
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1182
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:1498
void post_distributed_serialization_of_active_fe_indices()
DoFHandler(const bool hp_capability_enabled=false)
std::array< std::vector< active_fe_index_type >, dim+1 > hp_object_fe_indices
Definition: dof_handler.h:1486
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 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=[](const typename MeshType::active_cell_iterator &) { return true;})
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
std::unique_ptr< types::global_dof_index[]> indices
Definition: dof_handler.h:1359
void initialize_local_block_info()
static void reserve_space(DoFHandler< 1, spacedim > &dof_handler)
Definition: dof_handler.cc:260
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:451
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:783
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:301
const Triangulation< dim, spacedim > & get_triangulation() const
unsigned int coarsest_level
Definition: dof_handler.h:1343
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:368
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:923
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)
Definition: dof_handler.cc:988
void init(const unsigned int coarsest_level, const unsigned int finest_level, const unsigned int dofs_per_vertex)
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:1326
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:1763
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:371
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:827
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:343
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:1536
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:696
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:739
unsigned int boundary_id
Definition: types.h:129
size_type n_elements() const
Definition: index_set.h:1830
static ::ExceptionBase & ExcNoFESelected()
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
IteratorRange< cell_iterator > cell_iterators() const
static void reserve_space_faces(DoFHandler< dim, spacedim > &dof_handler)
types::global_dof_index n_locally_owned_dofs() const
std::vector< std::vector< offset_type > > cell_dof_cache_ptr
Definition: dof_handler.h:1459
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()