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