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