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