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