deal.II version GIT relicensing-2850-g1646a780ac 2025-03-17 15:10:00+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: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/config.h>
16
19#include <deal.II/base/mpi.templates.h>
21
22#include <deal.II/distributed/cell_data_transfer.templates.h>
26
29
32#include <deal.II/grid/tria.h>
36
37#include <algorithm>
38#include <memory>
39#include <set>
40#include <unordered_set>
41
43
44#ifndef DOXYGEN
45template <int dim, int spacedim>
48#endif
49
50namespace internal
51{
52 template <int dim, int spacedim>
53 std::string
54 policy_to_string(const ::internal::DoFHandlerImplementation::Policy::
55 PolicyBase<dim, spacedim> &policy)
56 {
57 std::string policy_name;
58 if (dynamic_cast<const typename ::internal::DoFHandlerImplementation::
59 Policy::Sequential<dim, spacedim> *>(&policy))
60 policy_name = "Policy::Sequential<";
61 else if (dynamic_cast<
62 const typename ::internal::DoFHandlerImplementation::
63 Policy::ParallelDistributed<dim, spacedim> *>(&policy))
64 policy_name = "Policy::ParallelDistributed<";
65 else if (dynamic_cast<
66 const typename ::internal::DoFHandlerImplementation::
67 Policy::ParallelShared<dim, spacedim> *>(&policy))
68 policy_name = "Policy::ParallelShared<";
69 else
71 policy_name += Utilities::int_to_string(dim) + "," +
72 Utilities::int_to_string(spacedim) + ">";
73 return policy_name;
74 }
75
76
77 namespace DoFHandlerImplementation
78 {
84 {
89 template <int spacedim>
90 static unsigned int
92 {
93 return std::min(static_cast<types::global_dof_index>(
94 3 * dof_handler.fe_collection.max_dofs_per_vertex() +
95 2 * dof_handler.fe_collection.max_dofs_per_line()),
96 dof_handler.n_dofs());
97 }
98
99 template <int spacedim>
100 static unsigned int
102 {
103 // get these numbers by drawing pictures
104 // and counting...
105 // example:
106 // | | |
107 // --x-----x--x--X--
108 // | | | |
109 // | x--x--x
110 // | | | |
111 // --x--x--*--x--x--
112 // | | | |
113 // x--x--x |
114 // | | | |
115 // --X--x--x-----x--
116 // | | |
117 // x = vertices connected with center vertex *;
118 // = total of 19
119 // (the X vertices are connected with * if
120 // the vertices adjacent to X are hanging
121 // nodes)
122 // count lines -> 28 (don't forget to count
123 // mother and children separately!)
124 types::global_dof_index max_couplings;
125 switch (dof_handler.tria->max_adjacent_cells())
126 {
127 case 4:
128 max_couplings =
129 19 * dof_handler.fe_collection.max_dofs_per_vertex() +
130 28 * dof_handler.fe_collection.max_dofs_per_line() +
131 8 * dof_handler.fe_collection.max_dofs_per_quad();
132 break;
133 case 5:
134 max_couplings =
135 21 * dof_handler.fe_collection.max_dofs_per_vertex() +
136 31 * dof_handler.fe_collection.max_dofs_per_line() +
137 9 * dof_handler.fe_collection.max_dofs_per_quad();
138 break;
139 case 6:
140 max_couplings =
141 28 * dof_handler.fe_collection.max_dofs_per_vertex() +
142 42 * dof_handler.fe_collection.max_dofs_per_line() +
143 12 * dof_handler.fe_collection.max_dofs_per_quad();
144 break;
145 case 7:
146 max_couplings =
147 30 * dof_handler.fe_collection.max_dofs_per_vertex() +
148 45 * dof_handler.fe_collection.max_dofs_per_line() +
149 13 * dof_handler.fe_collection.max_dofs_per_quad();
150 break;
151 case 8:
152 max_couplings =
153 37 * dof_handler.fe_collection.max_dofs_per_vertex() +
154 56 * dof_handler.fe_collection.max_dofs_per_line() +
155 16 * dof_handler.fe_collection.max_dofs_per_quad();
156 break;
157
158 // the following numbers are not based on actual counting but by
159 // extrapolating the number sequences from the previous ones (for
160 // example, for n_dofs_per_vertex(), the sequence above is 19, 21,
161 // 28, 30, 37, and is continued as follows):
162 case 9:
163 max_couplings =
164 39 * dof_handler.fe_collection.max_dofs_per_vertex() +
165 59 * dof_handler.fe_collection.max_dofs_per_line() +
166 17 * dof_handler.fe_collection.max_dofs_per_quad();
167 break;
168 case 10:
169 max_couplings =
170 46 * dof_handler.fe_collection.max_dofs_per_vertex() +
171 70 * dof_handler.fe_collection.max_dofs_per_line() +
172 20 * dof_handler.fe_collection.max_dofs_per_quad();
173 break;
174 case 11:
175 max_couplings =
176 48 * dof_handler.fe_collection.max_dofs_per_vertex() +
177 73 * dof_handler.fe_collection.max_dofs_per_line() +
178 21 * dof_handler.fe_collection.max_dofs_per_quad();
179 break;
180 case 12:
181 max_couplings =
182 55 * dof_handler.fe_collection.max_dofs_per_vertex() +
183 84 * dof_handler.fe_collection.max_dofs_per_line() +
184 24 * dof_handler.fe_collection.max_dofs_per_quad();
185 break;
186 case 13:
187 max_couplings =
188 57 * dof_handler.fe_collection.max_dofs_per_vertex() +
189 87 * dof_handler.fe_collection.max_dofs_per_line() +
190 25 * dof_handler.fe_collection.max_dofs_per_quad();
191 break;
192 case 14:
193 max_couplings =
194 63 * dof_handler.fe_collection.max_dofs_per_vertex() +
195 98 * dof_handler.fe_collection.max_dofs_per_line() +
196 28 * dof_handler.fe_collection.max_dofs_per_quad();
197 break;
198 case 15:
199 max_couplings =
200 65 * dof_handler.fe_collection.max_dofs_per_vertex() +
201 103 * dof_handler.fe_collection.max_dofs_per_line() +
202 29 * dof_handler.fe_collection.max_dofs_per_quad();
203 break;
204 case 16:
205 max_couplings =
206 72 * dof_handler.fe_collection.max_dofs_per_vertex() +
207 114 * dof_handler.fe_collection.max_dofs_per_line() +
208 32 * dof_handler.fe_collection.max_dofs_per_quad();
209 break;
210
211 default:
213 max_couplings = 0;
214 }
215 return std::min(max_couplings, dof_handler.n_dofs());
216 }
217
218 template <int spacedim>
219 static unsigned int
221 {
222 // TODO:[?] Invent significantly better estimates than the ones in this
223 // function
224
225 // doing the same thing here is a rather complicated thing, compared
226 // to the 2d case, since it is hard to draw pictures with several
227 // refined hexahedra :-) so I presently only give a coarse
228 // estimate for the case that at most 8 hexes meet at each vertex
229 //
230 // can anyone give better estimate here?
231 const unsigned int max_adjacent_cells =
232 dof_handler.tria->max_adjacent_cells();
233
234 types::global_dof_index max_couplings;
235 if (max_adjacent_cells <= 8)
236 max_couplings =
237 7 * 7 * 7 * dof_handler.fe_collection.max_dofs_per_vertex() +
238 7 * 6 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_line() +
239 9 * 4 * 7 * 3 * dof_handler.fe_collection.max_dofs_per_quad() +
240 27 * dof_handler.fe_collection.max_dofs_per_hex();
241 else
242 {
244 max_couplings = 0;
245 }
246
247 return std::min(max_couplings, dof_handler.n_dofs());
248 }
249
254 template <int dim, int spacedim>
255 static void
257 {
258 dof_handler.object_dof_indices.clear();
259 dof_handler.object_dof_indices.resize(dof_handler.tria->n_levels());
260 dof_handler.object_dof_indices.shrink_to_fit();
261
262 dof_handler.object_dof_ptr.clear();
263 dof_handler.object_dof_ptr.resize(dof_handler.tria->n_levels());
264 dof_handler.object_dof_ptr.shrink_to_fit();
265 }
266
270 template <int dim, int spacedim>
271 static void
273 const unsigned int n_inner_dofs_per_cell)
274 {
275 for (unsigned int i = 0; i < dof_handler.tria->n_levels(); ++i)
276 {
277 dof_handler.object_dof_ptr[i][dim].assign(
278 dof_handler.tria->n_raw_cells(i) + 1, 0);
279
280 for (const auto &cell :
281 dof_handler.tria->cell_iterators_on_level(i))
282 if (cell->is_active() && !cell->is_artificial())
283 dof_handler.object_dof_ptr[i][dim][cell->index() + 1] =
284 n_inner_dofs_per_cell;
285
286 for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i); ++j)
287 dof_handler.object_dof_ptr[i][dim][j + 1] +=
288 dof_handler.object_dof_ptr[i][dim][j];
289
290 dof_handler.object_dof_indices[i][dim].resize(
291 dof_handler.object_dof_ptr[i][dim].back(),
293 }
294 }
295
300 template <int dim, int spacedim, typename T>
301 static void
303 const unsigned int structdim,
304 const unsigned int n_raw_entities,
305 const T &cell_process)
306 {
307 if (dof_handler.tria->n_cells() == 0)
308 return;
309
310 dof_handler.object_dof_ptr[0][structdim].assign(n_raw_entities + 1, -1);
311 // determine for each entity the number of dofs
312 for (const auto &cell : dof_handler.tria->cell_iterators())
313 if (cell->is_active() && !cell->is_artificial())
314 cell_process(
315 cell,
316 [&](const unsigned int n_dofs_per_entity,
317 const unsigned int index) {
318 auto &n_dofs_per_entity_target =
319 dof_handler.object_dof_ptr[0][structdim][index + 1];
320
321 // make sure that either the entity has not been visited or
322 // the entity has the same number of dofs assigned
323 Assert((n_dofs_per_entity_target ==
324 static_cast<
326 -1) ||
327 n_dofs_per_entity_target == n_dofs_per_entity),
329
330 n_dofs_per_entity_target = n_dofs_per_entity;
331 });
332
333 // convert the absolute numbers to CRS
334 dof_handler.object_dof_ptr[0][structdim][0] = 0;
335 for (unsigned int i = 1; i < n_raw_entities + 1; ++i)
336 {
337 if (dof_handler.object_dof_ptr[0][structdim][i] ==
338 static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
339 -1))
340 dof_handler.object_dof_ptr[0][structdim][i] =
341 dof_handler.object_dof_ptr[0][structdim][i - 1];
342 else
343 dof_handler.object_dof_ptr[0][structdim][i] +=
344 dof_handler.object_dof_ptr[0][structdim][i - 1];
345 }
346
347 // allocate memory for indices
348 dof_handler.object_dof_indices[0][structdim].resize(
349 dof_handler.object_dof_ptr[0][structdim].back(),
351 }
352
359 template <int dim, int spacedim>
360 static void
362 {
363 reset_to_empty_objects(dof_handler);
364
365 const auto &fe = dof_handler.get_fe();
366
367 // cell
368 reserve_cells(dof_handler,
369 dim == 1 ? fe.n_dofs_per_line() :
370 (dim == 2 ? fe.n_dofs_per_quad(0) :
371 fe.n_dofs_per_hex()));
372
373 // vertices
374 reserve_subentities(dof_handler,
375 0,
376 dof_handler.tria->n_vertices(),
377 [&](const auto &cell, const auto &process) {
378 for (const auto vertex_index :
379 cell->vertex_indices())
380 process(fe.n_dofs_per_vertex(),
381 cell->vertex_index(vertex_index));
382 });
383
384 // lines
385 if (dim == 2 || dim == 3)
386 reserve_subentities(dof_handler,
387 1,
388 dof_handler.tria->n_raw_lines(),
389 [&](const auto &cell, const auto &process) {
390 for (const auto line_index :
391 cell->line_indices())
392 process(fe.n_dofs_per_line(),
393 cell->line(line_index)->index());
394 });
395
396 // quads
397 if (dim == 3)
398 reserve_subentities(dof_handler,
399 2,
400 dof_handler.tria->n_raw_quads(),
401 [&](const auto &cell, const auto &process) {
402 for (const auto face_index :
403 cell->face_indices())
404 process(fe.n_dofs_per_quad(face_index),
405 cell->face(face_index)->index());
406 });
407 }
408
409 template <int spacedim>
410 static void
412 {
413 Assert(dof_handler.get_triangulation().n_levels() > 0,
414 ExcMessage("Invalid triangulation"));
415 dof_handler.clear_mg_space();
416
417 const ::Triangulation<1, spacedim> &tria =
418 dof_handler.get_triangulation();
419 const unsigned int dofs_per_line =
420 dof_handler.get_fe().n_dofs_per_line();
421 const unsigned int n_levels = tria.n_levels();
422
423 for (unsigned int i = 0; i < n_levels; ++i)
424 {
425 dof_handler.mg_levels.emplace_back(
427 dof_handler.mg_levels.back()->dof_object.dofs =
428 std::vector<types::global_dof_index>(tria.n_raw_lines(i) *
429 dofs_per_line,
431 }
432
433 const unsigned int n_vertices = tria.n_vertices();
434
435 dof_handler.mg_vertex_dofs.resize(n_vertices);
436
437 std::vector<unsigned int> max_level(n_vertices, 0);
438 std::vector<unsigned int> min_level(n_vertices, n_levels);
439
440 for (typename ::Triangulation<1, spacedim>::cell_iterator cell =
441 tria.begin();
442 cell != tria.end();
443 ++cell)
444 {
445 const unsigned int level = cell->level();
446
447 for (const auto vertex : cell->vertex_indices())
448 {
449 const unsigned int vertex_index = cell->vertex_index(vertex);
450
451 if (min_level[vertex_index] > level)
452 min_level[vertex_index] = level;
453
454 if (max_level[vertex_index] < level)
455 max_level[vertex_index] = level;
456 }
457 }
458
459 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
460 if (tria.vertex_used(vertex))
461 {
462 Assert(min_level[vertex] < n_levels, ExcInternalError());
463 Assert(max_level[vertex] >= min_level[vertex],
465 dof_handler.mg_vertex_dofs[vertex].init(
466 min_level[vertex],
467 max_level[vertex],
468 dof_handler.get_fe().n_dofs_per_vertex());
469 }
470
471 else
472 {
473 Assert(min_level[vertex] == n_levels, ExcInternalError());
474 Assert(max_level[vertex] == 0, ExcInternalError());
475 dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
476 }
477 }
478
479 template <int spacedim>
480 static void
482 {
483 Assert(dof_handler.get_triangulation().n_levels() > 0,
484 ExcMessage("Invalid triangulation"));
485 dof_handler.clear_mg_space();
486
487 const ::FiniteElement<2, spacedim> &fe = dof_handler.get_fe();
488 const ::Triangulation<2, spacedim> &tria =
489 dof_handler.get_triangulation();
490 const unsigned int n_levels = tria.n_levels();
491
492 for (unsigned int i = 0; i < n_levels; ++i)
493 {
494 dof_handler.mg_levels.emplace_back(
495 std::make_unique<
497 dof_handler.mg_levels.back()->dof_object.dofs =
498 std::vector<types::global_dof_index>(
499 tria.n_raw_quads(i) *
500 fe.n_dofs_per_quad(0 /*note: in 2d there is only one quad*/),
502 }
503
504 dof_handler.mg_faces =
505 std::make_unique<internal::DoFHandlerImplementation::DoFFaces<2>>();
506 dof_handler.mg_faces->lines.dofs =
507 std::vector<types::global_dof_index>(tria.n_raw_lines() *
508 fe.n_dofs_per_line(),
510
511 const unsigned int n_vertices = tria.n_vertices();
512
513 dof_handler.mg_vertex_dofs.resize(n_vertices);
514
515 std::vector<unsigned int> max_level(n_vertices, 0);
516 std::vector<unsigned int> min_level(n_vertices, n_levels);
517
518 for (typename ::Triangulation<2, spacedim>::cell_iterator cell =
519 tria.begin();
520 cell != tria.end();
521 ++cell)
522 {
523 const unsigned int level = cell->level();
524
525 for (const auto vertex : cell->vertex_indices())
526 {
527 const unsigned int vertex_index = cell->vertex_index(vertex);
528
529 if (min_level[vertex_index] > level)
530 min_level[vertex_index] = level;
531
532 if (max_level[vertex_index] < level)
533 max_level[vertex_index] = level;
534 }
535 }
536
537 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
538 if (tria.vertex_used(vertex))
539 {
540 Assert(min_level[vertex] < n_levels, ExcInternalError());
541 Assert(max_level[vertex] >= min_level[vertex],
543 dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
544 max_level[vertex],
545 fe.n_dofs_per_vertex());
546 }
547
548 else
549 {
550 Assert(min_level[vertex] == n_levels, ExcInternalError());
551 Assert(max_level[vertex] == 0, ExcInternalError());
552 dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
553 }
554 }
555
556 template <int spacedim>
557 static void
559 {
560 Assert(dof_handler.get_triangulation().n_levels() > 0,
561 ExcMessage("Invalid triangulation"));
562 dof_handler.clear_mg_space();
563
564 const ::FiniteElement<3, spacedim> &fe = dof_handler.get_fe();
565 const ::Triangulation<3, spacedim> &tria =
566 dof_handler.get_triangulation();
567 const unsigned int n_levels = tria.n_levels();
568
569 for (unsigned int i = 0; i < n_levels; ++i)
570 {
571 dof_handler.mg_levels.emplace_back(
572 std::make_unique<
574 dof_handler.mg_levels.back()->dof_object.dofs =
575 std::vector<types::global_dof_index>(tria.n_raw_hexs(i) *
576 fe.n_dofs_per_hex(),
578 }
579
580 dof_handler.mg_faces =
581 std::make_unique<internal::DoFHandlerImplementation::DoFFaces<3>>();
582 dof_handler.mg_faces->lines.dofs =
583 std::vector<types::global_dof_index>(tria.n_raw_lines() *
584 fe.n_dofs_per_line(),
586
587 // TODO: the implementation makes the assumption that all faces have the
588 // same number of dofs
589 AssertDimension(fe.n_unique_faces(), 1);
590 dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index>(
591 tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/),
593
594 const unsigned int n_vertices = tria.n_vertices();
595
596 dof_handler.mg_vertex_dofs.resize(n_vertices);
597
598 std::vector<unsigned int> max_level(n_vertices, 0);
599 std::vector<unsigned int> min_level(n_vertices, n_levels);
600
601 for (typename ::Triangulation<3, spacedim>::cell_iterator cell =
602 tria.begin();
603 cell != tria.end();
604 ++cell)
605 {
606 const unsigned int level = cell->level();
607
608 for (const auto vertex : cell->vertex_indices())
609 {
610 const unsigned int vertex_index = cell->vertex_index(vertex);
611
612 if (min_level[vertex_index] > level)
613 min_level[vertex_index] = level;
614
615 if (max_level[vertex_index] < level)
616 max_level[vertex_index] = level;
617 }
618 }
619
620 for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
621 if (tria.vertex_used(vertex))
622 {
623 Assert(min_level[vertex] < n_levels, ExcInternalError());
624 Assert(max_level[vertex] >= min_level[vertex],
626 dof_handler.mg_vertex_dofs[vertex].init(min_level[vertex],
627 max_level[vertex],
628 fe.n_dofs_per_vertex());
629 }
630
631 else
632 {
633 Assert(min_level[vertex] == n_levels, ExcInternalError());
634 Assert(max_level[vertex] == 0, ExcInternalError());
635 dof_handler.mg_vertex_dofs[vertex].init(1, 0, 0);
636 }
637 }
638 };
639 } // namespace DoFHandlerImplementation
640
641
642
643 namespace hp
644 {
645 namespace DoFHandlerImplementation
646 {
652 {
658 template <int dim, int spacedim>
659 static void
661 DoFHandler<dim, spacedim> &dof_handler)
662 {
663 (void)dof_handler;
664 for (const auto &cell : dof_handler.active_cell_iterators())
665 if (cell->is_locally_owned())
666 Assert(
667 !cell->future_fe_index_set(),
669 "There shouldn't be any cells flagged for p-adaptation when partitioning."));
670 }
671
672
673
678 template <int dim, int spacedim>
679 static void
681 {
682 // The final step in all of the reserve_space() functions is to set
683 // up vertex dof information. since vertices are sequentially
684 // numbered, what we do first is to set up an array in which
685 // we record whether a vertex is associated with any of the
686 // given fe's, by setting a bit. in a later step, we then
687 // actually allocate memory for the required dofs
688 //
689 // in the following, we only need to consider vertices that are
690 // adjacent to either a locally owned or a ghost cell; we never
691 // store anything on vertices that are only surrounded by
692 // artificial cells. so figure out that subset of vertices
693 // first
694 std::vector<bool> locally_used_vertices(
695 dof_handler.tria->n_vertices(), false);
696 for (const auto &cell : dof_handler.active_cell_iterators())
697 if (!cell->is_artificial())
698 for (const auto v : cell->vertex_indices())
699 locally_used_vertices[cell->vertex_index(v)] = true;
700
701 std::vector<std::vector<bool>> vertex_fe_association(
702 dof_handler.fe_collection.size(),
703 std::vector<bool>(dof_handler.tria->n_vertices(), false));
704
705 for (const auto &cell : dof_handler.active_cell_iterators())
706 if (!cell->is_artificial())
707 for (const auto v : cell->vertex_indices())
708 vertex_fe_association[cell->active_fe_index()]
709 [cell->vertex_index(v)] = true;
710
711 // in debug mode, make sure that each vertex is associated
712 // with at least one FE (note that except for unused
713 // vertices, all vertices are actually active). this is of
714 // course only true for vertices that are part of either
715 // ghost or locally owned cells
716 if constexpr (running_in_debug_mode())
717 {
718 for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
719 if (locally_used_vertices[v] == true)
720 if (dof_handler.tria->vertex_used(v) == true)
721 {
722 unsigned int fe = 0;
723 for (; fe < dof_handler.fe_collection.size(); ++fe)
724 if (vertex_fe_association[fe][v] == true)
725 break;
726 Assert(fe != dof_handler.fe_collection.size(),
728 }
729 }
730
731 const unsigned int d = 0;
732 const unsigned int l = 0;
733
734 dof_handler.hp_object_fe_ptr[d].clear();
735 dof_handler.hp_object_fe_indices[d].clear();
736 dof_handler.object_dof_ptr[l][d].clear();
737 dof_handler.object_dof_indices[l][d].clear();
738
739 dof_handler.hp_object_fe_ptr[d].reserve(
740 dof_handler.tria->n_vertices() + 1);
741
742 unsigned int vertex_slots_needed = 0;
743 unsigned int fe_slots_needed = 0;
744
745 for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
746 {
747 dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
748
749 if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
750 {
751 for (unsigned int fe = 0;
752 fe < dof_handler.fe_collection.size();
753 ++fe)
754 if (vertex_fe_association[fe][v] == true)
755 {
756 ++fe_slots_needed;
757 vertex_slots_needed +=
758 dof_handler.get_fe(fe).n_dofs_per_vertex();
759 }
760 }
761 }
762
763 dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
764
765 dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
766 dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
767
768 dof_handler.object_dof_indices[l][d].reserve(vertex_slots_needed);
769
770 for (unsigned int v = 0; v < dof_handler.tria->n_vertices(); ++v)
771 if (dof_handler.tria->vertex_used(v) && locally_used_vertices[v])
772 {
773 for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
774 ++fe)
775 if (vertex_fe_association[fe][v] == true)
776 {
777 dof_handler.hp_object_fe_indices[d].push_back(fe);
778 dof_handler.object_dof_ptr[l][d].push_back(
779 dof_handler.object_dof_indices[l][d].size());
780
781 for (unsigned int i = 0;
782 i < dof_handler.get_fe(fe).n_dofs_per_vertex();
783 i++)
784 dof_handler.object_dof_indices[l][d].push_back(
786 }
787 }
788
789
790 dof_handler.object_dof_ptr[l][d].push_back(
791 dof_handler.object_dof_indices[l][d].size());
792
793 AssertDimension(vertex_slots_needed,
794 dof_handler.object_dof_indices[l][d].size());
795 AssertDimension(fe_slots_needed,
796 dof_handler.hp_object_fe_indices[d].size());
797 AssertDimension(fe_slots_needed + 1,
798 dof_handler.object_dof_ptr[l][d].size());
799 AssertDimension(dof_handler.tria->n_vertices() + 1,
800 dof_handler.hp_object_fe_ptr[d].size());
801
802 dof_handler.object_dof_indices[l][d].assign(
803 vertex_slots_needed, numbers::invalid_dof_index);
804 }
805
806
807
812 template <int dim, int spacedim>
813 static void
815 {
816 (void)dof_handler;
817 // count how much space we need on each level for the cell
818 // dofs and set the dof_*_offsets data. initially set the
819 // latter to an invalid index, and only later set it to
820 // something reasonable for active dof_handler.cells
821 //
822 // note that for dof_handler.cells, the situation is simpler
823 // than for other (lower dimensional) objects since exactly
824 // one finite element is used for it
825 for (unsigned int level = 0; level < dof_handler.tria->n_levels();
826 ++level)
827 {
828 dof_handler.object_dof_ptr[level][dim] =
829 std::vector<typename DoFHandler<dim, spacedim>::offset_type>(
830 dof_handler.tria->n_raw_cells(level),
831 static_cast<typename DoFHandler<dim, spacedim>::offset_type>(
832 -1));
833
834 types::global_dof_index next_free_dof = 0;
835 for (auto cell :
837 if (cell->is_active() && !cell->is_artificial())
838 {
839 dof_handler.object_dof_ptr[level][dim][cell->index()] =
840 next_free_dof;
841 next_free_dof +=
842 cell->get_fe().template n_dofs_per_object<dim>();
843 }
844
845 dof_handler.object_dof_indices[level][dim] =
846 std::vector<types::global_dof_index>(
847 next_free_dof, numbers::invalid_dof_index);
848 }
849 }
850
851
852
857 template <int dim, int spacedim>
858 static void
860 {
861 // FACE DOFS
862 //
863 // Count face dofs, then allocate as much space
864 // as we need and prime the linked list for faces (see the
865 // description in hp::DoFLevel) with the indices we will
866 // need. Note that our task is more complicated than for the
867 // cell case above since two adjacent cells may have different
868 // active FE indices, in which case we need to allocate
869 // *two* sets of face dofs for the same face. But they don't
870 // *have* to be different, and so we need to prepare for this
871 // as well.
872 //
873 // The way we do things is that we loop over all active cells (these
874 // are the only ones that have DoFs anyway) and all their faces. We
875 // note in the vector face_touched whether we have previously
876 // visited a face and if so skip it
877 {
878 std::vector<bool> face_touched(dim == 2 ?
879 dof_handler.tria->n_raw_lines() :
880 dof_handler.tria->n_raw_quads());
881
882 const unsigned int d = dim - 1;
883 const unsigned int l = 0;
884
885 dof_handler.hp_object_fe_ptr[d].clear();
886 dof_handler.hp_object_fe_indices[d].clear();
887 dof_handler.object_dof_ptr[l][d].clear();
888 dof_handler.object_dof_indices[l][d].clear();
889
890 dof_handler.hp_object_fe_ptr[d].resize(
891 dof_handler.tria->n_raw_faces() + 1);
892
893 // An array to hold how many slots (see the hp::DoFLevel
894 // class) we will have to store on each level
895 unsigned int n_face_slots = 0;
896
897 for (const auto &cell : dof_handler.active_cell_iterators())
898 if (!cell->is_artificial())
899 for (const auto face : cell->face_indices())
900 if (!face_touched[cell->face(face)->index()])
901 {
902 unsigned int fe_slots_needed = 0;
903
904 if (cell->at_boundary(face) ||
905 cell->face(face)->has_children() ||
906 cell->neighbor_is_coarser(face) ||
907 (!cell->at_boundary(face) &&
908 cell->neighbor(face)->is_artificial()) ||
909 (!cell->at_boundary(face) &&
910 !cell->neighbor(face)->is_artificial() &&
911 (cell->active_fe_index() ==
912 cell->neighbor(face)->active_fe_index())))
913 {
914 fe_slots_needed = 1;
915 n_face_slots +=
916 dof_handler.get_fe(cell->active_fe_index())
917 .template n_dofs_per_object<dim - 1>(face);
918 }
919 else
920 {
921 fe_slots_needed = 2;
922 n_face_slots +=
923 dof_handler.get_fe(cell->active_fe_index())
924 .template n_dofs_per_object<dim - 1>(face) +
925 dof_handler
926 .get_fe(cell->neighbor(face)->active_fe_index())
927 .template n_dofs_per_object<dim - 1>(
928 cell->neighbor_face_no(face));
929 }
930
931 // mark this face as visited
932 face_touched[cell->face(face)->index()] = true;
933
934 dof_handler
935 .hp_object_fe_ptr[d][cell->face(face)->index() + 1] =
936 fe_slots_needed;
937 }
938
939 for (unsigned int i = 1; i < dof_handler.hp_object_fe_ptr[d].size();
940 i++)
941 dof_handler.hp_object_fe_ptr[d][i] +=
942 dof_handler.hp_object_fe_ptr[d][i - 1];
943
944
945 dof_handler.hp_object_fe_indices[d].resize(
946 dof_handler.hp_object_fe_ptr[d].back());
947 dof_handler.object_dof_ptr[l][d].resize(
948 dof_handler.hp_object_fe_ptr[d].back() + 1);
949
950 dof_handler.object_dof_indices[l][d].reserve(n_face_slots);
951
952
953 // With the memory now allocated, loop over the
954 // dof_handler cells again and prime the _offset values as
955 // well as the fe_index fields
956 face_touched = std::vector<bool>(face_touched.size());
957
958 for (const auto &cell : dof_handler.active_cell_iterators())
959 if (!cell->is_artificial())
960 for (const auto face : cell->face_indices())
961 if (!face_touched[cell->face(face)->index()])
962 {
963 // Same decision tree as before
964 if (cell->at_boundary(face) ||
965 cell->face(face)->has_children() ||
966 cell->neighbor_is_coarser(face) ||
967 (!cell->at_boundary(face) &&
968 cell->neighbor(face)->is_artificial()) ||
969 (!cell->at_boundary(face) &&
970 !cell->neighbor(face)->is_artificial() &&
971 (cell->active_fe_index() ==
972 cell->neighbor(face)->active_fe_index())))
973 {
974 const types::fe_index fe = cell->active_fe_index();
975 const unsigned int n_dofs =
976 dof_handler.get_fe(fe)
977 .template n_dofs_per_object<dim - 1>(face);
978 const unsigned int offset =
979 dof_handler
980 .hp_object_fe_ptr[d][cell->face(face)->index()];
981
982 dof_handler.hp_object_fe_indices[d][offset] = fe;
983 dof_handler.object_dof_ptr[l][d][offset + 1] = n_dofs;
984
985 for (unsigned int i = 0; i < n_dofs; ++i)
986 dof_handler.object_dof_indices[l][d].push_back(
988 }
989 else
990 {
991 types::fe_index fe_1 = cell->active_fe_index();
992 unsigned int face_no_1 = face;
993 types::fe_index fe_2 =
994 cell->neighbor(face)->active_fe_index();
995 unsigned int face_no_2 = cell->neighbor_face_no(face);
996
997 if (fe_2 < fe_1)
998 {
999 std::swap(fe_1, fe_2);
1000 std::swap(face_no_1, face_no_2);
1001 }
1002
1003 const unsigned int n_dofs_1 =
1004 dof_handler.get_fe(fe_1)
1005 .template n_dofs_per_object<dim - 1>(face_no_1);
1006
1007 const unsigned int n_dofs_2 =
1008 dof_handler.get_fe(fe_2)
1009 .template n_dofs_per_object<dim - 1>(face_no_2);
1010
1011 const unsigned int offset =
1012 dof_handler
1013 .hp_object_fe_ptr[d][cell->face(face)->index()];
1014
1015 dof_handler.hp_object_fe_indices[d].push_back(
1016 cell->active_fe_index());
1017 dof_handler.object_dof_ptr[l][d].push_back(
1018 dof_handler.object_dof_indices[l][d].size());
1019
1020 dof_handler.hp_object_fe_indices[d][offset + 0] =
1021 fe_1;
1022 dof_handler.hp_object_fe_indices[d][offset + 1] =
1023 fe_2;
1024 dof_handler.object_dof_ptr[l][d][offset + 1] =
1025 n_dofs_1;
1026 dof_handler.object_dof_ptr[l][d][offset + 2] =
1027 n_dofs_2;
1028
1029
1030 for (unsigned int i = 0; i < n_dofs_1 + n_dofs_2; ++i)
1031 dof_handler.object_dof_indices[l][d].push_back(
1033 }
1034
1035 // mark this face as visited
1036 face_touched[cell->face(face)->index()] = true;
1037 }
1038
1039 for (unsigned int i = 1;
1040 i < dof_handler.object_dof_ptr[l][d].size();
1041 i++)
1042 dof_handler.object_dof_ptr[l][d][i] +=
1043 dof_handler.object_dof_ptr[l][d][i - 1];
1044 }
1045 }
1046
1047
1048
1055 template <int spacedim>
1056 static void
1058 {
1059 Assert(dof_handler.fe_collection.size() > 0,
1061 Assert(dof_handler.tria->n_levels() > 0,
1062 ExcMessage("The current Triangulation must not be empty."));
1063 Assert(dof_handler.tria->n_levels() ==
1064 dof_handler.hp_cell_future_fe_indices.size(),
1066
1068 reset_to_empty_objects(dof_handler);
1069
1071 tasks +=
1072 Threads::new_task(&reserve_space_cells<1, spacedim>, dof_handler);
1073 tasks += Threads::new_task(&reserve_space_vertices<1, spacedim>,
1074 dof_handler);
1075 tasks.join_all();
1076 }
1077
1078
1079
1080 template <int spacedim>
1081 static void
1083 {
1084 Assert(dof_handler.fe_collection.size() > 0,
1086 Assert(dof_handler.tria->n_levels() > 0,
1087 ExcMessage("The current Triangulation must not be empty."));
1088 Assert(dof_handler.tria->n_levels() ==
1089 dof_handler.hp_cell_future_fe_indices.size(),
1091
1093 reset_to_empty_objects(dof_handler);
1094
1096 tasks +=
1097 Threads::new_task(&reserve_space_cells<2, spacedim>, dof_handler);
1098 tasks +=
1099 Threads::new_task(&reserve_space_faces<2, spacedim>, dof_handler);
1100 tasks += Threads::new_task(&reserve_space_vertices<2, spacedim>,
1101 dof_handler);
1102 tasks.join_all();
1103 }
1104
1105
1106
1107 template <int spacedim>
1108 static void
1110 {
1111 Assert(dof_handler.fe_collection.size() > 0,
1113 Assert(dof_handler.tria->n_levels() > 0,
1114 ExcMessage("The current Triangulation must not be empty."));
1115 Assert(dof_handler.tria->n_levels() ==
1116 dof_handler.hp_cell_future_fe_indices.size(),
1118
1120 reset_to_empty_objects(dof_handler);
1121
1123 tasks +=
1124 Threads::new_task(&reserve_space_cells<3, spacedim>, dof_handler);
1125 tasks +=
1126 Threads::new_task(&reserve_space_faces<3, spacedim>, dof_handler);
1127 tasks += Threads::new_task(&reserve_space_vertices<3, spacedim>,
1128 dof_handler);
1129
1130 // While the tasks above are running, we can turn to line dofs
1131
1132 // the situation here is pretty much like with vertices:
1133 // there can be an arbitrary number of finite elements
1134 // associated with each line.
1135 //
1136 // the algorithm we use is somewhat similar to what we do in
1137 // reserve_space_vertices()
1138 {
1139 // what we do first is to set up an array in which we
1140 // record whether a line is associated with any of the
1141 // given fe's, by setting a bit. in a later step, we
1142 // then actually allocate memory for the required dofs
1143 std::vector<std::vector<bool>> line_fe_association(
1144 dof_handler.fe_collection.size(),
1145 std::vector<bool>(dof_handler.tria->n_raw_lines(), false));
1146
1147 for (const auto &cell : dof_handler.active_cell_iterators())
1148 if (!cell->is_artificial())
1149 for (const auto l : cell->line_indices())
1150 line_fe_association[cell->active_fe_index()]
1151 [cell->line_index(l)] = true;
1152
1153 // first check which of the lines is used at all,
1154 // i.e. is associated with a finite element. we do this
1155 // since not all lines may actually be used, in which
1156 // case we do not have to allocate any memory at all
1157 std::vector<bool> line_is_used(dof_handler.tria->n_raw_lines(),
1158 false);
1159 for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1160 ++line)
1161 for (unsigned int fe = 0; fe < dof_handler.fe_collection.size();
1162 ++fe)
1163 if (line_fe_association[fe][line] == true)
1164 {
1165 line_is_used[line] = true;
1166 break;
1167 }
1168
1169
1170
1171 const unsigned int d = 1;
1172 const unsigned int l = 0;
1173
1174 dof_handler.hp_object_fe_ptr[d].clear();
1175 dof_handler.hp_object_fe_indices[d].clear();
1176 dof_handler.object_dof_ptr[l][d].clear();
1177 dof_handler.object_dof_indices[l][d].clear();
1178
1179 dof_handler.hp_object_fe_ptr[d].reserve(
1180 dof_handler.tria->n_raw_lines() + 1);
1181
1182 unsigned int line_slots_needed = 0;
1183 unsigned int fe_slots_needed = 0;
1184
1185 for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1186 ++line)
1187 {
1188 dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1189
1190 if (line_is_used[line] == true)
1191 {
1192 for (unsigned int fe = 0;
1193 fe < dof_handler.fe_collection.size();
1194 ++fe)
1195 if (line_fe_association[fe][line] == true)
1196 {
1197 ++fe_slots_needed;
1198 line_slots_needed +=
1199 dof_handler.get_fe(fe).n_dofs_per_line();
1200 }
1201 }
1202 }
1203
1204 dof_handler.hp_object_fe_ptr[d].push_back(fe_slots_needed);
1205
1206 // make sure that all entries have been set
1207 AssertDimension(dof_handler.hp_object_fe_ptr[d].size(),
1208 dof_handler.tria->n_raw_lines() + 1);
1209
1210 dof_handler.hp_object_fe_indices[d].reserve(fe_slots_needed);
1211 dof_handler.object_dof_ptr[l][d].reserve(fe_slots_needed + 1);
1212
1213 dof_handler.object_dof_indices[l][d].reserve(line_slots_needed);
1214
1215 for (unsigned int line = 0; line < dof_handler.tria->n_raw_lines();
1216 ++line)
1217 if (line_is_used[line] == true)
1218 {
1219 for (unsigned int fe = 0;
1220 fe < dof_handler.fe_collection.size();
1221 ++fe)
1222 if (line_fe_association[fe][line] == true)
1223 {
1224 dof_handler.hp_object_fe_indices[d].push_back(fe);
1225 dof_handler.object_dof_ptr[l][d].push_back(
1226 dof_handler.object_dof_indices[l][d].size());
1227
1228 for (unsigned int i = 0;
1229 i < dof_handler.get_fe(fe).n_dofs_per_line();
1230 i++)
1231 dof_handler.object_dof_indices[l][d].push_back(
1233 }
1234 }
1235
1236 dof_handler.object_dof_ptr[l][d].push_back(
1237 dof_handler.object_dof_indices[l][d].size());
1238
1239 // make sure that all entries have been set
1240 AssertDimension(dof_handler.hp_object_fe_indices[d].size(),
1241 fe_slots_needed);
1242 AssertDimension(dof_handler.object_dof_ptr[l][d].size(),
1243 fe_slots_needed + 1);
1244 AssertDimension(dof_handler.object_dof_indices[l][d].size(),
1245 line_slots_needed);
1246 }
1247
1248 // Ensure that everything is done at this point.
1249 tasks.join_all();
1250 }
1251
1252
1253
1265 template <int dim, int spacedim>
1266 static void
1268 {
1269 Assert(
1270 dof_handler.hp_capability_enabled == true,
1272
1273 if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1274 dynamic_cast<
1275 const ::parallel::shared::Triangulation<dim, spacedim>
1276 *>(&dof_handler.get_triangulation()))
1277 {
1278 // we have a shared triangulation. in this case, every processor
1279 // knows about all cells, but every processor only has knowledge
1280 // about the active FE index on the cells it owns.
1281 //
1282 // we can create a complete set of active FE indices by letting
1283 // every processor create a vector of indices for all cells,
1284 // filling only those on the cells it owns and setting the indices
1285 // on the other cells to zero. then we add all of these vectors
1286 // up, and because every vector entry has exactly one processor
1287 // that owns it, the sum is correct
1288 std::vector<types::fe_index> active_fe_indices(
1289 tr->n_active_cells(), 0u);
1290 for (const auto &cell : dof_handler.active_cell_iterators())
1291 if (cell->is_locally_owned())
1292 active_fe_indices[cell->active_cell_index()] =
1293 cell->active_fe_index();
1294
1295 Utilities::MPI::sum(active_fe_indices,
1296 tr->get_mpi_communicator(),
1297 active_fe_indices);
1298
1299 // now go back and fill the active FE index on all other
1300 // cells. we would like to call cell->set_active_fe_index(),
1301 // but that function does not allow setting these indices on
1302 // non-locally_owned cells. so we have to work around the
1303 // issue a little bit by accessing the underlying data
1304 // structures directly
1305 for (const auto &cell : dof_handler.active_cell_iterators())
1306 if (!cell->is_locally_owned())
1307 dof_handler
1308 .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1309 active_fe_indices[cell->active_cell_index()];
1310 }
1311 else if (const ::parallel::
1312 DistributedTriangulationBase<dim, spacedim> *tr =
1313 dynamic_cast<
1314 const ::parallel::
1315 DistributedTriangulationBase<dim, spacedim> *>(
1316 &dof_handler.get_triangulation()))
1317 {
1318 // For completely distributed meshes, use the function that is
1319 // able to move data from locally owned cells on one processor to
1320 // the corresponding ghost cells on others. To this end, we need
1321 // to have functions that can pack and unpack the data we want to
1322 // transport -- namely, the single unsigned int active_fe_index
1323 // objects
1324 auto pack =
1325 [](
1327 &cell) -> types::fe_index {
1328 return cell->active_fe_index();
1329 };
1330
1331 auto unpack =
1332 [&dof_handler](
1334 &cell,
1335 const types::fe_index active_fe_index) -> void {
1336 // we would like to say
1337 // cell->set_active_fe_index(active_fe_index);
1338 // but this is not allowed on cells that are not
1339 // locally owned, and we are on a ghost cell
1340 dof_handler
1341 .hp_cell_active_fe_indices[cell->level()][cell->index()] =
1342 active_fe_index;
1343 };
1344
1347 DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1348 }
1349 else
1350 {
1351 // a sequential triangulation. there is nothing we need to do here
1352 Assert(
1353 (dynamic_cast<
1354 const ::parallel::TriangulationBase<dim, spacedim> *>(
1355 &dof_handler.get_triangulation()) == nullptr),
1357 }
1358 }
1359
1360
1361
1375 template <int dim, int spacedim>
1376 static void
1378 {
1379 Assert(
1380 dof_handler.hp_capability_enabled == true,
1382
1383 if (const ::parallel::shared::Triangulation<dim, spacedim> *tr =
1384 dynamic_cast<
1385 const ::parallel::shared::Triangulation<dim, spacedim>
1386 *>(&dof_handler.get_triangulation()))
1387 {
1388 std::vector<types::fe_index> future_fe_indices(
1389 tr->n_active_cells(), 0u);
1390 for (const auto &cell : dof_handler.active_cell_iterators() |
1392 future_fe_indices[cell->active_cell_index()] =
1393 dof_handler
1394 .hp_cell_future_fe_indices[cell->level()][cell->index()];
1395
1396 Utilities::MPI::sum(future_fe_indices,
1397 tr->get_mpi_communicator(),
1398 future_fe_indices);
1399
1400 for (const auto &cell : dof_handler.active_cell_iterators())
1401 if (!cell->is_locally_owned())
1402 dof_handler
1403 .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1404 future_fe_indices[cell->active_cell_index()];
1405 }
1406 else if (const ::parallel::
1407 DistributedTriangulationBase<dim, spacedim> *tr =
1408 dynamic_cast<
1409 const ::parallel::
1410 DistributedTriangulationBase<dim, spacedim> *>(
1411 &dof_handler.get_triangulation()))
1412 {
1413 auto pack =
1414 [&dof_handler](
1416 &cell) -> types::fe_index {
1417 return dof_handler
1418 .hp_cell_future_fe_indices[cell->level()][cell->index()];
1419 };
1420
1421 auto unpack =
1422 [&dof_handler](
1424 &cell,
1425 const types::fe_index future_fe_index) -> void {
1426 dof_handler
1427 .hp_cell_future_fe_indices[cell->level()][cell->index()] =
1428 future_fe_index;
1429 };
1430
1433 DoFHandler<dim, spacedim>>(dof_handler, pack, unpack);
1434 }
1435 else
1436 {
1437 Assert(
1438 (dynamic_cast<
1439 const ::parallel::TriangulationBase<dim, spacedim> *>(
1440 &dof_handler.get_triangulation()) == nullptr),
1442 }
1443 }
1444
1445
1446
1467 template <int dim, int spacedim>
1468 static void
1470 DoFHandler<dim, spacedim> &dof_handler)
1471 {
1472 const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1473
1474 for (const auto &cell : dof_handler.active_cell_iterators())
1475 if (cell->is_locally_owned())
1476 {
1477 if (cell->refine_flag_set())
1478 {
1479 // Store the active FE index of each cell that will be
1480 // refined to and distribute it later on its children.
1481 // Pick their future index if flagged for p-refinement.
1482 fe_transfer->refined_cells_fe_index.insert(
1483 {cell, cell->future_fe_index()});
1484 }
1485 else if (cell->coarsen_flag_set())
1486 {
1487 // From all cells that will be coarsened, determine their
1488 // parent and calculate its proper active FE index, so that
1489 // it can be set after refinement. But first, check if that
1490 // particular cell has a parent at all.
1491 Assert(cell->level() > 0, ExcInternalError());
1492 const auto &parent = cell->parent();
1493
1494 // Check if the active FE index for the current cell has
1495 // been determined already.
1496 if (fe_transfer->coarsened_cells_fe_index.find(parent) ==
1497 fe_transfer->coarsened_cells_fe_index.end())
1498 {
1499 // Find a suitable active FE index for the parent cell
1500 // based on the 'least dominant finite element' of its
1501 // children. Consider the childrens' hypothetical future
1502 // index when they have been flagged for p-refinement.
1503 if constexpr (library_build_mode ==
1505 {
1506 for (const auto &child : parent->child_iterators())
1507 Assert(child->is_active() &&
1508 child->coarsen_flag_set(),
1509 typename ::Triangulation<
1510 dim>::ExcInconsistentCoarseningFlags());
1511 }
1512
1513 const types::fe_index fe_index = ::internal::hp::
1514 DoFHandlerImplementation::Implementation::
1515 dominated_future_fe_on_children<dim, spacedim>(
1516 parent);
1517
1518 fe_transfer->coarsened_cells_fe_index.insert(
1519 {parent, fe_index});
1520 }
1521 }
1522 else
1523 {
1524 // No h-refinement is scheduled for this cell.
1525 // However, it may have p-refinement indicators, so we
1526 // choose a new active FE index based on its flags.
1527 if (cell->future_fe_index_set() == true)
1528 fe_transfer->persisting_cells_fe_index.insert(
1529 {cell, cell->future_fe_index()});
1530 }
1531 }
1532 }
1533
1534
1535
1540 template <int dim, int spacedim>
1541 static void
1543 DoFHandler<dim, spacedim> &dof_handler)
1544 {
1545 const auto &fe_transfer = dof_handler.active_fe_index_transfer;
1546
1547 // Set active FE indices on persisting cells.
1548 for (const auto &persist : fe_transfer->persisting_cells_fe_index)
1549 {
1550 const auto &cell = persist.first;
1551
1552 if (cell->is_locally_owned())
1553 {
1554 Assert(cell->is_active(), ExcInternalError());
1555 cell->set_active_fe_index(persist.second);
1556 }
1557 }
1558
1559 // Distribute active FE indices from all refined cells on their
1560 // respective children.
1561 for (const auto &refine : fe_transfer->refined_cells_fe_index)
1562 {
1563 const auto &parent = refine.first;
1564
1565 for (const auto &child : parent->child_iterators())
1566 if (child->is_locally_owned())
1567 {
1568 Assert(child->is_active(), ExcInternalError());
1569 child->set_active_fe_index(refine.second);
1570 }
1571 }
1572
1573 // Set active FE indices on coarsened cells that have been determined
1574 // before the actual coarsening happened.
1575 for (const auto &coarsen : fe_transfer->coarsened_cells_fe_index)
1576 {
1577 const auto &cell = coarsen.first;
1578
1579 if (cell->is_locally_owned())
1580 {
1581 Assert(cell->is_active(), ExcInternalError());
1582 cell->set_active_fe_index(coarsen.second);
1583 }
1584 }
1585 }
1586
1587
1598 template <int dim, int spacedim>
1599 static types::fe_index
1602 const std::vector<types::fe_index> &children_fe_indices,
1603 const ::hp::FECollection<dim, spacedim> &fe_collection)
1604 {
1605 Assert(!children_fe_indices.empty(), ExcInternalError());
1606
1607 // convert vector to set
1608 // TODO: Change set to types::fe_index
1609 const std::set<unsigned int> children_fe_indices_set(
1610 children_fe_indices.begin(), children_fe_indices.end());
1611
1612 const types::fe_index dominated_fe_index =
1613 fe_collection.find_dominated_fe_extended(children_fe_indices_set,
1614 /*codim=*/0);
1615
1616 Assert(dominated_fe_index != numbers::invalid_fe_index,
1618
1619 return dominated_fe_index;
1620 }
1621
1622
1630 template <int dim, int spacedim>
1631 static types::fe_index
1633 const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
1634 {
1635 Assert(
1636 !parent->is_active(),
1637 ExcMessage(
1638 "You ask for information on children of this cell which is only "
1639 "available for active cells. This cell has no children."));
1640
1641 const auto &dof_handler = parent->get_dof_handler();
1642 Assert(
1643 dof_handler.has_hp_capabilities(),
1645
1646 // TODO: Change set to types::fe_index
1647 std::set<unsigned int> future_fe_indices_children;
1648 for (const auto &child : parent->child_iterators())
1649 {
1650 Assert(
1651 child->is_active(),
1652 ExcMessage(
1653 "You ask for information on children of this cell which is only "
1654 "available for active cells. One of its children is not active."));
1655
1656 // Ghost siblings might occur on parallel Triangulation
1657 // objects. The public interface does not allow to access future
1658 // FE indices on ghost cells. However, we need this information
1659 // here and thus call the internal function that does not check
1660 // for cell ownership. This requires that future FE indices have
1661 // been communicated prior to calling this function.
1662 const types::fe_index future_fe_index_child =
1663 ::internal::DoFCellAccessorImplementation::
1664 Implementation::future_fe_index<dim, spacedim, false>(*child);
1665
1666 future_fe_indices_children.insert(future_fe_index_child);
1667 }
1668 Assert(!future_fe_indices_children.empty(), ExcInternalError());
1669
1670 const types::fe_index future_fe_index =
1671 dof_handler.fe_collection.find_dominated_fe_extended(
1672 future_fe_indices_children,
1673 /*codim=*/0);
1674
1675 Assert(future_fe_index != numbers::invalid_fe_index,
1677
1678 return future_fe_index;
1679 }
1680 };
1681
1682
1683
1687 template <int dim, int spacedim>
1688 void
1690 {
1691 Implementation::communicate_future_fe_indices<dim, spacedim>(
1692 dof_handler);
1693 }
1694
1695
1696
1700 template <int dim, int spacedim>
1701 unsigned int
1703 const typename DoFHandler<dim, spacedim>::cell_iterator &parent)
1704 {
1705 return Implementation::dominated_future_fe_on_children<dim, spacedim>(
1706 parent);
1707 }
1708 } // namespace DoFHandlerImplementation
1709 } // namespace hp
1710} // namespace internal
1711
1712#ifndef DOXYGEN
1713
1714template <int dim, int spacedim>
1717 : hp_capability_enabled(true)
1718 , tria(nullptr, typeid(*this).name())
1719 , mg_faces(nullptr)
1720{}
1721
1722
1723
1724template <int dim, int spacedim>
1727 : DoFHandler()
1728{
1729 reinit(tria);
1730}
1731
1732
1733
1734template <int dim, int spacedim>
1737{
1738 // unsubscribe all attachments to signals of the underlying triangulation
1739 for (auto &connection : this->tria_listeners)
1740 connection.disconnect();
1741 this->tria_listeners.clear();
1742
1743 for (auto &connection : this->tria_listeners_for_transfer)
1744 connection.disconnect();
1745 this->tria_listeners_for_transfer.clear();
1746
1747 // release allocated memory
1748 // virtual functions called in constructors and destructors never use the
1749 // override in a derived class
1750 // for clarity be explicit on which function is called
1752
1753 // also release the policy. this needs to happen before the
1754 // current object disappears because the policy objects
1755 // store references to the DoFhandler object they work on
1756 this->policy.reset();
1757}
1758
1759
1760
1761template <int dim, int spacedim>
1764{
1765 //
1766 // call destructor
1767 //
1768 // remove association with old triangulation
1769 for (auto &connection : this->tria_listeners)
1770 connection.disconnect();
1771 this->tria_listeners.clear();
1772
1773 for (auto &connection : this->tria_listeners_for_transfer)
1774 connection.disconnect();
1775 this->tria_listeners_for_transfer.clear();
1776
1777 // release allocated memory and policy
1779 this->policy.reset();
1780
1781 // reset the finite element collection
1782 this->fe_collection = hp::FECollection<dim, spacedim>();
1783
1784 //
1785 // call constructor
1786 //
1787 // establish connection to new triangulation
1788 this->tria = &tria;
1789 this->setup_policy();
1790
1791 // start in hp-mode and let distribute_dofs toggle it if necessary
1792 hp_capability_enabled = true;
1793 this->connect_to_triangulation_signals();
1794 this->create_active_fe_table();
1795}
1796
1797#endif
1798/*------------------------ Cell iterator functions ------------------------*/
1799#ifndef DOXYGEN
1800template <int dim, int spacedim>
1803 DoFHandler<dim, spacedim>::begin(const unsigned int level) const
1804{
1806 this->get_triangulation().begin(level);
1807 if (cell == this->get_triangulation().end(level))
1808 return end(level);
1809 return cell_iterator(*cell, this);
1810}
1811
1812
1813
1814template <int dim, int spacedim>
1817 DoFHandler<dim, spacedim>::begin_active(const unsigned int level) const
1818{
1819 // level is checked in begin
1820 cell_iterator i = begin(level);
1821 if (i.state() != IteratorState::valid)
1822 return i;
1823 while (i->has_children())
1824 if ((++i).state() != IteratorState::valid)
1825 return i;
1826 return i;
1827}
1828
1829
1830
1831template <int dim, int spacedim>
1835{
1836 return cell_iterator(&this->get_triangulation(), -1, -1, this);
1837}
1838
1839
1840
1841template <int dim, int spacedim>
1844 DoFHandler<dim, spacedim>::end(const unsigned int level) const
1845{
1847 this->get_triangulation().end(level);
1848 if (cell.state() != IteratorState::valid)
1849 return end();
1850 return cell_iterator(*cell, this);
1851}
1852
1853
1854
1855template <int dim, int spacedim>
1858 DoFHandler<dim, spacedim>::end_active(const unsigned int level) const
1859{
1861 this->get_triangulation().end_active(level);
1862 if (cell.state() != IteratorState::valid)
1863 return active_cell_iterator(end());
1864 return active_cell_iterator(*cell, this);
1865}
1866
1867
1868
1869template <int dim, int spacedim>
1872 DoFHandler<dim, spacedim>::begin_mg(const unsigned int level) const
1873{
1874 Assert(this->has_level_dofs(),
1875 ExcMessage("You can only iterate over mg "
1876 "levels if mg dofs got distributed."));
1878 this->get_triangulation().begin(level);
1879 if (cell == this->get_triangulation().end(level))
1880 return end_mg(level);
1881 return level_cell_iterator(*cell, this);
1882}
1883
1884
1885
1886template <int dim, int spacedim>
1889 DoFHandler<dim, spacedim>::end_mg(const unsigned int level) const
1890{
1891 Assert(this->has_level_dofs(),
1892 ExcMessage("You can only iterate over mg "
1893 "levels if mg dofs got distributed."));
1895 this->get_triangulation().end(level);
1896 if (cell.state() != IteratorState::valid)
1897 return end();
1898 return level_cell_iterator(*cell, this);
1899}
1900
1901
1902
1903template <int dim, int spacedim>
1907{
1908 return level_cell_iterator(&this->get_triangulation(), -1, -1, this);
1909}
1910
1911
1912
1913template <int dim, int spacedim>
1916 dim,
1917 spacedim>::cell_iterators() const
1918{
1920 begin(), end());
1921}
1922
1923
1924
1925template <int dim, int spacedim>
1928 active_cell_iterator> DoFHandler<dim, spacedim>::
1930{
1931 return IteratorRange<
1932 typename DoFHandler<dim, spacedim>::active_cell_iterator>(begin_active(),
1933 end());
1934}
1935
1936
1937
1938template <int dim, int spacedim>
1941 typename DoFHandler<dim, spacedim>::
1942 level_cell_iterator> DoFHandler<dim, spacedim>::mg_cell_iterators() const
1943{
1945 begin_mg(), end_mg());
1946}
1947
1948
1949
1950template <int dim, int spacedim>
1953 dim,
1954 spacedim>::cell_iterators_on_level(const unsigned int level) const
1955{
1957 begin(level), end(level));
1958}
1959
1960
1961
1962template <int dim, int spacedim>
1965 active_cell_iterator> DoFHandler<dim, spacedim>::
1966 active_cell_iterators_on_level(const unsigned int level) const
1967{
1968 return IteratorRange<
1970 begin_active(level), end_active(level));
1971}
1972
1973
1974
1975template <int dim, int spacedim>
1978 level_cell_iterator> DoFHandler<dim, spacedim>::
1979 mg_cell_iterators_on_level(const unsigned int level) const
1980{
1982 begin_mg(level), end_mg(level));
1983}
1984
1985
1986
1987//---------------------------------------------------------------------------
1988
1989
1990
1991template <int dim, int spacedim>
1994{
1995 Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
1996 ExcNotImplementedWithHP());
1997
1998 Assert(this->fe_collection.size() > 0, ExcNoFESelected());
1999
2000 std::unordered_set<types::global_dof_index> boundary_dofs;
2001 std::vector<types::global_dof_index> dofs_on_face;
2002 dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2003
2004 const IndexSet &owned_dofs = locally_owned_dofs();
2005
2006 // loop over all faces to check whether they are at a
2007 // boundary. note that we need not take special care of single
2008 // lines in 3d (using @p{cell->has_boundary_lines}), since we do
2009 // not support boundaries of dimension dim-2, and so every
2010 // boundary line is also part of a boundary face.
2011 for (const auto &cell : this->active_cell_iterators())
2012 if (cell->is_locally_owned() && cell->at_boundary())
2013 {
2014 for (const auto iface : cell->face_indices())
2015 {
2016 const auto face = cell->face(iface);
2017 if (face->at_boundary())
2018 {
2019 const unsigned int dofs_per_face =
2020 cell->get_fe().n_dofs_per_face(iface);
2021 dofs_on_face.resize(dofs_per_face);
2022
2023 face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2024 for (unsigned int i = 0; i < dofs_per_face; ++i)
2025 {
2026 const unsigned int global_idof_index = dofs_on_face[i];
2027 if (owned_dofs.is_element(global_idof_index))
2028 {
2029 boundary_dofs.insert(global_idof_index);
2030 }
2031 }
2032 }
2033 }
2034 }
2035 return boundary_dofs.size();
2036}
2037
2038
2039
2040template <int dim, int spacedim>
2043 const std::set<types::boundary_id> &boundary_ids) const
2044{
2045 Assert(!(dim == 2 && spacedim == 3) || hp_capability_enabled == false,
2046 ExcNotImplementedWithHP());
2047
2048 Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2049 Assert(boundary_ids.find(numbers::internal_face_boundary_id) ==
2050 boundary_ids.end(),
2051 ExcInvalidBoundaryIndicator());
2052
2053 // same as above, but with additional checks for set of boundary
2054 // indicators
2055 std::unordered_set<types::global_dof_index> boundary_dofs;
2056 std::vector<types::global_dof_index> dofs_on_face;
2057 dofs_on_face.reserve(this->get_fe_collection().max_dofs_per_face());
2058
2059 const IndexSet &owned_dofs = locally_owned_dofs();
2060
2061 for (const auto &cell : this->active_cell_iterators())
2062 if (cell->is_locally_owned() && cell->at_boundary())
2063 {
2064 for (const auto iface : cell->face_indices())
2065 {
2066 const auto face = cell->face(iface);
2067 const unsigned int boundary_id = face->boundary_id();
2068 if (face->at_boundary() &&
2069 (boundary_ids.find(boundary_id) != boundary_ids.end()))
2070 {
2071 const unsigned int dofs_per_face =
2072 cell->get_fe().n_dofs_per_face(iface);
2073 dofs_on_face.resize(dofs_per_face);
2074
2075 face->get_dof_indices(dofs_on_face, cell->active_fe_index());
2076 for (unsigned int i = 0; i < dofs_per_face; ++i)
2077 {
2078 const unsigned int global_idof_index = dofs_on_face[i];
2079 if (owned_dofs.is_element(global_idof_index))
2080 {
2081 boundary_dofs.insert(global_idof_index);
2082 }
2083 }
2084 }
2085 }
2086 }
2087 return boundary_dofs.size();
2088}
2089
2090
2091
2092template <int dim, int spacedim>
2095{
2096 std::size_t mem = MemoryConsumption::memory_consumption(this->tria) +
2097 MemoryConsumption::memory_consumption(this->fe_collection) +
2098 MemoryConsumption::memory_consumption(this->number_cache);
2099
2100 mem += MemoryConsumption::memory_consumption(object_dof_indices) +
2102 MemoryConsumption::memory_consumption(hp_object_fe_indices) +
2103 MemoryConsumption::memory_consumption(hp_object_fe_ptr) +
2104 MemoryConsumption::memory_consumption(hp_cell_active_fe_indices) +
2105 MemoryConsumption::memory_consumption(hp_cell_future_fe_indices);
2106
2107
2108 if (hp_capability_enabled)
2109 {
2110 // nothing to add
2111 }
2112 else
2113 {
2114 // collect size of multigrid data structures
2115
2116 mem += MemoryConsumption::memory_consumption(this->block_info_object);
2117
2118 for (unsigned int level = 0; level < this->mg_levels.size(); ++level)
2119 mem += this->mg_levels[level]->memory_consumption();
2120
2121 if (this->mg_faces != nullptr)
2122 mem += MemoryConsumption::memory_consumption(*this->mg_faces);
2123
2124 for (unsigned int i = 0; i < this->mg_vertex_dofs.size(); ++i)
2125 mem += sizeof(MGVertexDoFs) +
2126 (1 + this->mg_vertex_dofs[i].get_finest_level() -
2127 this->mg_vertex_dofs[i].get_coarsest_level()) *
2129 }
2130
2131 return mem;
2132}
2133
2134
2135
2136template <int dim, int spacedim>
2140{
2141 this->distribute_dofs(hp::FECollection<dim, spacedim>(fe));
2142}
2143
2144
2145
2146template <int dim, int spacedim>
2150{
2151 Assert(this->tria != nullptr,
2152 ExcMessage(
2153 "You need to set the Triangulation in the DoFHandler using reinit() "
2154 "or in the constructor before you can distribute DoFs."));
2155 Assert(this->tria->n_levels() > 0,
2156 ExcMessage("The Triangulation you are using is empty!"));
2157
2158 // verify size of provided FE collection
2159 Assert(ff.size() > 0, ExcMessage("The given hp::FECollection is empty!"));
2160 Assert((ff.size() <= std::numeric_limits<types::fe_index>::max()) &&
2162 ExcMessage("The given hp::FECollection contains more finite elements "
2163 "than the DoFHandler can cover with active FE indices."));
2164
2165 if constexpr (running_in_debug_mode())
2166 {
2167 // make sure that the provided FE collection is large enough to
2168 // cover all FE indices presently in use on the mesh
2169 if ((hp_cell_active_fe_indices.size() > 0) &&
2170 (hp_cell_future_fe_indices.size() > 0))
2171 {
2172 Assert(hp_capability_enabled, ExcInternalError());
2173
2174 for (const auto &cell : this->active_cell_iterators() |
2175 IteratorFilters::LocallyOwnedCell())
2176 {
2177 Assert(cell->active_fe_index() < ff.size(),
2178 ExcInvalidFEIndex(cell->active_fe_index(), ff.size()));
2179 Assert(cell->future_fe_index() < ff.size(),
2180 ExcInvalidFEIndex(cell->future_fe_index(), ff.size()));
2181 }
2182 }
2183 }
2184
2185 //
2186 // register the new finite element collection
2187 //
2188 // don't create a new object if the one we have is identical
2189 if (this->fe_collection != ff)
2190 {
2191 this->fe_collection = hp::FECollection<dim, spacedim>(ff);
2192
2193 const bool contains_multiple_fes = (this->fe_collection.size() > 1);
2194
2195 // disable hp-mode if only a single finite element has been registered
2196 if (hp_capability_enabled && !contains_multiple_fes)
2197 {
2198 hp_capability_enabled = false;
2199
2200 // unsubscribe connections to signals that are only relevant for
2201 // hp-mode, since we only have a single element here
2202 for (auto &connection : this->tria_listeners_for_transfer)
2203 connection.disconnect();
2204 this->tria_listeners_for_transfer.clear();
2205
2206 // release active and future finite element tables
2207 this->hp_cell_active_fe_indices.clear();
2208 this->hp_cell_active_fe_indices.shrink_to_fit();
2209 this->hp_cell_future_fe_indices.clear();
2210 this->hp_cell_future_fe_indices.shrink_to_fit();
2211 }
2212
2213 // re-enabling hp-mode is not permitted since the active and future FE
2214 // tables are no longer available
2216 hp_capability_enabled || !contains_multiple_fes,
2217 ExcMessage(
2218 "You cannot re-enable hp-capabilities after you registered a single "
2219 "finite element. Please call reinit() or create a new DoFHandler "
2220 "object instead."));
2221 }
2222
2223 //
2224 // enumerate all degrees of freedom
2225 //
2226 if (hp_capability_enabled)
2227 {
2228 // make sure every processor knows the active FE indices
2229 // on both its own cells and all ghost cells
2232 }
2233
2234 {
2235 // We would like to enumerate all dofs for shared::Triangulations. If an
2236 // underlying shared::Tria allows artificial cells, we need to restore the
2237 // true cell owners temporarily.
2238 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
2239 // the current set of subdomain ids, set subdomain ids to the "true" owner
2240 // of each cell upon construction of the TemporarilyRestoreSubdomainIds
2241 // object, and later restore these flags when it is destroyed.
2243 spacedim>
2244 subdomain_modifier(this->get_triangulation());
2245
2246 // Adjust size of levels to the triangulation. Note that we still have to
2247 // allocate space for all degrees of freedom on this mesh (including ghost
2248 // and cells that are entirely stored on different processors), though we
2249 // may not assign numbers to some of them (i.e. they will remain at
2250 // invalid_dof_index). We need to allocate the space because we will want
2251 // to be able to query the dof_indices on each cell, and simply be told
2252 // that we don't know them on some cell (i.e. get back invalid_dof_index)
2253 if (hp_capability_enabled)
2255 *this);
2256 else
2258 }
2259
2260 // hand the actual work over to the policy
2261 this->number_cache = this->policy->distribute_dofs();
2262
2263 // do some housekeeping: compress indices
2264 // if(hp_capability_enabled)
2265 // {
2266 // Threads::TaskGroup<> tg;
2267 // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2268 // tg += Threads::new_task(
2269 // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2270 // *this->levels_hp[level],
2271 // this->fe_collection);
2272 // tg.join_all();
2273 // }
2274
2275 // Initialize the block info object only if this is a sequential
2276 // triangulation. It doesn't work correctly yet if it is parallel and has not
2277 // yet been implemented for hp-mode.
2278 if (!hp_capability_enabled &&
2280 *>(&*this->tria) == nullptr)
2281 this->block_info_object.initialize(*this, false, true);
2282}
2283
2284
2285
2286template <int dim, int spacedim>
2289{
2290 AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2291
2292 Assert(
2293 this->object_dof_indices.size() > 0,
2294 ExcMessage(
2295 "Distribute active DoFs using distribute_dofs() before calling distribute_mg_dofs()."));
2296
2297 Assert(
2298 ((this->tria->get_mesh_smoothing() &
2301 ExcMessage(
2302 "The mesh smoothing requirement 'limit_level_difference_at_vertices' has to be set for using multigrid!"));
2303
2304 this->clear_mg_space();
2305
2307 this->mg_number_cache = this->policy->distribute_mg_dofs();
2308
2309 // initialize the block info object only if this is a sequential
2310 // triangulation. it doesn't work correctly yet if it is parallel
2311 if (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2312 &*this->tria) == nullptr)
2313 this->block_info_object.initialize(*this, true, false);
2314}
2315
2316
2317
2318template <int dim, int spacedim>
2321{
2322 AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2323
2324 this->block_info_object.initialize_local(*this);
2325}
2326
2327
2328
2329template <int dim, int spacedim>
2332{
2333 // decide whether we need a sequential or a parallel distributed policy
2334 if (dynamic_cast<const ::parallel::shared::Triangulation<dim, spacedim>
2335 *>(&this->get_triangulation()) != nullptr)
2336 this->policy = std::make_unique<internal::DoFHandlerImplementation::Policy::
2337 ParallelShared<dim, spacedim>>(*this);
2338 else if (dynamic_cast<
2339 const ::parallel::DistributedTriangulationBase<dim, spacedim>
2340 *>(&this->get_triangulation()) == nullptr)
2341 this->policy = std::make_unique<
2343 *this);
2344 else
2345 this->policy =
2346 std::make_unique<internal::DoFHandlerImplementation::Policy::
2347 ParallelDistributed<dim, spacedim>>(*this);
2348}
2349
2350
2351
2352template <int dim, int spacedim>
2355{
2356 // release memory
2357 this->clear_space();
2358 this->clear_mg_space();
2359}
2360
2361
2362
2363template <int dim, int spacedim>
2366{
2367 object_dof_indices.clear();
2368
2369 object_dof_ptr.clear();
2370
2371 this->number_cache.clear();
2372
2373 this->hp_cell_active_fe_indices.clear();
2374 this->hp_cell_future_fe_indices.clear();
2375}
2376
2377
2378
2379template <int dim, int spacedim>
2382{
2383 this->mg_levels.clear();
2384 this->mg_faces.reset();
2385
2386 std::vector<MGVertexDoFs> tmp;
2387
2388 std::swap(this->mg_vertex_dofs, tmp);
2389
2390 this->mg_number_cache.clear();
2391}
2392
2393
2394
2395template <int dim, int spacedim>
2398 const std::vector<types::global_dof_index> &new_numbers)
2399{
2400 if (hp_capability_enabled)
2401 {
2402 Assert(this->hp_cell_future_fe_indices.size() > 0,
2403 ExcMessage(
2404 "You need to distribute DoFs before you can renumber them."));
2405
2406 AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2407
2408 if constexpr (running_in_debug_mode())
2409 {
2410 // assert that the new indices are consecutively numbered if we are
2411 // working on a single processor. this doesn't need to
2412 // hold in the case of a parallel mesh since we map the interval
2413 // [0...n_dofs()) into itself but only globally, not on each processor
2414 if (this->n_locally_owned_dofs() == this->n_dofs())
2415 {
2416 std::vector<types::global_dof_index> tmp(new_numbers);
2417 std::sort(tmp.begin(), tmp.end());
2418 std::vector<types::global_dof_index>::const_iterator p =
2419 tmp.begin();
2421 for (; p != tmp.end(); ++p, ++i)
2422 Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2423 }
2424 else
2425 for (const auto new_number : new_numbers)
2426 Assert(
2427 new_number < this->n_dofs(),
2428 ExcMessage(
2429 "New DoF index is not less than the total number of dofs."));
2430 }
2431
2432 // uncompress the internal storage scheme of dofs on cells so that
2433 // we can access dofs in turns. uncompress in parallel, starting
2434 // with the most expensive levels (the highest ones)
2435 //{
2436 // Threads::TaskGroup<> tg;
2437 // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2438 // tg += Threads::new_task(
2439 // &::internal::hp::DoFLevel::uncompress_data<dim, spacedim>,
2440 // *this->levels_hp[level],
2441 // this->fe_collection);
2442 // tg.join_all();
2443 //}
2444
2445 // do the renumbering
2446 this->number_cache = this->policy->renumber_dofs(new_numbers);
2447
2448 // now re-compress the dof indices
2449 //{
2450 // Threads::TaskGroup<> tg;
2451 // for (int level = this->levels_hp.size() - 1; level >= 0; --level)
2452 // tg += Threads::new_task(
2453 // &::internal::hp::DoFLevel::compress_data<dim, spacedim>,
2454 // *this->levels_hp[level],
2455 // this->fe_collection);
2456 // tg.join_all();
2457 //}
2458 }
2459 else
2460 {
2461 Assert(this->object_dof_indices.size() > 0,
2462 ExcMessage(
2463 "You need to distribute DoFs before you can renumber them."));
2464
2465 if constexpr (running_in_debug_mode())
2466 {
2468 *>(&*this->tria) != nullptr)
2469 {
2470 Assert(new_numbers.size() == this->n_dofs() ||
2471 new_numbers.size() == this->n_locally_owned_dofs(),
2472 ExcMessage("Incorrect size of the input array."));
2473 }
2474 else if (dynamic_cast<
2476 *>(&*this->tria) != nullptr)
2477 {
2478 AssertDimension(new_numbers.size(), this->n_locally_owned_dofs());
2479 }
2480 else
2481 {
2482 AssertDimension(new_numbers.size(), this->n_dofs());
2483 }
2484
2485 // assert that the new indices are consecutively numbered if we are
2486 // working on a single processor. this doesn't need to
2487 // hold in the case of a parallel mesh since we map the interval
2488 // [0...n_dofs()) into itself but only globally, not on each processor
2489 if (this->n_locally_owned_dofs() == this->n_dofs())
2490 {
2491 std::vector<types::global_dof_index> tmp(new_numbers);
2492 std::sort(tmp.begin(), tmp.end());
2493 std::vector<types::global_dof_index>::const_iterator p =
2494 tmp.begin();
2496 for (; p != tmp.end(); ++p, ++i)
2497 Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2498 }
2499 else
2500 for (const auto new_number : new_numbers)
2501 Assert(
2502 new_number < this->n_dofs(),
2503 ExcMessage(
2504 "New DoF index is not less than the total number of dofs."));
2505 }
2506
2507 this->number_cache = this->policy->renumber_dofs(new_numbers);
2508 }
2509}
2510
2511
2512
2513template <int dim, int spacedim>
2516 const unsigned int level,
2517 const std::vector<types::global_dof_index> &new_numbers)
2518{
2519 AssertThrow(hp_capability_enabled == false, ExcNotImplementedWithHP());
2520
2521 Assert(
2522 this->mg_levels.size() > 0 && this->object_dof_indices.size() > 0,
2523 ExcMessage(
2524 "You need to distribute active and level DoFs before you can renumber level DoFs."));
2525 AssertIndexRange(level, this->get_triangulation().n_global_levels());
2526 AssertDimension(new_numbers.size(),
2527 this->locally_owned_mg_dofs(level).n_elements());
2528
2529 if constexpr (running_in_debug_mode())
2530 {
2531 // assert that the new indices are consecutively numbered if we are
2532 // working on a single processor. this doesn't need to hold in the case of
2533 // a parallel mesh since we map the interval [0...n_dofs(level)) into
2534 // itself but only globally, not on each processor
2535 if (this->n_locally_owned_dofs() == this->n_dofs())
2536 {
2537 std::vector<types::global_dof_index> tmp(new_numbers);
2538 std::sort(tmp.begin(), tmp.end());
2539 std::vector<types::global_dof_index>::const_iterator p = tmp.begin();
2541 for (; p != tmp.end(); ++p, ++i)
2542 Assert(*p == i, ExcNewNumbersNotConsecutive(i));
2543 }
2544 else
2545 for (const auto new_number : new_numbers)
2546 Assert(new_number < this->n_dofs(level),
2547 ExcMessage(
2548 "New DoF index is not less than the total number of dofs."));
2549 }
2550
2551 this->mg_number_cache[level] =
2552 this->policy->renumber_mg_dofs(level, new_numbers);
2553}
2554
2555
2556
2557template <int dim, int spacedim>
2560 const
2561{
2562 Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2563
2564 switch (dim)
2565 {
2566 case 1:
2567 return this->fe_collection.max_dofs_per_vertex();
2568 case 2:
2569 return (3 * this->fe_collection.max_dofs_per_vertex() +
2570 2 * this->fe_collection.max_dofs_per_line());
2571 case 3:
2572 // we need to take refinement of one boundary face into
2573 // consideration here; in fact, this function returns what
2574 // #max_coupling_between_dofs<2> returns
2575 //
2576 // we assume here, that only four faces meet at the boundary;
2577 // this assumption is not justified and needs to be fixed some
2578 // time. fortunately, omitting it for now does no harm since
2579 // the matrix will cry foul if its requirements are not
2580 // satisfied
2581 return (19 * this->fe_collection.max_dofs_per_vertex() +
2582 28 * this->fe_collection.max_dofs_per_line() +
2583 8 * this->fe_collection.max_dofs_per_quad());
2584 default:
2586 return 0;
2587 }
2588}
2589
2590
2591
2592template <int dim, int spacedim>
2595{
2596 Assert(this->fe_collection.size() > 0, ExcNoFESelected());
2599}
2600
2601
2602
2603template <int dim, int spacedim>
2606 const std::vector<types::fe_index> &active_fe_indices)
2607{
2608 Assert(active_fe_indices.size() == this->get_triangulation().n_active_cells(),
2609 ExcDimensionMismatch(active_fe_indices.size(),
2610 this->get_triangulation().n_active_cells()));
2611
2612 this->create_active_fe_table();
2613 // we could set the values directly, since they are stored as
2614 // protected data of this object, but for simplicity we use the
2615 // cell-wise access. this way we also have to pass some debug-mode
2616 // tests which we would have to duplicate ourselves otherwise
2617 for (const auto &cell : this->active_cell_iterators())
2618 if (cell->is_locally_owned())
2619 cell->set_active_fe_index(active_fe_indices[cell->active_cell_index()]);
2620}
2621
2622
2623
2624template <int dim, int spacedim>
2627 const std::vector<unsigned int> &active_fe_indices)
2628{
2629 set_active_fe_indices(std::vector<types::fe_index>(active_fe_indices.begin(),
2630 active_fe_indices.end()));
2631}
2632
2633
2634
2635template <int dim, int spacedim>
2637std::vector<types::fe_index> DoFHandler<dim, spacedim>::get_active_fe_indices()
2638 const
2639{
2640 std::vector<types::fe_index> active_fe_indices(
2641 this->get_triangulation().n_active_cells(), numbers::invalid_fe_index);
2642
2643 // we could try to extract the values directly, since they are
2644 // stored as protected data of this object, but for simplicity we
2645 // use the cell-wise access.
2646 for (const auto &cell : this->active_cell_iterators())
2647 if (!cell->is_artificial())
2648 active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
2649
2650 return active_fe_indices;
2651}
2652
2653
2654
2655template <int dim, int spacedim>
2658 std::vector<unsigned int> &active_fe_indices) const
2659{
2660 const std::vector<types::fe_index> indices = get_active_fe_indices();
2661
2662 active_fe_indices.assign(indices.begin(), indices.end());
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
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()
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:1905
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:2527
unsigned int n_vertices() const
unsigned int size() const
Definition collection.h:316
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:68
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:245
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int level
Definition grid_out.cc:4635
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)
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:745
void get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< unsigned int > &active_fe_indices)
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)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
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:466
Definition hp.h:117
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14909
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:263
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
constexpr types::boundary_id internal_face_boundary_id
Definition types.h:323
constexpr types::fe_index invalid_fe_index
Definition types.h:254
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
Definition types.h:68
unsigned int boundary_id
Definition types.h:153
boost::signals2::signal< void()> post_distributed_load
Definition tria.h:2516
boost::signals2::signal< void()> post_distributed_refinement
Definition tria.h:2472
boost::signals2::signal< void()> pre_refinement
Definition tria.h:2318
boost::signals2::signal< void()> create
Definition tria.h:2309
boost::signals2::signal< void()> clear
Definition tria.h:2388
boost::signals2::signal< void()> post_refinement
Definition tria.h:2325
boost::signals2::signal< void()> pre_partition
Definition tria.h:2333
boost::signals2::signal< void()> pre_distributed_repartition
Definition tria.h:2479
boost::signals2::signal< void()> post_p4est_refinement
Definition tria.h:2462
boost::signals2::signal< void()> post_distributed_repartition
Definition tria.h:2486
boost::signals2::signal< void()> post_distributed_save
Definition tria.h:2501
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)