deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
dof_handler_policy.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2010 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13
14#include <deal.II/base/geometry_info.h>
15#include <deal.II/base/memory_consumption.h>
16#include <deal.II/base/partitioner.h>
17#include <deal.II/base/thread_management.h>
18#include <deal.II/base/types.h>
19#include <deal.II/base/utilities.h>
20#include <deal.II/base/work_stream.h>
21
22#include <deal.II/distributed/shared_tria.h>
23#include <deal.II/distributed/tria.h>
24
25#include <deal.II/dofs/dof_accessor.h>
26#include <deal.II/dofs/dof_handler.h>
27#include <deal.II/dofs/dof_handler_policy.h>
28
29#include <deal.II/fe/fe.h>
30
31#include <deal.II/grid/grid_tools.h>
32#include <deal.II/grid/tria.h>
33#include <deal.II/grid/tria_iterator.h>
34
35#include <algorithm>
36#include <limits>
37#include <memory>
38#include <numeric>
39#include <set>
40
41DEAL_II_NAMESPACE_OPEN
42
43
44namespace internal
45{
46 namespace DoFHandlerImplementation
47 {
48 namespace Policy
49 {
50 namespace
51 {
58 const types::global_dof_index enumeration_dof_index =
59 numbers::invalid_dof_index - 1;
60
61
62 using DoFIdentities =
63 std::vector<std::pair<unsigned int, unsigned int>>;
64
65
76 template <int structdim, int dim, int spacedim>
77 const std::unique_ptr<DoFIdentities> &
78 ensure_existence_and_return_dof_identities(
79 const ::hp::FECollection<dim, spacedim> &fes,
80 const types::fe_index fe_index_1,
81 const types::fe_index fe_index_2,
82 std::unique_ptr<DoFIdentities> &identities,
83 const unsigned int face_no = numbers::invalid_unsigned_int,
84 const unsigned int face_no_neighbor = numbers::invalid_unsigned_int)
85 {
86 Assert(structdim == 2 ||
87 (face_no == numbers::invalid_unsigned_int &&
88 face_no_neighbor == numbers::invalid_unsigned_int),
89 ExcInternalError());
90 Assert((face_no == numbers::invalid_unsigned_int &&
91 face_no_neighbor == numbers::invalid_unsigned_int) ||
92 (face_no != numbers::invalid_unsigned_int &&
93 face_no_neighbor != numbers::invalid_unsigned_int),
94 ExcInternalError());
95
96 // see if we need to fill this entry, or whether it already
97 // exists
98 if (identities.get() == nullptr)
99 {
100 // TODO: Change to
101 // std::vector<std::map<types::fe_index, unsigned int>>
102 std::vector<std::map<unsigned int, unsigned int>>
103 complete_identities;
104
105 switch (structdim)
106 {
107 case 0:
108 {
109 // TODO: Change set to types::fe_index
110 complete_identities = fes.hp_vertex_dof_identities(
111 std::set<unsigned int>{fe_index_1, fe_index_2});
112 break;
113 }
114
115 case 1:
116 {
117 // TODO: Change set to types::fe_index
118 complete_identities = fes.hp_line_dof_identities(
119 std::set<unsigned int>{fe_index_1, fe_index_2});
120 break;
121 }
122
123 case 2:
124 {
125 // TODO: Change set to types::fe_index
126 std::pair<unsigned int, unsigned int> p1{fe_index_1,
127 face_no};
128 std::pair<unsigned int, unsigned int> p2{
129 fe_index_2, face_no_neighbor};
130
131 complete_identities = fes.hp_quad_dof_identities(
132 std::set<std::pair<unsigned int, unsigned int>>{p1,
133 p2});
134 break;
135 }
136
137 default:
138 DEAL_II_NOT_IMPLEMENTED();
139 }
140
141 if constexpr (running_in_debug_mode())
142 {
143 // Each entry of 'complete_identities' contains a set of
144 // pairs (fe_index,dof_index). Because we put in exactly
145 // two fe indices, we know that each entry of the outer
146 // vector needs to contain a set of exactly two such
147 // pairs. Check this. While there, also check that
148 // the two entries actually reference fe_index_1 and
149 // fe_index_2:
150 for (const auto &complete_identity : complete_identities)
151 {
152 Assert(complete_identity.size() == 2, ExcInternalError());
153 Assert(complete_identity.find(fe_index_1) !=
154 complete_identity.end(),
155 ExcInternalError());
156 Assert(complete_identity.find(fe_index_2) !=
157 complete_identity.end(),
158 ExcInternalError());
159 }
160 }
161
162 // Next reduce these sets of two pairs by removing the
163 // fe_index parts: We know which indices we have. But we
164 // have to make sure in which order we consider the
165 // pair, by considering whether the fe_index part we are
166 // throwing away matched fe_index_1 or fe_index_2. Fortunately,
167 // this is easy to do because we can ask the std::map for the
168 // dof_index that matches a given fe_index:
169 DoFIdentities reduced_identities;
170 for (const auto &complete_identity : complete_identities)
171 {
172 const unsigned int dof_index_1 =
173 complete_identity.at(fe_index_1);
174 const unsigned int dof_index_2 =
175 complete_identity.at(fe_index_2);
176
177 reduced_identities.emplace_back(dof_index_1, dof_index_2);
178 }
179
180 if constexpr (running_in_debug_mode())
181 {
182 // double check whether the newly created entries make
183 // any sense at all
184 for (const auto &identity : reduced_identities)
185 {
186 Assert(
187 identity.first <
188 fes[fe_index_1].template n_dofs_per_object<structdim>(
189 face_no),
190 ExcInternalError());
191 Assert(
192 identity.second <
193 fes[fe_index_2].template n_dofs_per_object<structdim>(
194 face_no_neighbor),
195 ExcInternalError());
196 }
197 }
198
199 identities =
200 std::make_unique<DoFIdentities>(std::move(reduced_identities));
201 }
202
203 return identities;
204 }
205 } // namespace
206
207
208
209 struct Implementation
210 {
211 /* -------------- distribute_dofs functionality ------------- */
212
217 template <int dim, int spacedim>
218 static std::map<types::global_dof_index, types::global_dof_index>
219 compute_vertex_dof_identities(
220 const DoFHandler<dim, spacedim> &dof_handler)
221 {
222 Assert(
223 dof_handler.hp_capability_enabled == true,
224 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
225
226 std::map<types::global_dof_index, types::global_dof_index>
227 dof_identities;
228
229 // Note: we may wish to have something here similar to what
230 // we do for lines and quads, namely that we only identify
231 // dofs for any FE towards the most dominating one. however,
232 // it is not clear whether this is actually necessary for
233 // vertices at all, I can't think of a finite element that
234 // would make that necessary...
235 ::Table<2, std::unique_ptr<DoFIdentities>>
236 vertex_dof_identities(dof_handler.get_fe_collection().size(),
237 dof_handler.get_fe_collection().size());
238
239 // loop over all vertices and see which one we need to work on
240 for (unsigned int vertex_index = 0;
241 vertex_index < dof_handler.get_triangulation().n_vertices();
242 ++vertex_index)
243 if (dof_handler.get_triangulation()
244 .get_used_vertices()[vertex_index] == true)
245 {
246 const unsigned int n_active_fe_indices =
247 ::internal::DoFAccessorImplementation::Implementation::
248 n_active_fe_indices(dof_handler,
249 0,
250 vertex_index,
251 std::integral_constant<int, 0>());
252
253 if (n_active_fe_indices > 1)
254 {
255 const std::set<types::fe_index> fe_indices =
256 ::internal::DoFAccessorImplementation::
257 Implementation::get_active_fe_indices(
258 dof_handler,
259 0,
260 vertex_index,
261 std::integral_constant<int, 0>());
262
263 // find out which is the most dominating finite
264 // element of the ones that are used on this vertex
265 // TODO: Change set to types::fe_index
266 types::fe_index most_dominating_fe_index =
267 dof_handler.get_fe_collection().find_dominating_fe(
268 {fe_indices.begin(), fe_indices.end()},
269 /*codim*/ dim);
270
271 // if we haven't found a dominating finite element,
272 // choose the very first one to be dominant
273 // TODO: Change assert to numbers::invalid_fe_index
274 if (most_dominating_fe_index == numbers::invalid_fe_index)
275 most_dominating_fe_index =
276 ::internal::DoFAccessorImplementation::
277 Implementation::nth_active_fe_index(
278 dof_handler,
279 0,
280 vertex_index,
281 0,
282 std::integral_constant<int, 0>());
283
284 // loop over the indices of all the finite
285 // elements that are not dominating, and
286 // identify their dofs to the most dominating
287 // one
288 for (const auto &other_fe_index : fe_indices)
289 if (other_fe_index != most_dominating_fe_index)
290 {
291 // make sure the entry in the equivalence
292 // table exists
293 const auto &identities =
294 *ensure_existence_and_return_dof_identities<0>(
295 dof_handler.get_fe_collection(),
296 most_dominating_fe_index,
297 other_fe_index,
298 vertex_dof_identities[most_dominating_fe_index]
299 [other_fe_index]);
300
301 // then loop through the identities we
302 // have. first get the global numbers of the
303 // dofs we want to identify and make sure they
304 // are not yet constrained to anything else,
305 // except for to each other. use the rule that
306 // we will always constrain the dof with the
307 // higher FE index to the one with the lower,
308 // to avoid circular reasoning.
309 for (const auto &identity : identities)
310 {
311 const types::global_dof_index primary_dof_index =
312 ::internal::DoFAccessorImplementation::
313 Implementation::get_dof_index(
314 dof_handler,
315 0,
316 vertex_index,
317 most_dominating_fe_index,
318 identity.first,
319 std::integral_constant<int, 0>());
320 const types::global_dof_index
321 dependent_dof_index =
322 ::internal::DoFAccessorImplementation::
323 Implementation::get_dof_index(
324 dof_handler,
325 0,
326 vertex_index,
327 other_fe_index,
328 identity.second,
329 std::integral_constant<int, 0>());
330
331 // on subdomain boundaries, we will
332 // encounter invalid DoFs on ghost cells,
333 // for which we have not yet distributed
334 // valid indices. depending on which finte
335 // element is dominating the other on this
336 // interface, we either have to constrain
337 // the valid to the invalid indices, or vice
338 // versa.
339 //
340 // we only store an identity if we are about
341 // to overwrite a valid DoF. we will skip
342 // constraining invalid DoFs for now, and
343 // consider them later in Phase 5.
344 if (dependent_dof_index !=
345 numbers::invalid_dof_index)
346 {
347 // if the DoF indices of both elements
348 // are already distributed, i.e., both
349 // of these 'fe_indices' are associated
350 // with a locally owned cell, then we
351 // should either not have a dof_identity
352 // yet, or it must come out here to be
353 // exactly as we had computed before
354 if (primary_dof_index !=
355 numbers::invalid_dof_index)
356 Assert(
357 (dof_identities.find(primary_dof_index) ==
358 dof_identities.end()) ||
359 (dof_identities[dependent_dof_index] ==
360 primary_dof_index),
361 ExcInternalError());
362
363 dof_identities[dependent_dof_index] =
364 primary_dof_index;
365 }
366 }
367 }
368 }
369 }
370
371 return dof_identities;
372 }
373
374
379 template <int spacedim>
380 static std::map<types::global_dof_index, types::global_dof_index>
381 compute_line_dof_identities(const DoFHandler<1, spacedim> &dof_handler)
382 {
383 (void)dof_handler;
384 Assert(dof_handler.hp_capability_enabled == true,
385 (typename DoFHandler<1, spacedim>::ExcOnlyAvailableWithHP()));
386
387 return std::map<types::global_dof_index, types::global_dof_index>();
388 }
389
390
391 template <int dim, int spacedim>
392 static std::map<types::global_dof_index, types::global_dof_index>
393 compute_line_dof_identities(
394 const DoFHandler<dim, spacedim> &dof_handler)
395 {
396 Assert(
397 dof_handler.hp_capability_enabled == true,
398 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
399
400 std::map<types::global_dof_index, types::global_dof_index>
401 dof_identities;
402
403 // An implementation of the algorithm described in the hp-paper,
404 // including the modification mentioned later in the "complications in
405 // 3-d" subsections
406 //
407 // as explained there, we do something only if there are exactly 2
408 // finite elements associated with an object. if there is only one,
409 // then there is nothing to do anyway, and if there are 3 or more,
410 // then we can get into trouble. note that this only happens for lines
411 // in 3d and higher, and for quads only in 4d and higher, so this
412 // isn't a particularly frequent case
413 //
414 // there is one case, however, that we would like to handle (see, for
415 // example, the hp/crash_15 testcase): if we have
416 // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
417 // then we should be able to handle this because we can simply unify
418 // *all* dofs, not only a some. so what we do is to first treat all
419 // pairs of finite elements that have *identical* dofs, and then only
420 // deal with those that are not identical of which we can handle at
421 // most 2
422 ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
423 dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
424
425 std::vector<bool> line_touched(
426 dof_handler.get_triangulation().n_raw_lines());
427 for (const auto &cell : dof_handler.active_cell_iterators())
428 for (const auto l : cell->line_indices())
429 if (!line_touched[cell->line(l)->index()])
430 {
431 const auto line = cell->line(l);
432 line_touched[line->index()] = true;
433
434 unsigned int unique_sets_of_dofs =
435 line->n_active_fe_indices();
436
437 // do a first loop over all sets of dofs and do identity
438 // uniquification
439 const unsigned int n_active_fe_indices =
440 line->n_active_fe_indices();
441 for (unsigned int f = 0; f < n_active_fe_indices; ++f)
442 for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
443 {
444 const types::fe_index fe_index_1 =
445 line->nth_active_fe_index(f),
446 fe_index_2 =
447 line->nth_active_fe_index(g);
448
449 // as described in the hp-paper, we only unify on lines
450 // when there are at most two different FE objects
451 // assigned on it.
452 // however, more than two 'active_fe_indices' can be
453 // attached that still fulfill the above criterion,
454 // i.e. when two different FiniteElement objects are
455 // assigned to neighboring cells that map their degrees
456 // of freedom one-to-one.
457 // we cannot verify with certainty if two dofs each of
458 // separate FiniteElement objects actually map
459 // one-to-one. however, checking for the number of
460 // 'dofs_per_line' turned out to be a reasonable
461 // approach, that also works for e.g. two different
462 // FE_Q objects of the same order, from which one is
463 // enhanced by a bubble function that is zero on the
464 // boundary.
465 if ((dof_handler.get_fe(fe_index_1).n_dofs_per_line() ==
466 dof_handler.get_fe(fe_index_2)
467 .n_dofs_per_line()) &&
468 (dof_handler.get_fe(fe_index_1).n_dofs_per_line() >
469 0))
470 {
471 // the number of dofs per line is identical
472 const unsigned int dofs_per_line =
473 dof_handler.get_fe(fe_index_1).n_dofs_per_line();
474
475 const auto &identities =
476 *ensure_existence_and_return_dof_identities<1>(
477 dof_handler.get_fe_collection(),
478 fe_index_1,
479 fe_index_2,
480 line_dof_identities[fe_index_1][fe_index_2]);
481 // see if these sets of dofs are identical. the
482 // first condition for this is that indeed there are
483 // n identities
484 if (identities.size() == dofs_per_line)
485 {
486 unsigned int i = 0;
487 for (; i < dofs_per_line; ++i)
488 if (identities[i] != std::pair{i, i})
489 // not an identity
490 break;
491
492 if (i == dofs_per_line)
493 {
494 // The line dofs (i.e., the ones interior to
495 // a line) of these two finite elements are
496 // identical. Note that there could be
497 // situations when one element still
498 // dominates another, e.g.: FE_Q(2) x
499 // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
500
501 --unique_sets_of_dofs;
502
503 // determine which one of both finite
504 // elements is the dominating one.
505 const std::set<types::fe_index> fe_indices{
506 fe_index_1, fe_index_2};
507
508 // TODO: Change set to types::fe_index
509 types::fe_index dominating_fe_index =
510 dof_handler.get_fe_collection()
511 .find_dominating_fe({fe_indices.begin(),
512 fe_indices.end()},
513 /*codim=*/dim - 1);
514 types::fe_index other_fe_index =
515 numbers::invalid_fe_index;
516
517 if (dominating_fe_index !=
518 numbers::invalid_fe_index)
519 other_fe_index =
520 (dominating_fe_index == fe_index_1) ?
521 fe_index_2 :
522 fe_index_1;
523 else
524 {
525 // if we haven't found a dominating
526 // finite element, choose the one with
527 // the lower index to be dominating
528 dominating_fe_index = fe_index_1;
529 other_fe_index = fe_index_2;
530 }
531
532 for (unsigned int j = 0; j < dofs_per_line;
533 ++j)
534 {
535 const types::global_dof_index
536 primary_dof_index = line->dof_index(
537 j, dominating_fe_index);
538 const types::global_dof_index
539 dependent_dof_index =
540 line->dof_index(j, other_fe_index);
541
542 // on subdomain boundaries, we will
543 // encounter invalid DoFs on ghost
544 // cells, for which we have not yet
545 // distributed valid indices. depending
546 // on which finte element is dominating
547 // the other on this interface, we
548 // either have to constrain the valid to
549 // the invalid indices, or vice versa.
550 //
551 // we only store an identity if we are
552 // about to overwrite a valid DoF. we
553 // will skip constraining invalid DoFs
554 // for now, and consider them later in
555 // Phase 5.
556 if (dependent_dof_index !=
557 numbers::invalid_dof_index)
558 {
559 if (primary_dof_index !=
560 numbers::invalid_dof_index)
561 {
562 // if primary dof was already
563 // constrained, constrain to
564 // that one, otherwise constrain
565 // dependent to primary
566 if (dof_identities.find(
567 primary_dof_index) !=
568 dof_identities.end())
569 {
570 // if the DoF indices of
571 // both elements are already
572 // distributed, i.e., both
573 // of these 'fe_indices' are
574 // associated with a locally
575 // owned cell, then we
576 // should either not have a
577 // dof_identity yet, or it
578 // must come out here to be
579 // exactly as we had
580 // computed before
581 Assert(
582 dof_identities.find(
583 dof_identities
584 [primary_dof_index]) ==
585 dof_identities.end(),
586 ExcInternalError());
587
588 dof_identities
589 [dependent_dof_index] =
590 dof_identities
591 [primary_dof_index];
592 }
593 else
594 {
595 // see comment above for an
596 // explanation of this
597 // assertion
598 Assert(
599 (dof_identities.find(
600 primary_dof_index) ==
601 dof_identities.end()) ||
602 (dof_identities
603 [dependent_dof_index] ==
604 primary_dof_index),
605 ExcInternalError());
606
607 dof_identities
608 [dependent_dof_index] =
609 primary_dof_index;
610 }
611 }
612 else
613 {
614 // set dependent_dof to
615 // primary_dof_index, which is
616 // invalid
617 dof_identities
618 [dependent_dof_index] =
619 numbers::invalid_dof_index;
620 }
621 }
622 }
623 }
624 }
625 }
626 }
627
628 // if at this point, there is only one unique set of dofs
629 // left, then we have taken care of everything above. if there
630 // are two, then we need to deal with them here. if there are
631 // more, then we punt, as described in the paper (and
632 // mentioned above)
633 // TODO: The check for 'dim==2' was inserted by intuition. It
634 // fixes
635 // the previous problems with @ref step_27 "step-27" in 3d. But an
636 // explanation for this is still required, and what we do here
637 // is not what we describe in the paper!.
638 if ((unique_sets_of_dofs == 2) && (dim == 2))
639 {
640 const std::set<types::fe_index> fe_indices =
641 line->get_active_fe_indices();
642
643 // find out which is the most dominating finite element of
644 // the ones that are used on this line
645 // TODO: Change set to types::fe_index
646 const types::fe_index most_dominating_fe_index =
647 dof_handler.get_fe_collection().find_dominating_fe(
648 {fe_indices.begin(), fe_indices.end()},
649 /*codim=*/dim - 1);
650
651 // if we found the most dominating element, then use this
652 // to eliminate some of the degrees of freedom by
653 // identification. otherwise, the code that computes
654 // hanging node constraints will have to deal with it by
655 // computing appropriate constraints along this face/edge
656 if (most_dominating_fe_index != numbers::invalid_fe_index)
657 {
658 // loop over the indices of all the finite elements
659 // that are not dominating, and identify their dofs to
660 // the most dominating one
661 for (const auto &other_fe_index : fe_indices)
662 if (other_fe_index != most_dominating_fe_index)
663 {
664 const auto &identities =
665 *ensure_existence_and_return_dof_identities<
666 1>(dof_handler.get_fe_collection(),
667 most_dominating_fe_index,
668 other_fe_index,
669 line_dof_identities
670 [most_dominating_fe_index]
671 [other_fe_index]);
672
673 for (const auto &identity : identities)
674 {
675 const types::global_dof_index
676 primary_dof_index = line->dof_index(
677 identity.first,
678 most_dominating_fe_index);
679 const types::global_dof_index
680 dependent_dof_index =
681 line->dof_index(identity.second,
682 other_fe_index);
683
684 // on subdomain boundaries, we will
685 // encounter invalid DoFs on ghost cells,
686 // for which we have not yet distributed
687 // valid indices. depending on which finte
688 // element is dominating the other on this
689 // interface, we either have to constrain
690 // the valid to the invalid indices, or vice
691 // versa.
692 //
693 // we only store an identity if we are about
694 // to overwrite a valid DoF. we will skip
695 // constraining invalid DoFs for now, and
696 // consider them later in Phase 5.
697 if (dependent_dof_index !=
698 numbers::invalid_dof_index)
699 {
700 // if the DoF indices of both elements
701 // are already distributed, i.e., both
702 // of these 'fe_indices' are associated
703 // with a locally owned cell, then we
704 // should either not have a dof_identity
705 // yet, or it must come out here to be
706 // exactly as we had computed before
707 if (primary_dof_index !=
708 numbers::invalid_dof_index)
709 Assert((dof_identities.find(
710 primary_dof_index) ==
711 dof_identities.end()) ||
712 (dof_identities
713 [dependent_dof_index] ==
714 primary_dof_index),
715 ExcInternalError());
716
717 dof_identities[dependent_dof_index] =
718 primary_dof_index;
719 }
720 }
721 }
722 }
723 }
724 }
725
726 return dof_identities;
727 }
728
729
730
735 template <int dim, int spacedim>
736 static std::map<types::global_dof_index, types::global_dof_index>
737 compute_quad_dof_identities(
738 const DoFHandler<dim, spacedim> &dof_handler)
739 {
740 (void)dof_handler;
741 Assert(
742 dof_handler.hp_capability_enabled == true,
743 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
744
745 // this function should only be called for dim<3 where there are
746 // no quad dof identities. for dim==3, the specialization below should
747 // take care of it
748 Assert(dim < 3, ExcInternalError());
749
750 return std::map<types::global_dof_index, types::global_dof_index>();
751 }
752
753
754 template <int spacedim>
755 static std::map<types::global_dof_index, types::global_dof_index>
756 compute_quad_dof_identities(const DoFHandler<3, spacedim> &dof_handler)
757 {
758 Assert(dof_handler.hp_capability_enabled == true,
759 (typename DoFHandler<3, spacedim>::ExcOnlyAvailableWithHP()));
760
761 const int dim = 3;
762
763 std::map<types::global_dof_index, types::global_dof_index>
764 dof_identities;
765
766 // An implementation of the algorithm described in the hp-
767 // paper, including the modification mentioned later in the
768 // "complications in 3-d" subsections
769 //
770 // as explained there, we do something only if there are
771 // exactly 2 finite elements associated with an object. if
772 // there is only one, then there is nothing to do anyway,
773 // and if there are 3 or more, then we can get into
774 // trouble. note that this only happens for lines in 3d and
775 // higher, and for quads only in 4d and higher, so this
776 // isn't a particularly frequent case
777 ::Table<3, std::unique_ptr<DoFIdentities>> quad_dof_identities(
778 dof_handler.fe_collection.size(),
779 dof_handler.fe_collection.size(),
780 2 /*triangle (0) or quadrilateral (1)*/);
781
782 std::vector<bool> quad_touched(
783 dof_handler.get_triangulation().n_raw_quads());
784 for (const auto &cell : dof_handler.active_cell_iterators())
785 for (const auto q : cell->face_indices())
786 if (!quad_touched[cell->quad(q)->index()] &&
787 (cell->quad(q)->n_active_fe_indices() == 2))
788 {
789 const auto quad = cell->quad(q);
790 quad_touched[quad->index()] = true;
791
792 const std::set<types::fe_index> fe_indices =
793 quad->get_active_fe_indices();
794
795 // find out which is the most dominating finite
796 // element of the ones that are used on this quad
797 // TODO: Change set to types::fe_index
798 const types::fe_index most_dominating_fe_index =
799 dof_handler.get_fe_collection().find_dominating_fe(
800 {fe_indices.begin(), fe_indices.end()},
801 /*codim=*/dim - 2);
802
803 // check if this cell is the dominating one and get the face
804 // indices
805 const bool this_cell_is_dominating =
806 cell->active_fe_index() == most_dominating_fe_index;
807
808 const unsigned int most_dominating_fe_index_face_no =
809 this_cell_is_dominating ? q : cell->neighbor_face_no(q);
810
811 const unsigned int other_fe_index_face_no =
812 this_cell_is_dominating ? cell->neighbor_face_no(q) : q;
813
814 // if we found the most dominating element, then use
815 // this to eliminate some of the degrees of freedom
816 // by identification. otherwise, the code that
817 // computes hanging node constraints will have to
818 // deal with it by computing appropriate constraints
819 // along this face/edge
820 if (most_dominating_fe_index != numbers::invalid_fe_index)
821 {
822 // loop over the indices of all the finite
823 // elements that are not dominating, and
824 // identify their dofs to the most dominating
825 // one
826 for (const auto &other_fe_index : fe_indices)
827 if (other_fe_index != most_dominating_fe_index)
828 {
829 const auto &identities =
830 *ensure_existence_and_return_dof_identities<2>(
831 dof_handler.get_fe_collection(),
832 most_dominating_fe_index,
833 other_fe_index,
834 quad_dof_identities
835 [most_dominating_fe_index][other_fe_index]
836 [cell->quad(q)->reference_cell() ==
837 ReferenceCells::Quadrilateral],
838 most_dominating_fe_index_face_no,
839 other_fe_index_face_no);
840
841 for (const auto &identity : identities)
842 {
843 const types::global_dof_index
844 primary_dof_index =
845 quad->dof_index(identity.first,
846 most_dominating_fe_index);
847 const types::global_dof_index
848 dependent_dof_index =
849 quad->dof_index(identity.second,
850 other_fe_index);
851
852 // we only store an identity if we are about to
853 // overwrite a valid degree of freedom. we will
854 // skip invalid degrees of freedom (that are
855 // associated with ghost cells) for now, and
856 // consider them later in phase 5.
857 if (dependent_dof_index !=
858 numbers::invalid_dof_index)
859 {
860 // if the DoF indices of both elements are
861 // already distributed, i.e., both of these
862 // 'fe_indices' are associated with a
863 // locally owned cell, then we should either
864 // not have a dof_identity yet, or it must
865 // come out here to be exactly as we had
866 // computed before
867 if (primary_dof_index !=
868 numbers::invalid_dof_index)
869 Assert((dof_identities.find(
870 primary_dof_index) ==
871 dof_identities.end()) ||
872 (dof_identities
873 [dependent_dof_index] ==
874 primary_dof_index),
875 ExcInternalError());
876
877 dof_identities[dependent_dof_index] =
878 primary_dof_index;
879 }
880 }
881 }
882 }
883 }
884
885 return dof_identities;
886 }
887
888
889
894 template <int dim, int spacedim>
895 static void
896 compute_dof_identities(std::vector<std::map<types::global_dof_index,
897 types::global_dof_index>>
898 &all_constrained_indices,
899 const DoFHandler<dim, spacedim> &dof_handler)
900 {
901 if (dof_handler.hp_capability_enabled == false)
902 return;
903
904 Assert(all_constrained_indices.size() == dim, ExcInternalError());
905
906 Threads::TaskGroup<> tasks;
907
908 unsigned int i = 0;
909 tasks += Threads::new_task([&, i]() {
910 all_constrained_indices[i] =
911 compute_vertex_dof_identities(dof_handler);
912 });
913
914 if (dim > 1)
915 {
916 ++i;
917 tasks += Threads::new_task([&, i]() {
918 all_constrained_indices[i] =
919 compute_line_dof_identities(dof_handler);
920 });
921 }
922
923 if (dim > 2)
924 {
925 ++i;
926 tasks += Threads::new_task([&, i]() {
927 all_constrained_indices[i] =
928 compute_quad_dof_identities(dof_handler);
929 });
930 }
931
932 tasks.join_all();
933 }
934
935
936
956 static types::global_dof_index
957 enumerate_dof_indices_for_renumbering(
958 std::vector<types::global_dof_index> &new_dof_indices,
959 const std::vector<
960 std::map<types::global_dof_index, types::global_dof_index>>
961 &all_constrained_indices,
962 const types::global_dof_index start_dof_index)
963 {
964 // first preset the new DoF indices that are identities
965 for (const auto &constrained_dof_indices : all_constrained_indices)
966 for (const auto &p : constrained_dof_indices)
967 if (new_dof_indices[p.first] != numbers::invalid_dof_index)
968 {
969 Assert(new_dof_indices[p.first] == enumeration_dof_index,
970 ExcInternalError());
971
972 new_dof_indices[p.first] = p.second;
973 }
974
975 // then enumerate the rest
976 types::global_dof_index next_free_dof = start_dof_index;
977 for (auto &new_dof_index : new_dof_indices)
978 if (new_dof_index == enumeration_dof_index)
979 new_dof_index = next_free_dof++;
980
981 // then loop over all those that are constrained and record the
982 // new dof number for those
983 for (const auto &constrained_dof_indices : all_constrained_indices)
984 for (const auto &p : constrained_dof_indices)
985 if (new_dof_indices[p.first] != numbers::invalid_dof_index)
986 {
987 Assert(new_dof_indices[p.first] != enumeration_dof_index,
988 ExcInternalError());
989
990 if (p.second != numbers::invalid_dof_index)
991 new_dof_indices[p.first] = new_dof_indices[p.second];
992 }
993
994 for (const types::global_dof_index new_dof_index : new_dof_indices)
995 {
996 (void)new_dof_index;
997 Assert(new_dof_index != enumeration_dof_index,
998 ExcInternalError());
999 Assert(new_dof_index < next_free_dof ||
1000 new_dof_index == numbers::invalid_dof_index,
1001 ExcInternalError());
1002 }
1003
1004 return next_free_dof;
1005 }
1006
1007
1008
1017 template <int dim, int spacedim>
1018 static types::global_dof_index
1019 unify_dof_indices(
1020 const DoFHandler<dim, spacedim> &dof_handler,
1021 const types::global_dof_index n_dofs_before_identification,
1022 const bool check_validity)
1023 {
1024 if (dof_handler.hp_capability_enabled == false)
1025 return n_dofs_before_identification;
1026
1027 std::vector<
1028 std::map<types::global_dof_index, types::global_dof_index>>
1029 all_constrained_indices(dim);
1030 compute_dof_identities(all_constrained_indices, dof_handler);
1031
1032 std::vector<::types::global_dof_index> renumbering(
1033 n_dofs_before_identification, enumeration_dof_index);
1034 const types::global_dof_index n_dofs =
1035 enumerate_dof_indices_for_renumbering(renumbering,
1036 all_constrained_indices,
1037 0);
1038
1039 renumber_dofs(renumbering, IndexSet(0), dof_handler, check_validity);
1040
1041 return n_dofs;
1042 }
1043
1044
1045
1050 template <int dim, int spacedim>
1051 static void
1052 merge_invalid_vertex_dofs_on_ghost_interfaces(
1053 DoFHandler<dim, spacedim> &dof_handler)
1054 {
1055 Assert(
1056 dof_handler.hp_capability_enabled == true,
1057 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
1058
1059 // Note: we may wish to have something here similar to what
1060 // we do for lines and quads, namely that we only identify
1061 // dofs for any FE towards the most dominating one. however,
1062 // it is not clear whether this is actually necessary for
1063 // vertices at all, I can't think of a finite element that
1064 // would make that necessary...
1065 ::Table<2, std::unique_ptr<DoFIdentities>>
1066 vertex_dof_identities(dof_handler.get_fe_collection().size(),
1067 dof_handler.get_fe_collection().size());
1068
1069 // mark all vertices on ghost cells to identify those cells that we
1070 // have already treated
1071 std::vector<bool> include_vertex(
1072 dof_handler.get_triangulation().n_vertices(), false);
1073 if (dynamic_cast<const ::parallel::
1074 DistributedTriangulationBase<dim, spacedim> *>(
1075 &dof_handler.get_triangulation()) != nullptr)
1076 for (const auto &cell : dof_handler.active_cell_iterators())
1077 if (cell->is_ghost())
1078 for (const unsigned int v : cell->vertex_indices())
1079 include_vertex[cell->vertex_index(v)] = true;
1080
1081 // loop over all vertices and see which one we need to work on
1082 for (unsigned int vertex_index = 0;
1083 vertex_index < dof_handler.get_triangulation().n_vertices();
1084 ++vertex_index)
1085 if ((dof_handler.get_triangulation()
1086 .get_used_vertices()[vertex_index] == true) &&
1087 (include_vertex[vertex_index] == true))
1088 {
1089 const unsigned int n_active_fe_indices =
1090 ::internal::DoFAccessorImplementation::Implementation::
1091 n_active_fe_indices(dof_handler,
1092 0,
1093 vertex_index,
1094 std::integral_constant<int, 0>());
1095
1096 if (n_active_fe_indices > 1)
1097 {
1098 const std::set<types::fe_index> fe_indices =
1099 ::internal::DoFAccessorImplementation::
1100 Implementation::get_active_fe_indices(
1101 dof_handler,
1102 0,
1103 vertex_index,
1104 std::integral_constant<int, 0>());
1105
1106 // find out which is the most dominating finite
1107 // element of the ones that are used on this vertex
1108 // TODO: Change set to types::fe_index
1109 types::fe_index most_dominating_fe_index =
1110 dof_handler.get_fe_collection().find_dominating_fe(
1111 {fe_indices.begin(), fe_indices.end()},
1112 /*codim=*/dim);
1113
1114 // if we haven't found a dominating finite element,
1115 // choose the very first one to be dominant similar
1116 // to compute_vertex_dof_identities()
1117 if (most_dominating_fe_index == numbers::invalid_fe_index)
1118 most_dominating_fe_index =
1119 ::internal::DoFAccessorImplementation::
1120 Implementation::nth_active_fe_index(
1121 dof_handler,
1122 0,
1123 vertex_index,
1124 0,
1125 std::integral_constant<int, 0>());
1126
1127 // loop over the indices of all the finite
1128 // elements that are not dominating, and
1129 // identify their dofs to the most dominating
1130 // one
1131 for (const auto &other_fe_index : fe_indices)
1132 if (other_fe_index != most_dominating_fe_index)
1133 {
1134 // make sure the entry in the equivalence
1135 // table exists
1136 const auto &identities =
1137 *ensure_existence_and_return_dof_identities<0>(
1138 dof_handler.get_fe_collection(),
1139 most_dominating_fe_index,
1140 other_fe_index,
1141 vertex_dof_identities[most_dominating_fe_index]
1142 [other_fe_index]);
1143
1144 // then loop through the identities we
1145 // have. first get the global numbers of the
1146 // dofs we want to identify and make sure they
1147 // are not yet constrained to anything else,
1148 // except for to each other. use the rule that
1149 // we will always constrain the dof with the
1150 // higher FE index to the one with the lower,
1151 // to avoid circular reasoning.
1152 for (const auto &identity : identities)
1153 {
1154 const types::global_dof_index primary_dof_index =
1155 ::internal::DoFAccessorImplementation::
1156 Implementation::get_dof_index(
1157 dof_handler,
1158 0,
1159 vertex_index,
1160 most_dominating_fe_index,
1161 identity.first,
1162 std::integral_constant<int, 0>());
1163 const types::global_dof_index
1164 dependent_dof_index =
1165 ::internal::DoFAccessorImplementation::
1166 Implementation::get_dof_index(
1167 dof_handler,
1168 0,
1169 vertex_index,
1170 other_fe_index,
1171 identity.second,
1172 std::integral_constant<int, 0>());
1173
1174 // check if we are on an interface between
1175 // a locally owned and a ghost cell on which
1176 // we need to work on.
1177 //
1178 // all degrees of freedom belonging to
1179 // dominating FE indices or to a processor
1180 // with a higher rank have been set at this
1181 // point (either in Phase 2, or after the
1182 // first ghost exchange in Phase 5). thus,
1183 // we only have to set the indices of
1184 // degrees of freedom that have been
1185 // previously flagged invalid.
1186 if ((dependent_dof_index ==
1187 numbers::invalid_dof_index) &&
1188 (primary_dof_index !=
1189 numbers::invalid_dof_index))
1190 ::internal::DoFAccessorImplementation::
1191 Implementation::set_dof_index(
1192 dof_handler,
1193 0,
1194 vertex_index,
1195 other_fe_index,
1196 identity.second,
1197 std::integral_constant<int, 0>(),
1198 primary_dof_index);
1199 }
1200 }
1201 }
1202 }
1203 }
1204
1205
1206
1211 template <int spacedim>
1212 static void
1213 merge_invalid_line_dofs_on_ghost_interfaces(
1214 DoFHandler<1, spacedim> &dof_handler)
1215 {
1216 (void)dof_handler;
1217 Assert(dof_handler.hp_capability_enabled == true,
1218 (typename DoFHandler<1, spacedim>::ExcOnlyAvailableWithHP()));
1219 }
1220
1221
1222 template <int dim, int spacedim>
1223 static void
1224 merge_invalid_line_dofs_on_ghost_interfaces(
1225 DoFHandler<dim, spacedim> &dof_handler)
1226 {
1227 Assert(
1228 dof_handler.hp_capability_enabled == true,
1229 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
1230
1231 // mark all lines on ghost cells
1232 std::vector<bool> line_marked(
1233 dof_handler.get_triangulation().n_raw_lines());
1234 for (const auto &cell : dof_handler.active_cell_iterators())
1235 if (cell->is_ghost())
1236 for (const auto l : cell->line_indices())
1237 line_marked[cell->line(l)->index()] = true;
1238
1239 // An implementation of the algorithm described in the hp-paper,
1240 // including the modification mentioned later in the "complications in
1241 // 3-d" subsections
1242 //
1243 // as explained there, we do something only if there are exactly 2
1244 // finite elements associated with an object. if there is only one,
1245 // then there is nothing to do anyway, and if there are 3 or more,
1246 // then we can get into trouble. note that this only happens for lines
1247 // in 3d and higher, and for quads only in 4d and higher, so this
1248 // isn't a particularly frequent case
1249 //
1250 // there is one case, however, that we would like to handle (see, for
1251 // example, the hp/crash_15 testcase): if we have
1252 // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
1253 // then we should be able to handle this because we can simply unify
1254 // *all* dofs, not only a some. so what we do is to first treat all
1255 // pairs of finite elements that have *identical* dofs, and then only
1256 // deal with those that are not identical of which we can handle at
1257 // most 2
1258 ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
1259 dof_handler.fe_collection.size(), dof_handler.fe_collection.size());
1260
1261 for (const auto &cell : dof_handler.active_cell_iterators())
1262 for (const auto l : cell->line_indices())
1263 if ((cell->is_locally_owned()) &&
1264 line_marked[cell->line(l)->index()])
1265 {
1266 const auto line = cell->line(l);
1267 line_marked[line->index()] = false;
1268
1269 unsigned int unique_sets_of_dofs =
1270 line->n_active_fe_indices();
1271
1272 // do a first loop over all sets of dofs and do identity
1273 // uniquification
1274 const unsigned int n_active_fe_indices =
1275 line->n_active_fe_indices();
1276 for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1277 for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
1278 {
1279 const types::fe_index fe_index_1 =
1280 line->nth_active_fe_index(f),
1281 fe_index_2 =
1282 line->nth_active_fe_index(g);
1283
1284 if ((dof_handler.get_fe(fe_index_1).n_dofs_per_line() ==
1285 dof_handler.get_fe(fe_index_2)
1286 .n_dofs_per_line()) &&
1287 (dof_handler.get_fe(fe_index_1).n_dofs_per_line() >
1288 0))
1289 {
1290 // the number of dofs per line is identical
1291 const unsigned int dofs_per_line =
1292 dof_handler.get_fe(fe_index_1).n_dofs_per_line();
1293
1294 const auto &identities =
1295 *ensure_existence_and_return_dof_identities<1>(
1296 dof_handler.get_fe_collection(),
1297 fe_index_1,
1298 fe_index_2,
1299 line_dof_identities[fe_index_1][fe_index_2]);
1300 // see if these sets of dofs are identical. the
1301 // first condition for this is that indeed there are
1302 // n identities
1303 if (identities.size() == dofs_per_line)
1304 {
1305 unsigned int i = 0;
1306 for (; i < dofs_per_line; ++i)
1307 if ((identities[i].first != i) &&
1308 (identities[i].second != i))
1309 // not an identity
1310 break;
1311
1312 if (i == dofs_per_line)
1313 {
1314 // The line dofs (i.e., the ones interior to
1315 // a line) of these two finite elements are
1316 // identical. Note that there could be
1317 // situations when one element still
1318 // dominates another, e.g.: FE_Q(2) x
1319 // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)
1320
1321 --unique_sets_of_dofs;
1322
1323 // determine which one of both finite
1324 // elements is the dominating one.
1325 const std::set<types::fe_index> fe_indices{
1326 fe_index_1, fe_index_2};
1327
1328 // TODO: Change set to types::fe_index
1329 types::fe_index dominating_fe_index =
1330 dof_handler.get_fe_collection()
1331 .find_dominating_fe({fe_indices.begin(),
1332 fe_indices.end()},
1333 /*codim*/ dim - 1);
1334 types::fe_index other_fe_index =
1335 numbers::invalid_fe_index;
1336
1337 if (dominating_fe_index !=
1338 numbers::invalid_fe_index)
1339 other_fe_index =
1340 (dominating_fe_index == fe_index_1) ?
1341 fe_index_2 :
1342 fe_index_1;
1343 else
1344 {
1345 // if we haven't found a dominating
1346 // finite element, choose the one with
1347 // the lower index to be dominating
1348 dominating_fe_index = fe_index_1;
1349 other_fe_index = fe_index_2;
1350 }
1351
1352 for (unsigned int j = 0; j < dofs_per_line;
1353 ++j)
1354 {
1355 const types::global_dof_index
1356 primary_dof_index = line->dof_index(
1357 j, dominating_fe_index);
1358 const types::global_dof_index
1359 dependent_dof_index =
1360 line->dof_index(j, other_fe_index);
1361
1362 // check if we are on an interface
1363 // between a locally owned and a ghost
1364 // cell on which we need to work on.
1365 //
1366 // all degrees of freedom belonging to
1367 // dominating fe_indices or to a
1368 // processor with a higher rank have
1369 // been set at this point (either in
1370 // Phase 2, or after the first ghost
1371 // exchange in Phase 5). thus, we only
1372 // have to set the indices of degrees
1373 // of freedom that have been previously
1374 // flagged invalid.
1375 if ((dependent_dof_index ==
1376 numbers::invalid_dof_index) &&
1377 (primary_dof_index !=
1378 numbers::invalid_dof_index))
1379 line->set_dof_index(j,
1380 primary_dof_index,
1381 other_fe_index);
1382 }
1383 }
1384 }
1385 }
1386 }
1387
1388 // if at this point, there is only one unique set of dofs
1389 // left, then we have taken care of everything above. if there
1390 // are two, then we need to deal with them here. if there are
1391 // more, then we punt, as described in the paper (and
1392 // mentioned above)
1393 // TODO: The check for 'dim==2' was inserted by intuition. It
1394 // fixes
1395 // the previous problems with @ref step_27 "step-27" in 3d. But an
1396 // explanation for this is still required, and what we do here
1397 // is not what we describe in the paper!.
1398 if ((unique_sets_of_dofs == 2) && (dim == 2))
1399 {
1400 const std::set<types::fe_index> fe_indices =
1401 line->get_active_fe_indices();
1402
1403 // find out which is the most dominating finite element of
1404 // the ones that are used on this line
1405 // TODO: Change set to types::fe_index
1406 const types::fe_index most_dominating_fe_index =
1407 dof_handler.get_fe_collection().find_dominating_fe(
1408 {fe_indices.begin(), fe_indices.end()},
1409 /*codim=*/dim - 1);
1410
1411 // if we found the most dominating element, then use this
1412 // to eliminate some of the degrees of freedom by
1413 // identification. otherwise, the code that computes
1414 // hanging node constraints will have to deal with it by
1415 // computing appropriate constraints along this face/edge
1416 if (most_dominating_fe_index != numbers::invalid_fe_index)
1417 {
1418 // loop over the indices of all the finite elements
1419 // that are not dominating, and identify their dofs to
1420 // the most dominating one
1421 for (const auto &other_fe_index : fe_indices)
1422 if (other_fe_index != most_dominating_fe_index)
1423 {
1424 const auto &identities =
1425 *ensure_existence_and_return_dof_identities<
1426 1>(dof_handler.get_fe_collection(),
1427 most_dominating_fe_index,
1428 other_fe_index,
1429 line_dof_identities
1430 [most_dominating_fe_index]
1431 [other_fe_index]);
1432
1433 for (const auto &identity : identities)
1434 {
1435 const types::global_dof_index
1436 primary_dof_index = line->dof_index(
1437 identity.first,
1438 most_dominating_fe_index);
1439 const types::global_dof_index
1440 dependent_dof_index =
1441 line->dof_index(identity.second,
1442 other_fe_index);
1443
1444 // check if we are on an interface between
1445 // a locally owned and a ghost cell on which
1446 // we need to work on.
1447 //
1448 // all degrees of freedom belonging to
1449 // dominating FE indices or to a processor
1450 // with a higher rank have been set at this
1451 // point (either in Phase 2, or after the
1452 // first ghost exchange in Phase 5). thus,
1453 // we only have to set the indices of
1454 // degrees of freedom that have been
1455 // previously flagged invalid.
1456 if ((dependent_dof_index ==
1457 numbers::invalid_dof_index) &&
1458 (primary_dof_index !=
1459 numbers::invalid_dof_index))
1460 line->set_dof_index(identity.second,
1461 primary_dof_index,
1462 other_fe_index);
1463 }
1464 }
1465 }
1466 }
1467 }
1468 }
1469
1470
1471
1476 template <int dim, int spacedim>
1477 static void
1478 merge_invalid_quad_dofs_on_ghost_interfaces(
1479 DoFHandler<dim, spacedim> &dof_handler)
1480 {
1481 (void)dof_handler;
1482 Assert(
1483 dof_handler.hp_capability_enabled == true,
1484 (typename DoFHandler<dim, spacedim>::ExcOnlyAvailableWithHP()));
1485
1486 // this function should only be called for dim<3 where there are
1487 // no quad dof identities. for dim>=3, the specialization below should
1488 // take care of it
1489 Assert(dim < 3, ExcInternalError());
1490 }
1491
1492
1493 template <int spacedim>
1494 static void
1495 merge_invalid_quad_dofs_on_ghost_interfaces(
1496 DoFHandler<3, spacedim> &dof_handler)
1497 {
1498 Assert(dof_handler.hp_capability_enabled == true,
1499 (typename DoFHandler<3, spacedim>::ExcOnlyAvailableWithHP()));
1500
1501 const int dim = 3;
1502
1503 // mark all quads on ghost cells
1504 std::vector<bool> quad_marked(
1505 dof_handler.get_triangulation().n_raw_quads());
1506 for (const auto &cell : dof_handler.active_cell_iterators())
1507 if (cell->is_ghost())
1508 for (const auto q : cell->face_indices())
1509 quad_marked[cell->quad(q)->index()] = true;
1510
1511 // An implementation of the algorithm described in the hp-
1512 // paper, including the modification mentioned later in the
1513 // "complications in 3-d" subsections
1514 //
1515 // as explained there, we do something only if there are
1516 // exactly 2 finite elements associated with an object. if
1517 // there is only one, then there is nothing to do anyway,
1518 // and if there are 3 or more, then we can get into
1519 // trouble. note that this only happens for lines in 3d and
1520 // higher, and for quads only in 4d and higher, so this
1521 // isn't a particularly frequent case
1522 ::Table<3, std::unique_ptr<DoFIdentities>> quad_dof_identities(
1523 dof_handler.fe_collection.size(),
1524 dof_handler.fe_collection.size(),
1525 2 /*triangle (0) or quadrilateral (1)*/);
1526
1527 for (const auto &cell : dof_handler.active_cell_iterators())
1528 for (const auto q : cell->face_indices())
1529 if ((cell->is_locally_owned()) &&
1530 quad_marked[cell->quad(q)->index()] &&
1531 (cell->quad(q)->n_active_fe_indices() == 2))
1532 {
1533 const auto quad = cell->quad(q);
1534 quad_marked[quad->index()] = false;
1535
1536 const std::set<types::fe_index> fe_indices =
1537 quad->get_active_fe_indices();
1538
1539 // find out which is the most dominating finite
1540 // element of the ones that are used on this quad
1541 // TODO: Change set to types::fe_index
1542 const types::fe_index most_dominating_fe_index =
1543 dof_handler.get_fe_collection().find_dominating_fe(
1544 {fe_indices.begin(), fe_indices.end()},
1545 /*codim=*/dim - 2);
1546
1547 // check if this cell is the dominating one and get the face
1548 // indices
1549 const bool this_cell_is_dominating =
1550 cell->active_fe_index() == most_dominating_fe_index;
1551
1552 const unsigned int most_dominating_fe_index_face_no =
1553 this_cell_is_dominating ? q : cell->neighbor_face_no(q);
1554
1555 const unsigned int other_fe_index_face_no =
1556 this_cell_is_dominating ? cell->neighbor_face_no(q) : q;
1557 // if we found the most dominating element, then use
1558 // this to eliminate some of the degrees of freedom
1559 // by identification. otherwise, the code that
1560 // computes hanging node constraints will have to
1561 // deal with it by computing appropriate constraints
1562 // along this face/edge
1563 if (most_dominating_fe_index != numbers::invalid_fe_index)
1564 {
1565 // loop over the indices of all the finite
1566 // elements that are not dominating, and
1567 // identify their dofs to the most dominating
1568 // one
1569 for (const auto &other_fe_index : fe_indices)
1570 if (other_fe_index != most_dominating_fe_index)
1571 {
1572 const auto &identities =
1573 *ensure_existence_and_return_dof_identities<2>(
1574 dof_handler.get_fe_collection(),
1575 most_dominating_fe_index,
1576 other_fe_index,
1577 quad_dof_identities
1578 [most_dominating_fe_index][other_fe_index]
1579 [cell->quad(q)->reference_cell() ==
1580 ReferenceCells::Quadrilateral],
1581 most_dominating_fe_index_face_no,
1582 other_fe_index_face_no);
1583
1584 for (const auto &identity : identities)
1585 {
1586 const types::global_dof_index
1587 primary_dof_index =
1588 quad->dof_index(identity.first,
1589 most_dominating_fe_index);
1590 const types::global_dof_index
1591 dependent_dof_index =
1592 quad->dof_index(identity.second,
1593 other_fe_index);
1594
1595 // check if we are on an interface between
1596 // a locally owned and a ghost cell on which
1597 // we need to work on.
1598 //
1599 // all degrees of freedom belonging to
1600 // dominating FE indices or to a processor with
1601 // a higher rank have been set at this point
1602 // (either in Phase 2, or after the first ghost
1603 // exchange in Phase 5). thus, we only have to
1604 // set the indices of degrees of freedom that
1605 // have been previously flagged invalid.
1606 if ((dependent_dof_index ==
1607 numbers::invalid_dof_index) &&
1608 (primary_dof_index !=
1609 numbers::invalid_dof_index))
1610 quad->set_dof_index(identity.second,
1611 primary_dof_index,
1612 other_fe_index);
1613 }
1614 }
1615 }
1616 }
1617 }
1618
1619
1620
1633 template <int dim, int spacedim>
1634 static void
1635 merge_invalid_dof_indices_on_ghost_interfaces(
1636 DoFHandler<dim, spacedim> &dof_handler)
1637 {
1638 if (dof_handler.hp_capability_enabled == false)
1639 return;
1640
1641 {
1642 Threads::TaskGroup<> tasks;
1643
1644 tasks += Threads::new_task([&]() {
1645 merge_invalid_vertex_dofs_on_ghost_interfaces(dof_handler);
1646 });
1647
1648 if (dim > 1)
1649 {
1650 tasks += Threads::new_task([&]() {
1651 merge_invalid_line_dofs_on_ghost_interfaces(dof_handler);
1652 });
1653 }
1654
1655 if (dim > 2)
1656 {
1657 tasks += Threads::new_task([&]() {
1658 merge_invalid_quad_dofs_on_ghost_interfaces(dof_handler);
1659 });
1660 }
1661
1662 tasks.join_all();
1663 }
1664 }
1665
1666
1667
1674 template <int dim, int spacedim>
1675 static types::global_dof_index
1676 distribute_dofs(const types::subdomain_id subdomain_id,
1677 DoFHandler<dim, spacedim> &dof_handler)
1678 {
1679 Assert(dof_handler.get_triangulation().n_levels() > 0,
1680 ExcMessage("Empty triangulation"));
1681
1682 // distribute dofs on all cells excluding artificial ones
1683 types::global_dof_index next_free_dof = 0;
1684
1685 for (auto cell : dof_handler.active_cell_iterators())
1686 if (!cell->is_artificial() &&
1687 ((subdomain_id == numbers::invalid_subdomain_id) ||
1688 (cell->subdomain_id() == subdomain_id)))
1689 {
1690 // feed the process_dof_indices function with an empty type
1691 // `std::tuple<>`, as we do not want to retrieve any DoF
1692 // indices here and rather modify the stored ones
1693 DoFAccessorImplementation::Implementation::process_dof_indices(
1694 *cell,
1695 std::make_tuple(),
1696 cell->active_fe_index(),
1697 DoFAccessorImplementation::Implementation::
1698 DoFIndexProcessor<dim, spacedim>(),
1699 [&next_free_dof](auto &stored_index, auto) {
1700 if (stored_index == numbers::invalid_dof_index)
1701 {
1702 stored_index = next_free_dof;
1703 Assert(
1704 next_free_dof !=
1705 std::numeric_limits<types::global_dof_index>::max(),
1706 ExcMessage(
1707 "You have reached the maximal number of degrees of "
1708 "freedom that can be stored in the chosen data "
1709 "type. In practice, this can only happen if you "
1710 "are using 32-bit data types. You will have to "
1711 "re-compile deal.II with the "
1712 "`DEAL_II_WITH_64BIT_INDICES' flag set to `ON'."));
1713 ++next_free_dof;
1714 }
1715 },
1716 false);
1717 }
1718
1719 return next_free_dof;
1720 }
1721
1722
1723
1737 template <int dim, int spacedim>
1738 static void
1739 invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering(
1740 std::vector<types::global_dof_index> &renumbering,
1741 const types::subdomain_id subdomain_id,
1742 const DoFHandler<dim, spacedim> &dof_handler)
1743 {
1744 std::vector<types::global_dof_index> local_dof_indices;
1745
1746 for (const auto &cell : dof_handler.active_cell_iterators())
1747 if (cell->is_ghost() && (cell->subdomain_id() < subdomain_id))
1748 {
1749 // we found a neighboring ghost cell whose subdomain
1750 // is "stronger" than our own subdomain
1751
1752 // delete all dofs that live there and that we have
1753 // previously assigned a number to (i.e. the ones on
1754 // the interface); make sure to not use the cache
1755 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1756 internal::DoFAccessorImplementation::Implementation::
1757 get_dof_indices(*cell,
1758 local_dof_indices,
1759 cell->active_fe_index());
1760 for (const auto &local_dof_index : local_dof_indices)
1761 if (local_dof_index != numbers::invalid_dof_index)
1762 renumbering[local_dof_index] = numbers::invalid_dof_index;
1763 }
1764 }
1765
1766
1767
1768 /* -------------- distribute_mg_dofs functionality ------------- */
1769
1770
1771
1772 template <int dim, int spacedim>
1773 static types::global_dof_index
1774 distribute_dofs_on_level(const types::subdomain_id level_subdomain_id,
1775 DoFHandler<dim, spacedim> &dof_handler,
1776 const unsigned int level)
1777 {
1778 Assert(dof_handler.hp_capability_enabled == false,
1779 ExcInternalError());
1780
1781 const ::Triangulation<dim, spacedim> &tria =
1782 dof_handler.get_triangulation();
1783 Assert(tria.n_levels() > 0, ExcMessage("Empty triangulation"));
1784 if (level >= tria.n_levels())
1785 return 0; // this is allowed for multigrid
1786
1787 types::global_dof_index next_free_dof = 0;
1788
1789 for (auto cell : dof_handler.cell_iterators_on_level(level))
1790 if ((level_subdomain_id == numbers::invalid_subdomain_id) ||
1791 (cell->level_subdomain_id() == level_subdomain_id))
1792 {
1793 DoFAccessorImplementation::Implementation::process_dof_indices(
1794 *cell,
1795 std::make_tuple(),
1796 0,
1797 DoFAccessorImplementation::Implementation::
1798 MGDoFIndexProcessor<dim, spacedim>(level),
1799 [&next_free_dof](auto &stored_index, auto) {
1800 if (stored_index == numbers::invalid_dof_index)
1801 {
1802 stored_index = next_free_dof;
1803 Assert(
1804 next_free_dof !=
1805 std::numeric_limits<types::global_dof_index>::max(),
1806 ExcMessage(
1807 "You have reached the maximal number of degrees of "
1808 "freedom that can be stored in the chosen data "
1809 "type. In practice, this can only happen if you "
1810 "are using 32-bit data types. You will have to "
1811 "re-compile deal.II with the "
1812 "`DEAL_II_WITH_64BIT_INDICES' flag set to `ON'."));
1813 ++next_free_dof;
1814 }
1815 },
1816 true);
1817 }
1818
1819 return next_free_dof;
1820 }
1821
1822
1823
1824 /* --------------------- renumber_dofs functionality ---------------- */
1825
1826
1834 template <int dim, int spacedim>
1835 static void
1836 renumber_face_dofs(
1837 const std::vector<types::global_dof_index> &new_numbers,
1838 const IndexSet &indices_we_care_about,
1839 DoFHandler<dim, spacedim> &dof_handler)
1840 {
1841 for (unsigned int d = 1; d < dim; ++d)
1842 for (auto &i : dof_handler.object_dof_indices[0][d])
1843 if (i != numbers::invalid_dof_index)
1844 i = ((indices_we_care_about.size() == 0) ?
1845 new_numbers[i] :
1846 new_numbers[indices_we_care_about.index_within_set(i)]);
1847 }
1848
1849
1850
1851 template <int dim, int spacedim>
1852 static void
1853 renumber_vertex_dofs(
1854 const std::vector<types::global_dof_index> &new_numbers,
1855 const IndexSet &indices_we_care_about,
1856 DoFHandler<dim, spacedim> &dof_handler,
1857 const bool check_validity)
1858 {
1859 if (dof_handler.hp_capability_enabled == false)
1860 {
1861 // we can not use cell iterators in this function since then
1862 // we would renumber the dofs on the interface of two cells
1863 // more than once. Anyway, this way it's not only more
1864 // correct but also faster; note, however, that dof numbers
1865 // may be invalid_dof_index, namely when the appropriate
1866 // vertex/line/etc is unused
1867 for (std::vector<types::global_dof_index>::iterator i =
1868 dof_handler.object_dof_indices[0][0].begin();
1869 i != dof_handler.object_dof_indices[0][0].end();
1870 ++i)
1871 if (*i != numbers::invalid_dof_index)
1872 *i =
1873 (indices_we_care_about.size() == 0) ?
1874 (new_numbers[*i]) :
1875 (new_numbers[indices_we_care_about.index_within_set(*i)]);
1876 else if (check_validity)
1877 // if index is invalid_dof_index: check if this one
1878 // really is unused
1879 Assert(dof_handler.get_triangulation().vertex_used(
1880 (i - dof_handler.object_dof_indices[0][0].begin()) /
1881 dof_handler.get_fe().n_dofs_per_vertex()) == false,
1882 ExcInternalError());
1883 return;
1884 }
1885
1886
1887 for (unsigned int vertex_index = 0;
1888 vertex_index < dof_handler.get_triangulation().n_vertices();
1889 ++vertex_index)
1890 {
1891 const unsigned int n_active_fe_indices =
1892 ::internal::DoFAccessorImplementation::Implementation::
1893 n_active_fe_indices(dof_handler,
1894 0,
1895 vertex_index,
1896 std::integral_constant<int, 0>());
1897
1898 // if this vertex is unused, then we really ought not to have
1899 // allocated any space for it, i.e., n_active_fe_indices should be
1900 // zero, and there is no space to actually store dof indices for
1901 // this vertex
1902 if (dof_handler.get_triangulation().vertex_used(vertex_index) ==
1903 false)
1904 Assert(n_active_fe_indices == 0, ExcInternalError());
1905
1906 // otherwise the vertex is used; it may still not hold any dof
1907 // indices if it is located on an artificial cell and not adjacent
1908 // to a ghost cell, but in that case there is simply nothing for
1909 // us to do
1910 for (unsigned int f = 0; f < n_active_fe_indices; ++f)
1911 {
1912 const types::fe_index fe_index =
1913 ::internal::DoFAccessorImplementation::
1914 Implementation::nth_active_fe_index(
1915 dof_handler,
1916 0,
1917 vertex_index,
1918 f,
1919 std::integral_constant<int, 0>());
1920
1921 for (unsigned int d = 0;
1922 d < dof_handler.get_fe(fe_index).n_dofs_per_vertex();
1923 ++d)
1924 {
1925 const types::global_dof_index old_dof_index =
1926 ::internal::DoFAccessorImplementation::
1927 Implementation::get_dof_index(
1928 dof_handler,
1929 0,
1930 vertex_index,
1931 fe_index,
1932 d,
1933 std::integral_constant<int, 0>());
1934
1935 // if check_validity was set, then we are to verify that
1936 // the previous indices were all valid. this really should
1937 // be the case: we allocated space for these vertex dofs,
1938 // i.e., at least one adjacent cell has a valid
1939 // active FE index, so there are DoFs that really live
1940 // on this vertex. if check_validity is set, then we
1941 // must make sure that they have been set to something
1942 // useful
1943 if (check_validity)
1944 Assert(old_dof_index != numbers::invalid_dof_index,
1945 ExcInternalError());
1946
1947 if (old_dof_index != numbers::invalid_dof_index)
1948 {
1949 // In the following blocks, we first check whether
1950 // we were given an IndexSet of DoFs to touch. If not
1951 // (the first 'if' case here), then we are in the
1952 // sequential case and are allowed to touch all DoFs.
1953 //
1954 // If yes (the 'else' case), then we need to
1955 // distinguish whether the DoF whose number we want to
1956 // touch is in fact locally owned (i.e., is in the
1957 // index set) and then we can actually assign it a new
1958 // number; otherwise, we have encountered a
1959 // non-locally owned DoF for which we don't know the
1960 // new number yet and so set it to an invalid index.
1961 // This will later be fixed up after the first ghost
1962 // exchange phase when we unify hp-DoFs on neighboring
1963 // cells.
1964 if (indices_we_care_about.size() == 0)
1965 ::internal::DoFAccessorImplementation::
1966 Implementation::set_dof_index(
1967 dof_handler,
1968 0,
1969 vertex_index,
1970 fe_index,
1971 d,
1972 std::integral_constant<int, 0>(),
1973 new_numbers[old_dof_index]);
1974 else
1975 {
1976 if (indices_we_care_about.is_element(
1977 old_dof_index))
1978 ::internal::DoFAccessorImplementation::
1979 Implementation::set_dof_index(
1980 dof_handler,
1981 0,
1982 vertex_index,
1983 fe_index,
1984 d,
1985 std::integral_constant<int, 0>(),
1986 new_numbers[indices_we_care_about
1987 .index_within_set(
1988 old_dof_index)]);
1989 else
1990 ::internal::DoFAccessorImplementation::
1991 Implementation::set_dof_index(
1992 dof_handler,
1993 0,
1994 vertex_index,
1995 fe_index,
1996 d,
1997 std::integral_constant<int, 0>(),
1998 numbers::invalid_dof_index);
1999 }
2000 }
2001 }
2002 }
2003 }
2004 }
2005
2006
2007
2008 template <int dim, int spacedim>
2009 static void
2010 renumber_cell_dofs(
2011 const std::vector<types::global_dof_index> &new_numbers,
2012 const IndexSet &indices_we_care_about,
2013 DoFHandler<dim, spacedim> &dof_handler)
2014 {
2015 if (dof_handler.hp_capability_enabled == false)
2016 {
2017 for (unsigned int level = 0;
2018 level < dof_handler.object_dof_indices.size();
2019 ++level)
2020 for (auto &i : dof_handler.object_dof_indices[level][dim])
2021 if (i != numbers::invalid_dof_index)
2022 i = ((indices_we_care_about.size() == 0) ?
2023 new_numbers[i] :
2024 new_numbers[indices_we_care_about.index_within_set(
2025 i)]);
2026 return;
2027 }
2028
2029 for (const auto &cell : dof_handler.active_cell_iterators())
2030 if (!cell->is_artificial())
2031 {
2032 const types::fe_index fe_index = cell->active_fe_index();
2033
2034 for (unsigned int d = 0;
2035 d < dof_handler.get_fe(fe_index)
2036 .template n_dofs_per_object<dim>();
2037 ++d)
2038 {
2039 const types::global_dof_index old_dof_index =
2040 cell->dof_index(d, fe_index);
2041 if (old_dof_index != numbers::invalid_dof_index)
2042 {
2043 // In the following blocks, we first check whether
2044 // we were given an IndexSet of DoFs to touch. If not
2045 // (the first 'if' case here), then we are in the
2046 // sequential case and are allowed to touch all DoFs.
2047 //
2048 // If yes (the 'else' case), then we need to distinguish
2049 // whether the DoF whose number we want to touch is in
2050 // fact locally owned (i.e., is in the index set) and
2051 // then we can actually assign it a new number;
2052 // otherwise, we have encountered a non-locally owned
2053 // DoF for which we don't know the new number yet and so
2054 // set it to an invalid index. This will later be fixed
2055 // up after the first ghost exchange phase when we unify
2056 // hp-DoFs on neighboring cells.
2057 if (indices_we_care_about.size() == 0)
2058 cell->set_dof_index(d,
2059 new_numbers[old_dof_index],
2060 fe_index);
2061 else
2062 {
2063 if (indices_we_care_about.is_element(old_dof_index))
2064 cell->set_dof_index(
2065 d,
2066 new_numbers[indices_we_care_about
2067 .index_within_set(old_dof_index)],
2068 fe_index);
2069 else
2070 cell->set_dof_index(d,
2071 numbers::invalid_dof_index,
2072 fe_index);
2073 }
2074 }
2075 }
2076 }
2077 }
2078
2079
2080
2081 template <int spacedim>
2082 static void
2083 renumber_face_dofs(
2084 const std::vector<types::global_dof_index> & /*new_numbers*/,
2085 const IndexSet & /*indices_we_care_about*/,
2086 DoFHandler<1, spacedim> & /*dof_handler*/)
2087 {
2088 // nothing to do in 1d since there are no separate faces -- we've
2089 // already taken care of this when dealing with the vertices
2090 }
2091
2092
2093
2094 template <int spacedim>
2095 static void
2096 renumber_face_dofs(
2097 const std::vector<types::global_dof_index> &new_numbers,
2098 const IndexSet &indices_we_care_about,
2099 DoFHandler<2, spacedim> &dof_handler)
2100 {
2101 const unsigned int dim = 2;
2102
2103 if (dof_handler.hp_capability_enabled == false)
2104 {
2105 for (unsigned int d = 1; d < dim; ++d)
2106 for (auto &i : dof_handler.object_dof_indices[0][d])
2107 if (i != numbers::invalid_dof_index)
2108 i = ((indices_we_care_about.size() == 0) ?
2109 new_numbers[i] :
2110 new_numbers[indices_we_care_about.index_within_set(
2111 i)]);
2112 return;
2113 }
2114
2115 // deal with DoFs on lines
2116 {
2117 std::vector<bool> line_touched(
2118 dof_handler.get_triangulation().n_raw_lines());
2119 for (const auto &cell : dof_handler.active_cell_iterators())
2120 if (!cell->is_artificial())
2121 for (const auto l : cell->line_indices())
2122 if (!line_touched[cell->line(l)->index()])
2123 {
2124 const auto line = cell->line(l);
2125 line_touched[line->index()] = true;
2126
2127 const unsigned int n_active_fe_indices =
2128 line->n_active_fe_indices();
2129
2130 for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2131 {
2132 const types::fe_index fe_index =
2133 line->nth_active_fe_index(f);
2134
2135 for (unsigned int d = 0;
2136 d <
2137 dof_handler.get_fe(fe_index).n_dofs_per_line();
2138 ++d)
2139 {
2140 const types::global_dof_index old_dof_index =
2141 line->dof_index(d, fe_index);
2142 if (old_dof_index != numbers::invalid_dof_index)
2143 {
2144 // In the following blocks, we first check
2145 // whether we were given an IndexSet of DoFs
2146 // to touch. If not (the first 'if' case
2147 // here), then we are in the sequential case
2148 // and are allowed to touch all DoFs.
2149 //
2150 // If yes (the 'else' case), then we need to
2151 // distinguish whether the DoF whose number we
2152 // want to touch is in fact locally owned
2153 // (i.e., is in the index set) and then we can
2154 // actually assign it a new number; otherwise,
2155 // we have encountered a non-locally owned DoF
2156 // for which we don't know the new number yet
2157 // and so set it to an invalid index. This
2158 // will later be fixed up after the first
2159 // ghost exchange phase when we unify hp-DoFs
2160 // on neighboring cells.
2161 if (indices_we_care_about.size() == 0)
2162 line->set_dof_index(
2163 d, new_numbers[old_dof_index], fe_index);
2164 else
2165 {
2166 if (indices_we_care_about.is_element(
2167 old_dof_index))
2168 line->set_dof_index(
2169 d,
2170 new_numbers[indices_we_care_about
2171 .index_within_set(
2172 old_dof_index)],
2173 fe_index);
2174 else
2175 line->set_dof_index(
2176 d,
2177 numbers::invalid_dof_index,
2178 fe_index);
2179 }
2180 }
2181 }
2182 }
2183 }
2184 }
2185 }
2186
2187
2188
2189 template <int spacedim>
2190 static void
2191 renumber_face_dofs(
2192 const std::vector<types::global_dof_index> &new_numbers,
2193 const IndexSet &indices_we_care_about,
2194 DoFHandler<3, spacedim> &dof_handler)
2195 {
2196 const unsigned int dim = 3;
2197
2198 if (dof_handler.hp_capability_enabled == false)
2199 {
2200 for (unsigned int d = 1; d < dim; ++d)
2201 for (auto &i : dof_handler.object_dof_indices[0][d])
2202 if (i != numbers::invalid_dof_index)
2203 i = ((indices_we_care_about.size() == 0) ?
2204 new_numbers[i] :
2205 new_numbers[indices_we_care_about.index_within_set(
2206 i)]);
2207 return;
2208 }
2209
2210 // deal with DoFs on lines
2211 {
2212 std::vector<bool> line_touched(
2213 dof_handler.get_triangulation().n_raw_lines());
2214 for (const auto &cell : dof_handler.active_cell_iterators())
2215 if (!cell->is_artificial())
2216 for (const auto l : cell->line_indices())
2217 if (!line_touched[cell->line(l)->index()])
2218 {
2219 const auto line = cell->line(l);
2220 line_touched[line->index()] = true;
2221
2222 const unsigned int n_active_fe_indices =
2223 line->n_active_fe_indices();
2224
2225 for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2226 {
2227 const types::fe_index fe_index =
2228 line->nth_active_fe_index(f);
2229
2230 for (unsigned int d = 0;
2231 d <
2232 dof_handler.get_fe(fe_index).n_dofs_per_line();
2233 ++d)
2234 {
2235 const types::global_dof_index old_dof_index =
2236 line->dof_index(d, fe_index);
2237 if (old_dof_index != numbers::invalid_dof_index)
2238 {
2239 // In the following blocks, we first check
2240 // whether we were given an IndexSet of DoFs
2241 // to touch. If not (the first 'if' case
2242 // here), then we are in the sequential case
2243 // and are allowed to touch all DoFs.
2244 //
2245 // If yes (the 'else' case), then we need to
2246 // distinguish whether the DoF whose number we
2247 // want to touch is in fact locally owned
2248 // (i.e., is in the index set) and then we can
2249 // actually assign it a new number; otherwise,
2250 // we have encountered a non-locally owned DoF
2251 // for which we don't know the new number yet
2252 // and so set it to an invalid index. This
2253 // will later be fixed up after the first
2254 // ghost exchange phase when we unify hp-DoFs
2255 // on neighboring cells.
2256 if (indices_we_care_about.size() == 0)
2257 line->set_dof_index(
2258 d, new_numbers[old_dof_index], fe_index);
2259 else if (indices_we_care_about.is_element(
2260 old_dof_index))
2261 line->set_dof_index(
2262 d,
2263 new_numbers[indices_we_care_about
2264 .index_within_set(
2265 old_dof_index)],
2266 fe_index);
2267 else
2268 line->set_dof_index(
2269 d, numbers::invalid_dof_index, fe_index);
2270 }
2271 }
2272 }
2273 }
2274 }
2275
2276 // then deal with dofs on quads
2277 {
2278 std::vector<bool> quad_touched(
2279 dof_handler.get_triangulation().n_raw_quads());
2280 for (const auto &cell : dof_handler.active_cell_iterators())
2281 if (!cell->is_artificial())
2282 for (const auto q : cell->face_indices())
2283 if (!quad_touched[cell->quad(q)->index()])
2284 {
2285 const auto quad = cell->quad(q);
2286 quad_touched[quad->index()] = true;
2287
2288 const unsigned int n_active_fe_indices =
2289 quad->n_active_fe_indices();
2290
2291 for (unsigned int f = 0; f < n_active_fe_indices; ++f)
2292 {
2293 const types::fe_index fe_index =
2294 quad->nth_active_fe_index(f);
2295
2296 // figure out on which side of the face we are on
2297 const unsigned int face_no =
2298 cell->active_fe_index() == fe_index ?
2299 q :
2300 cell->neighbor_face_no(q);
2301
2302 for (unsigned int d = 0;
2303 d < dof_handler.get_fe(fe_index).n_dofs_per_quad(
2304 face_no);
2305 ++d)
2306 {
2307 const types::global_dof_index old_dof_index =
2308 quad->dof_index(d, fe_index);
2309 if (old_dof_index != numbers::invalid_dof_index)
2310 {
2311 // In the following blocks, we first check
2312 // whether we were given an IndexSet of DoFs
2313 // to touch. If not (the first 'if' case
2314 // here), then we are in the sequential case
2315 // and are allowed to touch all DoFs.
2316 //
2317 // If yes (the 'else' case), then we need to
2318 // distinguish whether the DoF whose number we
2319 // want to touch is in fact locally owned
2320 // (i.e., is in the index set) and then we can
2321 // actually assign it a new number; otherwise,
2322 // we have encountered a non-locally owned DoF
2323 // for which we don't know the new number yet
2324 // and so set it to an invalid index. This
2325 // will later be fixed up after the first
2326 // ghost exchange phase when we unify hp-DoFs
2327 // on neighboring cells.
2328 if (indices_we_care_about.size() == 0)
2329 quad->set_dof_index(
2330 d, new_numbers[old_dof_index], fe_index);
2331 else
2332 {
2333 if (indices_we_care_about.is_element(
2334 old_dof_index))
2335 quad->set_dof_index(
2336 d,
2337 new_numbers[indices_we_care_about
2338 .index_within_set(
2339 old_dof_index)],
2340 fe_index);
2341 else
2342 quad->set_dof_index(
2343 d,
2344 numbers::invalid_dof_index,
2345 fe_index);
2346 }
2347 }
2348 }
2349 }
2350 }
2351 }
2352 }
2353
2354
2355
2367 template <int dim, int space_dim>
2368 static void
2369 renumber_dofs(const std::vector<types::global_dof_index> &new_numbers,
2370 const IndexSet &indices_we_care_about,
2371 const DoFHandler<dim, space_dim> &dof_handler,
2372 const bool check_validity)
2373 {
2374 if (dim == 1)
2375 Assert(indices_we_care_about == IndexSet(0), ExcNotImplemented());
2376
2377 // renumber DoF indices on vertices, cells, and faces. this
2378 // can be done in parallel because the respective functions
2379 // work on separate data structures
2380 Threads::TaskGroup<> tasks;
2381 tasks += Threads::new_task([&]() {
2382 renumber_vertex_dofs(new_numbers,
2383 indices_we_care_about,
2384 const_cast<DoFHandler<dim, space_dim> &>(
2385 dof_handler),
2386 check_validity);
2387 });
2388 tasks += Threads::new_task([&]() {
2389 renumber_face_dofs(new_numbers,
2390 indices_we_care_about,
2391 const_cast<DoFHandler<dim, space_dim> &>(
2392 dof_handler));
2393 });
2394 tasks += Threads::new_task([&]() {
2395 renumber_cell_dofs(new_numbers,
2396 indices_we_care_about,
2397 const_cast<DoFHandler<dim, space_dim> &>(
2398 dof_handler));
2399 });
2400 tasks.join_all();
2401 }
2402
2403
2404
2405 /* --------------------- renumber_mg_dofs functionality ----------------
2406 */
2407
2415 template <int dim, int spacedim>
2416 static void
2417 renumber_vertex_mg_dofs(
2418 const std::vector<::types::global_dof_index> &new_numbers,
2419 const IndexSet &indices_we_care_about,
2420 DoFHandler<dim, spacedim> &dof_handler,
2421 const unsigned int level)
2422 {
2423 Assert(level < dof_handler.get_triangulation().n_levels(),
2424 ExcInternalError());
2425
2426 for (auto i = dof_handler.mg_vertex_dofs.begin();
2427 i != dof_handler.mg_vertex_dofs.end();
2428 ++i)
2429 // if the present vertex lives on the current level
2430 if ((i->get_coarsest_level() <= level) &&
2431 (i->get_finest_level() >= level))
2432 for (unsigned int d = 0;
2433 d < dof_handler.get_fe().n_dofs_per_vertex();
2434 ++d)
2435 {
2436 const ::types::global_dof_index idx =
2437 i->access_index(level,
2438 d,
2439 dof_handler.get_fe().n_dofs_per_vertex());
2440
2441 if (idx != numbers::invalid_dof_index)
2442 {
2443 Assert(indices_we_care_about.size() > 0 ?
2444 indices_we_care_about.is_element(idx) :
2445 (idx < new_numbers.size()),
2446 ExcInternalError());
2447 i->access_index(
2448 level, d, dof_handler.get_fe().n_dofs_per_vertex()) =
2449 (indices_we_care_about.size() == 0) ?
2450 new_numbers[idx] :
2451 new_numbers[indices_we_care_about.index_within_set(
2452 idx)];
2453 }
2454 }
2455 }
2456
2457
2458
2466 template <int dim, int spacedim>
2467 static void
2468 renumber_cell_mg_dofs(
2469 const std::vector<::types::global_dof_index> &new_numbers,
2470 const IndexSet &indices_we_care_about,
2471 DoFHandler<dim, spacedim> &dof_handler,
2472 const unsigned int level)
2473 {
2474 for (std::vector<types::global_dof_index>::iterator i =
2475 dof_handler.mg_levels[level]->dof_object.dofs.begin();
2476 i != dof_handler.mg_levels[level]->dof_object.dofs.end();
2477 ++i)
2478 {
2479 if (*i != numbers::invalid_dof_index)
2480 {
2481 Assert((indices_we_care_about.size() > 0 ?
2482 indices_we_care_about.is_element(*i) :
2483 (*i < new_numbers.size())),
2484 ExcInternalError());
2485 *i =
2486 (indices_we_care_about.size() == 0) ?
2487 (new_numbers[*i]) :
2488 (new_numbers[indices_we_care_about.index_within_set(*i)]);
2489 }
2490 }
2491 }
2492
2493
2494
2502 template <int spacedim>
2503 static void
2504 renumber_face_mg_dofs(
2505 const std::vector<types::global_dof_index> & /*new_numbers*/,
2506 const IndexSet & /*indices_we_care_about*/,
2507 DoFHandler<1, spacedim> & /*dof_handler*/,
2508 const unsigned int /*level*/,
2509 const bool /*check_validity*/)
2510 {
2511 // nothing to do in 1d because there are no separate faces
2512 }
2513
2514
2515
2516 template <int dim, int spacedim>
2517 static void
2518 renumber_face_mg_dofs(
2519 const std::vector<::types::global_dof_index> &new_numbers,
2520 const IndexSet &indices_we_care_about,
2521 DoFHandler<dim, spacedim> &dof_handler,
2522 const unsigned int level,
2523 const bool check_validity)
2524 {
2525 const unsigned int dofs_per_line =
2526 dof_handler.get_fe().n_dofs_per_line();
2527 if (dofs_per_line > 0 ||
2528 (dim > 2 && dof_handler.get_fe().max_dofs_per_quad() > 0))
2529 {
2530 // visit all lines/quads adjacent to cells of the current level
2531 // exactly once, as those lines/quads logically belong to the same
2532 // level as the cell, at least for isotropic refinement
2533 std::vector<bool> line_touched(
2534 dof_handler.get_triangulation().n_raw_lines());
2535 std::vector<bool> quad_touched(
2536 dim > 2 ? dof_handler.get_triangulation().n_raw_quads() : 0);
2537 for (const auto &cell :
2538 dof_handler.cell_iterators_on_level(level))
2539 if (cell->level_subdomain_id() !=
2540 numbers::artificial_subdomain_id)
2541 {
2542 // lines
2543 if (dofs_per_line > 0)
2544 {
2545 const auto line_indices =
2546 internal::TriaAccessorImplementation::Implementation::
2547 get_line_indices_of_cell(*cell);
2548 for (const auto line : cell->line_indices())
2549 {
2550 if (!line_touched[line_indices[line]])
2551 {
2552 line_touched[line_indices[line]] = true;
2553 ::types::global_dof_index *indices =
2554 &internal::DoFAccessorImplementation::
2555 Implementation::get_mg_dof_index(
2556 dof_handler,
2557 dof_handler.mg_levels[level],
2558 dof_handler.mg_faces,
2559 line_indices[line],
2560 0,
2561 0,
2562 std::integral_constant<int, 1>());
2563 for (unsigned int d = 0; d < dofs_per_line; ++d)
2564 {
2565 if (check_validity)
2566 Assert(indices[d] !=
2567 numbers::invalid_dof_index,
2568 ExcInternalError());
2569
2570 if (indices[d] !=
2571 numbers::invalid_dof_index)
2572 indices[d] =
2573 (indices_we_care_about.size() == 0) ?
2574 new_numbers[indices[d]] :
2575 new_numbers[indices_we_care_about
2576 .index_within_set(
2577 indices[d])];
2578 }
2579 }
2580 }
2581 }
2582
2583 // quads
2584 if (dim > 2)
2585 for (const auto quad : cell->face_indices())
2586 if (!quad_touched[cell->quad(quad)->index()])
2587 {
2588 quad_touched[cell->quad(quad)->index()] = true;
2589 const unsigned int dofs_per_quad =
2590 dof_handler.get_fe().n_dofs_per_quad(quad);
2591 if (dofs_per_quad > 0)
2592 {
2593 ::types::global_dof_index *indices =
2594 &internal::DoFAccessorImplementation::
2595 Implementation::get_mg_dof_index(
2596 dof_handler,
2597 dof_handler.mg_levels[level],
2598 dof_handler.mg_faces,
2599 cell->quad(quad)->index(),
2600 0,
2601 0,
2602 std::integral_constant<int, 2>());
2603 for (unsigned int d = 0; d < dofs_per_quad; ++d)
2604 {
2605 if (check_validity)
2606 Assert(indices[d] !=
2607 numbers::invalid_dof_index,
2608 ExcInternalError());
2609
2610 if (indices[d] !=
2611 numbers::invalid_dof_index)
2612 indices[d] =
2613 (indices_we_care_about.size() == 0) ?
2614 new_numbers[indices[d]] :
2615 new_numbers[indices_we_care_about
2616 .index_within_set(
2617 indices[d])];
2618 }
2619 }
2620 }
2621 }
2622 }
2623 }
2624
2625
2626
2627 template <int dim, int spacedim>
2628 static void
2629 renumber_mg_dofs(
2630 const std::vector<::types::global_dof_index> &new_numbers,
2631 const IndexSet &indices_we_care_about,
2632 DoFHandler<dim, spacedim> &dof_handler,
2633 const unsigned int level,
2634 const bool check_validity)
2635 {
2636 Assert(
2637 dof_handler.hp_capability_enabled == false,
2638 (typename DoFHandler<dim, spacedim>::ExcNotImplementedWithHP()));
2639
2640 Assert(level < dof_handler.get_triangulation().n_global_levels(),
2641 ExcInternalError());
2642
2643 // renumber DoF indices on vertices, cells, and faces. this
2644 // can be done in parallel because the respective functions
2645 // work on separate data structures
2646 Threads::TaskGroup<> tasks;
2647 tasks += Threads::new_task([&]() {
2648 renumber_vertex_mg_dofs(new_numbers,
2649 indices_we_care_about,
2650 dof_handler,
2651 level);
2652 });
2653 tasks += Threads::new_task([&]() {
2654 renumber_face_mg_dofs(new_numbers,
2655 indices_we_care_about,
2656 dof_handler,
2657 level,
2658 check_validity);
2659 });
2660 tasks += Threads::new_task([&]() {
2661 renumber_cell_mg_dofs(new_numbers,
2662 indices_we_care_about,
2663 dof_handler,
2664 level);
2665 });
2666 tasks.join_all();
2667 }
2668 };
2669
2670
2671
2672 /* --------------------- class Sequential ---------------- */
2673
2674
2675
2676 template <int dim, int spacedim>
2677 Sequential<dim, spacedim>::Sequential(
2678 DoFHandler<dim, spacedim> &dof_handler)
2679 : dof_handler(&dof_handler)
2680 {}
2681
2682
2683
2684 template <int dim, int spacedim>
2685 NumberCache
2686 Sequential<dim, spacedim>::distribute_dofs() const
2687 {
2688 const types::global_dof_index n_initial_dofs =
2689 Implementation::distribute_dofs(numbers::invalid_subdomain_id,
2690 *dof_handler);
2691
2692 const types::global_dof_index n_dofs =
2693 Implementation::unify_dof_indices(*dof_handler,
2694 n_initial_dofs,
2695 /*check_validity=*/true);
2696
2697 // return a sequential, complete index set
2698 return NumberCache(n_dofs);
2699 }
2700
2701
2702
2703 template <int dim, int spacedim>
2704 std::vector<NumberCache>
2705 Sequential<dim, spacedim>::distribute_mg_dofs() const
2706 {
2707 std::vector<NumberCache> number_caches;
2708 number_caches.reserve(dof_handler->get_triangulation().n_levels());
2709 for (unsigned int level = 0;
2710 level < dof_handler->get_triangulation().n_levels();
2711 ++level)
2712 {
2713 // first distribute dofs on this level
2714 const types::global_dof_index n_level_dofs =
2715 Implementation::distribute_dofs_on_level(
2716 numbers::invalid_subdomain_id, *dof_handler, level);
2717
2718 // then add a complete, sequential index set
2719 number_caches.emplace_back(n_level_dofs);
2720 }
2721
2722 return number_caches;
2723 }
2724
2725
2726
2727 template <int dim, int spacedim>
2728 NumberCache
2729 Sequential<dim, spacedim>::renumber_dofs(
2730 const std::vector<types::global_dof_index> &new_numbers) const
2731 {
2732 Implementation::renumber_dofs(new_numbers,
2733 IndexSet(0),
2734 *dof_handler,
2735 /*check_validity=*/true);
2736
2737 // return a sequential, complete index set. take into account that the
2738 // number of DoF indices may in fact be smaller than there were before
2739 // if some previously separately numbered dofs have been identified.
2740 // this is, for example, what we do when the DoFHandler has hp-
2741 // capabilities enabled: it first enumerates all DoFs on cells
2742 // independently, and then unifies some located at vertices or faces;
2743 // this leaves us with fewer DoFs than there were before, so use the
2744 // largest index as the one to determine the size of the index space
2745 if (new_numbers.empty())
2746 return NumberCache();
2747 else
2748 return NumberCache(
2749 *std::max_element(new_numbers.begin(), new_numbers.end()) + 1);
2750 }
2751
2752
2753
2754 template <int dim, int spacedim>
2755 NumberCache
2756 Sequential<dim, spacedim>::renumber_mg_dofs(
2757 const unsigned int level,
2758 const std::vector<types::global_dof_index> &new_numbers) const
2759 {
2760 Implementation::renumber_mg_dofs(
2761 new_numbers, IndexSet(0), *dof_handler, level, true);
2762
2763 // return a sequential, complete index set
2764 return NumberCache(new_numbers.size());
2765 }
2766
2767
2768 /* --------------------- class ParallelShared ---------------- */
2769
2770
2771 template <int dim, int spacedim>
2772 ParallelShared<dim, spacedim>::ParallelShared(
2773 DoFHandler<dim, spacedim> &dof_handler)
2774 : dof_handler(&dof_handler)
2775 {}
2776
2777
2778
2779 namespace
2780 {
2789 template <int dim, int spacedim>
2790 std::vector<types::subdomain_id>
2791 get_dof_subdomain_association(
2792 const DoFHandler<dim, spacedim> &dof_handler,
2793 const types::global_dof_index n_dofs,
2794 const unsigned int n_procs)
2795 {
2796 (void)n_procs;
2797 std::vector<types::subdomain_id> subdomain_association(
2798 n_dofs, numbers::invalid_subdomain_id);
2799 std::vector<types::global_dof_index> local_dof_indices;
2800 local_dof_indices.reserve(
2801 dof_handler.get_fe_collection().max_dofs_per_cell());
2802
2803 // loop over all cells and record which subdomain a DoF belongs to.
2804 // give to the smaller subdomain_id in case it is on an interface
2805 for (const auto &cell : dof_handler.active_cell_iterators())
2806 {
2807 // get the owner of the cell; note that we have made sure above
2808 // that all cells are either locally owned or ghosts (not
2809 // artificial), so this call will always yield the true owner;
2810 // note that the cache is not assigned yet, so we must bypass it
2811 const types::subdomain_id subdomain_id = cell->subdomain_id();
2812 const unsigned int dofs_per_cell =
2813 cell->get_fe().n_dofs_per_cell();
2814 local_dof_indices.resize(dofs_per_cell);
2815 internal::DoFAccessorImplementation::Implementation::
2816 get_dof_indices(*cell,
2817 local_dof_indices,
2818 cell->active_fe_index());
2819
2820 // set subdomain ids. if dofs already have their values set then
2821 // they must be on partition interfaces. In that case assign them
2822 // to the processor with the smaller subdomain id.
2823 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2824 if (subdomain_association[local_dof_indices[i]] ==
2825 numbers::invalid_subdomain_id)
2826 subdomain_association[local_dof_indices[i]] = subdomain_id;
2827 else if (subdomain_association[local_dof_indices[i]] >
2828 subdomain_id)
2829 {
2830 subdomain_association[local_dof_indices[i]] = subdomain_id;
2831 }
2832 }
2833
2834 Assert(std::find(subdomain_association.begin(),
2835 subdomain_association.end(),
2836 numbers::invalid_subdomain_id) ==
2837 subdomain_association.end(),
2838 ExcInternalError());
2839
2840 Assert(*std::max_element(subdomain_association.begin(),
2841 subdomain_association.end()) < n_procs,
2842 ExcInternalError());
2843
2844 return subdomain_association;
2845 }
2846
2847
2854 template <int dim, int spacedim>
2855 std::vector<types::subdomain_id>
2856 get_dof_level_subdomain_association(
2857 const DoFHandler<dim, spacedim> &dof_handler,
2858 const types::global_dof_index n_dofs_on_level,
2859 const unsigned int n_procs,
2860 const unsigned int level)
2861 {
2862 (void)n_procs;
2863 std::vector<types::subdomain_id> level_subdomain_association(
2864 n_dofs_on_level, numbers::invalid_subdomain_id);
2865 std::vector<types::global_dof_index> local_dof_indices;
2866 local_dof_indices.reserve(
2867 dof_handler.get_fe_collection().max_dofs_per_cell());
2868
2869 // loop over all cells and record which subdomain a DoF belongs to.
2870 // interface goes to processor with smaller subdomain id
2871 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
2872 {
2873 // get the owner of the cell; note that we have made sure above
2874 // that all cells are either locally owned or ghosts (not
2875 // artificial), so this call will always yield the true owner
2876 const types::subdomain_id level_subdomain_id =
2877 cell->level_subdomain_id();
2878 const unsigned int dofs_per_cell =
2879 cell->get_fe().n_dofs_per_cell();
2880 local_dof_indices.resize(dofs_per_cell);
2881 cell->get_mg_dof_indices(local_dof_indices);
2882
2883 // set level subdomain ids. if dofs already have their values set
2884 // then they must be on partition interfaces. In that case assign
2885 // them to the processor with the smaller subdomain id.
2886 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2887 if (level_subdomain_association[local_dof_indices[i]] ==
2888 numbers::invalid_subdomain_id)
2889 level_subdomain_association[local_dof_indices[i]] =
2890 level_subdomain_id;
2891 else if (level_subdomain_association[local_dof_indices[i]] >
2892 level_subdomain_id)
2893 {
2894 level_subdomain_association[local_dof_indices[i]] =
2895 level_subdomain_id;
2896 }
2897 }
2898
2899 Assert(std::find(level_subdomain_association.begin(),
2900 level_subdomain_association.end(),
2901 numbers::invalid_subdomain_id) ==
2902 level_subdomain_association.end(),
2903 ExcInternalError());
2904
2905 Assert(*std::max_element(level_subdomain_association.begin(),
2906 level_subdomain_association.end()) < n_procs,
2907 ExcInternalError());
2908
2909 return level_subdomain_association;
2910 }
2911 } // namespace
2912
2913
2914
2915 template <int dim, int spacedim>
2916 NumberCache
2917 ParallelShared<dim, spacedim>::distribute_dofs() const
2918 {
2919 const ::parallel::shared::Triangulation<dim, spacedim> *tr =
2920 (dynamic_cast<
2921 const ::parallel::shared::Triangulation<dim, spacedim> *>(
2922 &this->dof_handler->get_triangulation()));
2923 Assert(tr != nullptr, ExcInternalError());
2924
2925 const unsigned int n_procs =
2926 Utilities::MPI::n_mpi_processes(tr->get_mpi_communicator());
2927
2928 // If an underlying shared::Tria allows artificial cells, we need to
2929 // restore the true cell owners temporarily.
2930 // We use the TemporarilyRestoreSubdomainIds class for this purpose: we
2931 // save the current set of subdomain ids, set subdomain ids to the
2932 // "true" owner of each cell upon construction of the
2933 // TemporarilyRestoreSubdomainIds object, and later restore these flags
2934 // when it is destroyed.
2935 const internal::parallel::shared::
2936 TemporarilyRestoreSubdomainIds<dim, spacedim>
2937 subdomain_modifier(*tr);
2938
2939 // first let the sequential algorithm do its magic. it is going to
2940 // enumerate DoFs on all cells, regardless of owner
2941 const types::global_dof_index n_initial_dofs =
2942 Implementation::distribute_dofs(numbers::invalid_subdomain_id,
2943 *this->dof_handler);
2944
2945 const types::global_dof_index n_dofs =
2946 Implementation::unify_dof_indices(*this->dof_handler,
2947 n_initial_dofs,
2948 /*check_validity=*/true);
2949
2950 // then re-enumerate them based on their subdomain association.
2951 // for this, we first have to identify for each current DoF
2952 // index which subdomain they belong to. ideally, we would
2953 // like to call DoFRenumbering::subdomain_wise(), but
2954 // because the NumberCache of the current DoFHandler is not
2955 // fully set up yet, we can't quite do that. also, that
2956 // function has to deal with other kinds of triangulations as
2957 // well, whereas we here know what kind of triangulation
2958 // we have and can simplify the code accordingly
2959 std::vector<types::global_dof_index> new_dof_indices(
2960 n_dofs, enumeration_dof_index);
2961 {
2962 // first get the association of each dof with a subdomain and
2963 // determine the total number of subdomain ids used
2964 const std::vector<types::subdomain_id> subdomain_association =
2965 get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
2966
2967 // then renumber the subdomains by first looking at those belonging
2968 // to subdomain 0, then those of subdomain 1, etc. note that the
2969 // algorithm is stable, i.e. if two dofs i,j have i<j and belong to
2970 // the same subdomain, then they will be in this order also after
2971 // reordering
2972 types::global_dof_index next_free_index = 0;
2973 for (types::subdomain_id subdomain = 0; subdomain < n_procs;
2974 ++subdomain)
2975 for (types::global_dof_index i = 0; i < n_dofs; ++i)
2976 if (subdomain_association[i] == subdomain)
2977 {
2978 Assert(new_dof_indices[i] == enumeration_dof_index,
2979 ExcInternalError());
2980 new_dof_indices[i] = next_free_index;
2981 ++next_free_index;
2982 }
2983
2984 // we should have numbered all dofs
2985 Assert(next_free_index == n_dofs, ExcInternalError());
2986 Assert(std::find(new_dof_indices.begin(),
2987 new_dof_indices.end(),
2988 enumeration_dof_index) == new_dof_indices.end(),
2989 ExcInternalError());
2990 }
2991 // finally do the renumbering. we can use the sequential
2992 // version of the function because we do things on all
2993 // cells and all cells have their subdomain ids and DoFs
2994 // correctly set
2995 Implementation::renumber_dofs(new_dof_indices,
2996 IndexSet(0),
2997 *this->dof_handler,
2998 /*check_validity=*/true);
2999
3000 // update the number cache. for this, we first have to find the
3001 // subdomain association for each DoF again following renumbering, from
3002 // which we can then compute the IndexSets of locally owned DoFs for all
3003 // processors. all other fields then follow from this
3004 //
3005 // given the way we enumerate degrees of freedom, the locally owned
3006 // ranges must all be contiguous and consecutive. this makes filling
3007 // the IndexSets cheap. an assertion at the top verifies that this
3008 // assumption is true
3009 const std::vector<types::subdomain_id> subdomain_association =
3010 get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);
3011
3012 for (types::global_dof_index i = 1; i < n_dofs; ++i)
3013 Assert(subdomain_association[i] >= subdomain_association[i - 1],
3014 ExcInternalError());
3015
3016 std::vector<IndexSet> locally_owned_dofs_per_processor(
3017 n_procs, IndexSet(n_dofs));
3018 {
3019 // we know that the set of subdomain indices is contiguous from
3020 // the assertion above; find the start and end index for each
3021 // processor, taking into account that sometimes a processor
3022 // may not in fact have any DoFs at all. we do the latter
3023 // by just identifying contiguous ranges of subdomain_ids
3024 // and filling IndexSets for those subdomains; subdomains
3025 // that don't appear will lead to IndexSets that are simply
3026 // never touched and remain empty as initialized above.
3027 types::global_dof_index start_index = 0;
3028 types::global_dof_index end_index = 0;
3029 while (start_index < n_dofs)
3030 {
3031 while ((end_index < n_dofs) &&
3032 (subdomain_association[end_index] ==
3033 subdomain_association[start_index]))
3034 ++end_index;
3035
3036 // we've now identified a range of same indices. set that
3037 // range in the corresponding IndexSet
3038 if (end_index > start_index)
3039 {
3040 const types::subdomain_id subdomain_owner =
3041 subdomain_association[start_index];
3042 locally_owned_dofs_per_processor[subdomain_owner].add_range(
3043 start_index, end_index);
3044 }
3045
3046 // then move on to thinking about the next range
3047 start_index = end_index;
3048 }
3049 }
3050
3051 // return a NumberCache object made up from the sets of locally
3052 // owned DoFs
3053 return NumberCache(
3054 locally_owned_dofs_per_processor,
3055 this->dof_handler->get_triangulation().locally_owned_subdomain());
3056 }
3057
3058
3059
3060 template <int dim, int spacedim>
3061 std::vector<NumberCache>
3062 ParallelShared<dim, spacedim>::distribute_mg_dofs() const
3063 {
3064 const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3065 (dynamic_cast<
3066 const ::parallel::shared::Triangulation<dim, spacedim> *>(
3067 &this->dof_handler->get_triangulation()));
3068 Assert(tr != nullptr, ExcInternalError());
3069
3070 AssertThrow((tr->is_multilevel_hierarchy_constructed()),
3071 ExcMessage(
3072 "Multigrid DoFs can only be distributed on a parallel "
3073 "Triangulation if the flag construct_multigrid_hierarchy "
3074 "is set in the constructor."));
3075
3076 const unsigned int n_procs =
3077 Utilities::MPI::n_mpi_processes(tr->get_mpi_communicator());
3078 const unsigned int n_levels = tr->n_global_levels();
3079
3080 std::vector<NumberCache> number_caches;
3081 number_caches.reserve(n_levels);
3082
3083 // We create an index set for each level
3084 for (unsigned int lvl = 0; lvl < n_levels; ++lvl)
3085 {
3086 // If the underlying shared::Tria allows artificial cells,
3087 // then save the current set of level subdomain ids, and set
3088 // subdomain ids to the "true" owner of each cell. we later
3089 // restore these flags
3090 // Note: "allows_artificial_cells" is currently enforced for
3091 // MG computations.
3092 std::vector<types::subdomain_id> saved_level_subdomain_ids;
3093 saved_level_subdomain_ids.resize(tr->n_cells(lvl));
3094 {
3095 typename ::parallel::shared::Triangulation<dim, spacedim>::
3096 cell_iterator cell =
3097 this->dof_handler->get_triangulation().begin(
3098 lvl),
3099 endc =
3100 this->dof_handler->get_triangulation().end(lvl);
3101
3102 const std::vector<types::subdomain_id> &true_level_subdomain_ids =
3103 tr->get_true_level_subdomain_ids_of_cells(lvl);
3104
3105 for (unsigned int index = 0; cell != endc; ++cell, ++index)
3106 {
3107 saved_level_subdomain_ids[index] = cell->level_subdomain_id();
3108 cell->set_level_subdomain_id(true_level_subdomain_ids[index]);
3109 }
3110 }
3111
3112 // Next let the sequential algorithm do its magic. it is going to
3113 // enumerate DoFs on all cells on the given level, regardless of
3114 // owner
3115 const types::global_dof_index n_dofs_on_level =
3116 Implementation::distribute_dofs_on_level(
3117 numbers::invalid_subdomain_id, *this->dof_handler, lvl);
3118
3119 // then re-enumerate them based on their level subdomain
3120 // association. for this, we first have to identify for each current
3121 // DoF index which subdomain they belong to. ideally, we would like
3122 // to call DoFRenumbering::subdomain_wise(), but because the
3123 // NumberCache of the current DoFHandler is not fully set up yet, we
3124 // can't quite do that. also, that function has to deal with other
3125 // kinds of triangulations as well, whereas we here know what kind
3126 // of triangulation we have and can simplify the code accordingly
3127 std::vector<types::global_dof_index> new_dof_indices(
3128 n_dofs_on_level, numbers::invalid_dof_index);
3129 {
3130 // first get the association of each dof with a subdomain and
3131 // determine the total number of subdomain ids used
3132 const std::vector<types::subdomain_id>
3133 level_subdomain_association =
3134 get_dof_level_subdomain_association(*this->dof_handler,
3135 n_dofs_on_level,
3136 n_procs,
3137 lvl);
3138
3139 // then renumber the subdomains by first looking at those
3140 // belonging to subdomain 0, then those of subdomain 1, etc. note
3141 // that the algorithm is stable, i.e. if two dofs i,j have i<j and
3142 // belong to the same subdomain, then they will be in this order
3143 // also after reordering
3144 types::global_dof_index next_free_index = 0;
3145 for (types::subdomain_id level_subdomain = 0;
3146 level_subdomain < n_procs;
3147 ++level_subdomain)
3148 for (types::global_dof_index i = 0; i < n_dofs_on_level; ++i)
3149 if (level_subdomain_association[i] == level_subdomain)
3150 {
3151 Assert(new_dof_indices[i] == numbers::invalid_dof_index,
3152 ExcInternalError());
3153 new_dof_indices[i] = next_free_index;
3154 ++next_free_index;
3155 }
3156
3157 // we should have numbered all dofs
3158 Assert(next_free_index == n_dofs_on_level, ExcInternalError());
3159 Assert(std::find(new_dof_indices.begin(),
3160 new_dof_indices.end(),
3161 numbers::invalid_dof_index) ==
3162 new_dof_indices.end(),
3163 ExcInternalError());
3164 }
3165
3166 // finally do the renumbering. we can use the sequential
3167 // version of the function because we do things on all
3168 // cells and all cells have their subdomain ids and DoFs
3169 // correctly set
3170 Implementation::renumber_mg_dofs(
3171 new_dof_indices, IndexSet(0), *this->dof_handler, lvl, true);
3172
3173 // update the number cache. for this, we first have to find the
3174 // level subdomain association for each DoF again following
3175 // renumbering, from which we can then compute the IndexSets of
3176 // locally owned DoFs for all processors. all other fields then
3177 // follow from this
3178 //
3179 // given the way we enumerate degrees of freedom, the locally owned
3180 // ranges must all be contiguous and consecutive. this makes filling
3181 // the IndexSets cheap. an assertion at the top verifies that this
3182 // assumption is true
3183 const std::vector<types::subdomain_id> level_subdomain_association =
3184 get_dof_level_subdomain_association(*this->dof_handler,
3185 n_dofs_on_level,
3186 n_procs,
3187 lvl);
3188
3189 for (types::global_dof_index i = 1; i < n_dofs_on_level; ++i)
3190 Assert(level_subdomain_association[i] >=
3191 level_subdomain_association[i - 1],
3192 ExcInternalError());
3193
3194 std::vector<IndexSet> locally_owned_dofs_per_processor(
3195 n_procs, IndexSet(n_dofs_on_level));
3196 {
3197 // we know that the set of subdomain indices is contiguous from
3198 // the assertion above; find the start and end index for each
3199 // processor, taking into account that sometimes a processor
3200 // may not in fact have any DoFs at all. we do the latter
3201 // by just identifying contiguous ranges of level_subdomain_ids
3202 // and filling IndexSets for those subdomains; subdomains
3203 // that don't appear will lead to IndexSets that are simply
3204 // never touched and remain empty as initialized above.
3205 unsigned int start_index = 0;
3206 unsigned int end_index = 0;
3207 while (start_index < n_dofs_on_level)
3208 {
3209 while ((end_index) < n_dofs_on_level &&
3210 (level_subdomain_association[end_index] ==
3211 level_subdomain_association[start_index]))
3212 ++end_index;
3213
3214 // we've now identified a range of same indices. set that
3215 // range in the corresponding IndexSet
3216 if (end_index > start_index)
3217 {
3218 const unsigned int level_subdomain_owner =
3219 level_subdomain_association[start_index];
3220 locally_owned_dofs_per_processor[level_subdomain_owner]
3221 .add_range(start_index, end_index);
3222 }
3223
3224 // then move on to thinking about the next range
3225 start_index = end_index;
3226 }
3227 }
3228
3229 // finally, restore current level subdomain ids
3230 {
3231 typename ::parallel::shared::Triangulation<dim, spacedim>::
3232 cell_iterator cell =
3233 this->dof_handler->get_triangulation().begin(
3234 lvl),
3235 endc =
3236 this->dof_handler->get_triangulation().end(lvl);
3237
3238 for (unsigned int index = 0; cell != endc; ++cell, ++index)
3239 cell->set_level_subdomain_id(saved_level_subdomain_ids[index]);
3240
3241 // add NumberCache for current level
3242 number_caches.emplace_back(
3243 NumberCache(locally_owned_dofs_per_processor,
3244 this->dof_handler->get_triangulation()
3245 .locally_owned_subdomain()));
3246 }
3247 }
3248
3249 return number_caches;
3250 }
3251
3252
3253
3254 template <int dim, int spacedim>
3255 NumberCache
3256 ParallelShared<dim, spacedim>::renumber_dofs(
3257 const std::vector<types::global_dof_index> &new_numbers) const
3258 {
3259#ifndef DEAL_II_WITH_MPI
3260 (void)new_numbers;
3261 DEAL_II_NOT_IMPLEMENTED();
3262 return NumberCache();
3263#else
3264 // Similar to distribute_dofs() we need to have a special treatment in
3265 // case artificial cells are present.
3266 const ::parallel::shared::Triangulation<dim, spacedim> *tr =
3267 (dynamic_cast<
3268 const ::parallel::shared::Triangulation<dim, spacedim> *>(
3269 &this->dof_handler->get_triangulation()));
3270 Assert(tr != nullptr, ExcInternalError());
3271
3272 // Set subdomain IDs to the "true" owner of each cell.
3273 const internal::parallel::shared::
3274 TemporarilyRestoreSubdomainIds<dim, spacedim>
3275 subdomain_modifier(*tr);
3276
3277 std::vector<types::global_dof_index> global_gathered_numbers(
3278 this->dof_handler->n_dofs(), 0);
3279 // as we call DoFRenumbering::subdomain_wise(*dof_handler) from
3280 // distribute_dofs(), we need to support sequential-like input.
3281 // Distributed-like input from, for example, component_wise renumbering
3282 // is also supported.
3283 const bool uses_sequential_numbering =
3284 new_numbers.size() == this->dof_handler->n_dofs();
3285 const bool all_use_sequential_numbering =
3286 Utilities::MPI::logical_and(uses_sequential_numbering,
3287 tr->get_mpi_communicator());
3288 if (all_use_sequential_numbering)
3289 {
3290 global_gathered_numbers = new_numbers;
3291 }
3292 else
3293 {
3294 Assert(new_numbers.size() ==
3295 this->dof_handler->locally_owned_dofs().n_elements(),
3296 ExcInternalError());
3297 const unsigned int n_cpu =
3298 Utilities::MPI::n_mpi_processes(tr->get_mpi_communicator());
3299 std::vector<types::global_dof_index> gathered_new_numbers(
3300 this->dof_handler->n_dofs(), 0);
3301 Assert(Utilities::MPI::this_mpi_process(
3302 tr->get_mpi_communicator()) ==
3303 this->dof_handler->get_triangulation()
3304 .locally_owned_subdomain(),
3305 ExcInternalError());
3306
3307 // gather new numbers among processors into one vector
3308 {
3309 std::vector<types::global_dof_index> new_numbers_copy(
3310 new_numbers);
3311
3312 // store the number of elements that are to be received from each
3313 // process
3314 std::vector<int> rcounts(n_cpu);
3315
3316 types::global_dof_index shift = 0;
3317 // set rcounts based on new_numbers:
3318 int cur_count = new_numbers_copy.size();
3319 int ierr = MPI_Allgather(&cur_count,
3320 1,
3321 MPI_INT,
3322 rcounts.data(),
3323 1,
3324 MPI_INT,
3325 tr->get_mpi_communicator());
3326 AssertThrowMPI(ierr);
3327
3328 // compute the displacements (relative to recvbuf)
3329 // at which to place the incoming data from process i
3330 std::vector<int> displacements(n_cpu);
3331 for (unsigned int i = 0; i < n_cpu; ++i)
3332 {
3333 displacements[i] = shift;
3334 shift += rcounts[i];
3335 }
3336 Assert(new_numbers_copy.size() ==
3337 static_cast<unsigned int>(
3338 rcounts[Utilities::MPI::this_mpi_process(
3339 tr->get_mpi_communicator())]),
3340 ExcInternalError());
3341 ierr = MPI_Allgatherv(
3342 new_numbers_copy.data(),
3343 new_numbers_copy.size(),
3344 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
3345 gathered_new_numbers.data(),
3346 rcounts.data(),
3347 displacements.data(),
3348 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
3349 tr->get_mpi_communicator());
3350 AssertThrowMPI(ierr);
3351 }
3352
3353 // put new numbers according to the current
3354 // locally_owned_dofs_per_processor IndexSets
3355 types::global_dof_index shift = 0;
3356 // Make sure we have one-to-one relation between old and new DoFs by
3357 // counting:
3358 std::vector<unsigned int> count_rename_from(
3359 this->dof_handler->n_dofs(), 0);
3360 std::vector<unsigned int> count_rename_to(
3361 this->dof_handler->n_dofs(), 0);
3362
3363 const std::vector<IndexSet> locally_owned_dofs_per_processor =
3364 Utilities::MPI::all_gather(
3365 tr->get_mpi_communicator(),
3366 this->dof_handler->locally_owned_dofs());
3367
3368 for (unsigned int i = 0; i < n_cpu; ++i)
3369 {
3370 const IndexSet &iset = locally_owned_dofs_per_processor[i];
3371 for (types::global_dof_index ind = 0; ind < iset.n_elements();
3372 ind++)
3373 {
3374 const types::global_dof_index target =
3375 iset.nth_index_in_set(ind);
3376 const types::global_dof_index value =
3377 gathered_new_numbers[shift + ind];
3378 Assert(target < this->dof_handler->n_dofs(),
3379 ExcInternalError());
3380 Assert(value < this->dof_handler->n_dofs(),
3381 ExcInternalError());
3382 global_gathered_numbers[target] = value;
3383 count_rename_to[target]++;
3384 count_rename_from[value]++;
3385 }
3386 shift += iset.n_elements();
3387 }
3388
3389 Assert(*std::max_element(count_rename_from.begin(),
3390 count_rename_from.end()) == 1,
3391 ExcInternalError());
3392 Assert(*std::min_element(count_rename_from.begin(),
3393 count_rename_from.end()) == 1,
3394 ExcInternalError());
3395 Assert((*std::max_element(count_rename_to.begin(),
3396 count_rename_to.end())) == 1,
3397 ExcInternalError());
3398 Assert((*std::min_element(count_rename_to.begin(),
3399 count_rename_to.end())) == 1,
3400 ExcInternalError());
3401 }
3402
3403 // let the sequential algorithm do its magic; ignore the
3404 // return type, but reconstruct the number cache based on
3405 // which DoFs each process owns
3406 Implementation::renumber_dofs(global_gathered_numbers,
3407 IndexSet(0),
3408 *this->dof_handler,
3409 /*check_validity=*/true);
3410
3411 const NumberCache number_cache(
3412 DoFTools::locally_owned_dofs_per_subdomain(*this->dof_handler),
3413 this->dof_handler->get_triangulation().locally_owned_subdomain());
3414
3415 return number_cache;
3416#endif
3417 }
3418
3419
3420
3421 template <int dim, int spacedim>
3422 NumberCache
3423 ParallelShared<dim, spacedim>::renumber_mg_dofs(
3424 const unsigned int /*level*/,
3425 const std::vector<types::global_dof_index> & /*new_numbers*/) const
3426 {
3427 // multigrid is not currently implemented for shared triangulations
3428 DEAL_II_NOT_IMPLEMENTED();
3429
3430 return {};
3431 }
3432
3433
3434
3435 /* --------------------- class ParallelDistributed ---------------- */
3436
3437#ifdef DEAL_II_WITH_MPI
3438
3439 namespace
3440 {
3441 template <int dim, int spacedim>
3442 void
3443 communicate_mg_ghost_cells(DoFHandler<dim, spacedim> &dof_handler,
3444 std::vector<std::vector<bool>> &cell_marked)
3445 {
3446 const auto pack = [](const auto &cell) {
3447 // why would somebody request a cell that is not ours?
3448 Assert(cell->is_locally_owned_on_level(), ExcInternalError());
3449
3450 std::vector<::types::global_dof_index> data(
3451 cell->get_fe().n_dofs_per_cell());
3452 cell->get_mg_dof_indices(data);
3453
3454 return data;
3455 };
3456
3457 const auto unpack = [&cell_marked](const auto &cell,
3458 const auto &dofs) {
3459 Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3460 ExcInternalError());
3461
3462 Assert(cell->level_subdomain_id() !=
3463 ::numbers::artificial_subdomain_id,
3464 ExcInternalError());
3465
3466 bool complete = true;
3467 DoFAccessorImplementation::Implementation::process_dof_indices(
3468 *cell,
3469 dofs,
3470 0,
3471 DoFAccessorImplementation::Implementation::
3472 MGDoFIndexProcessor<dim, spacedim>(cell->level()),
3473
3474 // Intel ICC 18 and earlier for some reason believe that
3475 // numbers::invalid_dof_index is not a valid object
3476 // inside the lambda function. Fix this by creating a
3477 // local variable initialized by the global one.
3478 //
3479 // Intel ICC 19 and earlier have trouble with our Assert
3480 // macros inside the lambda function. We disable the macro
3481 // for these compilers.
3482 [&complete, invalid_dof_index = numbers::invalid_dof_index](
3483 auto &stored_index, auto received_index) {
3484 if (*received_index != invalid_dof_index)
3485 {
3486# if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
3487 Assert((stored_index == invalid_dof_index) ||
3488 (stored_index == *received_index),
3489 ExcInternalError());
3490# endif
3491 stored_index = *received_index;
3492 }
3493 else
3494 complete = false;
3495 },
3496 true);
3497
3498 if (!complete)
3499 {
3500 // We should have the cell already marked
3501 Assert(cell_marked[cell->level()][cell->index()],
3502 ExcInternalError());
3503 }
3504 else
3505 cell_marked[cell->level()][cell->index()] = false;
3506 };
3507
3508 const auto filter = [&cell_marked](const auto &cell) {
3509 return cell_marked[cell->level()][cell->index()];
3510 };
3511
3512 GridTools::exchange_cell_data_to_level_ghosts<
3513 std::vector<types::global_dof_index>,
3514 DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter);
3515 }
3516
3517
3518
3537 template <int dim, int spacedim>
3538 void
3539 communicate_dof_indices_on_marked_cells(
3540 const DoFHandler<dim, spacedim> &dof_handler,
3541 std::vector<bool> &cell_marked)
3542 {
3543# ifndef DEAL_II_WITH_MPI
3544 (void)dof_handler;
3545 DEAL_II_NOT_IMPLEMENTED();
3546# else
3547
3548 // define functions that pack data on cells that are ghost cells
3549 // somewhere else, and unpack data on cells where we get information
3550 // from elsewhere
3551 const auto pack = [](const auto &cell) {
3552 Assert(cell->is_locally_owned(), ExcInternalError());
3553
3554 std::vector<::types::global_dof_index> data(
3555 cell->get_fe().n_dofs_per_cell());
3556
3557 // bypass the cache which is not filled yet
3558 internal::DoFAccessorImplementation::Implementation::
3559 get_dof_indices(*cell, data, cell->active_fe_index());
3560
3561 return data;
3562 };
3563
3564 const auto unpack = [&cell_marked](const auto &cell,
3565 const auto &dofs) {
3566 Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3567 ExcInternalError());
3568
3569 Assert(cell->is_ghost(), ExcInternalError());
3570
3571 // Use a combined read/set function on the entities of the dof
3572 // indices to speed things up against get_dof_indices +
3573 // set_dof_indices
3574 bool complete = true;
3575 DoFAccessorImplementation::Implementation::process_dof_indices(
3576 *cell,
3577 dofs,
3578 cell->active_fe_index(),
3579 DoFAccessorImplementation::Implementation::
3580 DoFIndexProcessor<dim, spacedim>(),
3581
3582 // Intel ICC 18 and earlier for some reason believe that
3583 // numbers::invalid_dof_index is not a valid object
3584 // inside the lambda function. Fix this by creating a
3585 // local variable initialized by the global one.
3586 //
3587 // Intel ICC 19 and earlier have trouble with our Assert
3588 // macros inside the lambda function. We disable the macro
3589 // for these compilers.
3590 [&complete, invalid_dof_index = numbers::invalid_dof_index](
3591 auto &stored_index, const auto received_index) {
3592 if (*received_index != invalid_dof_index)
3593 {
3594# if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
3595 Assert((stored_index == invalid_dof_index) ||
3596 (stored_index == *received_index),
3597 ExcInternalError());
3598# endif
3599 stored_index = *received_index;
3600 }
3601 else
3602 complete = false;
3603 },
3604 false);
3605
3606 if (!complete)
3607 {
3608 // We should have the cell already marked
3609 Assert(cell_marked[cell->active_cell_index()],
3610 ExcInternalError());
3611 }
3612 else
3613 cell_marked[cell->active_cell_index()] = false;
3614 };
3615
3616 const auto filter = [&cell_marked](const auto &cell) {
3617 return cell_marked[cell->active_cell_index()];
3618 };
3619
3620 GridTools::exchange_cell_data_to_ghosts<
3621 std::vector<types::global_dof_index>,
3622 DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter);
3623# endif
3624 }
3625
3626
3627
3628 } // namespace
3629
3630#endif // DEAL_II_WITH_MPI
3631
3632
3633
3634 template <int dim, int spacedim>
3635 ParallelDistributed<dim, spacedim>::ParallelDistributed(
3636 DoFHandler<dim, spacedim> &dof_handler)
3637 : dof_handler(&dof_handler)
3638 {}
3639
3640
3641
3642 template <int dim, int spacedim>
3643 NumberCache
3644 ParallelDistributed<dim, spacedim>::distribute_dofs() const
3645 {
3646#ifndef DEAL_II_WITH_MPI
3647 DEAL_II_NOT_IMPLEMENTED();
3648 return NumberCache();
3649#else
3650
3651 ::parallel::DistributedTriangulationBase<dim, spacedim>
3652 *triangulation =
3653 (dynamic_cast<
3654 ::parallel::DistributedTriangulationBase<dim, spacedim> *>(
3655 const_cast<::Triangulation<dim, spacedim> *>(
3656 &dof_handler->get_triangulation())));
3657 Assert(triangulation != nullptr, ExcInternalError());
3658
3659 const types::subdomain_id subdomain_id =
3660 triangulation->locally_owned_subdomain();
3661
3662
3663 /*
3664 The following algorithm has a number of stages that are all
3665 documented in the paper that describes the parallel::distributed
3666 functionality:
3667
3668 1/ locally enumerate dofs on locally owned cells
3669 2/ eliminate dof duplicates on all cells.
3670 un-numerate those that are on interfaces with ghost
3671 cells and that we don't own based on the tie-breaking
3672 criterion. unify dofs afterwards.
3673 3/ unify dofs and re-enumerate the remaining valid ones.
3674 the end result is that we only enumerate locally owned
3675 DoFs
3676 4/ shift indices so that each processor has a unique
3677 range of indices
3678 5/ for all locally owned cells that are ghost
3679 cells somewhere else, send our own DoF indices
3680 to the appropriate set of other processors.
3681 overwrite invalid DoF indices on ghost interfaces
3682 with the corresponding valid ones that we now know.
3683 6/ send DoF indices again to get the correct indices
3684 on ghost cells that we may not have known earlier
3685 */
3686
3687 // --------- Phase 1: enumerate dofs on locally owned cells
3688 const types::global_dof_index n_initial_local_dofs =
3689 Implementation::distribute_dofs(subdomain_id, *dof_handler);
3690
3691 // --------- Phase 2: eliminate dof duplicates on all cells:
3692 // - un-numerate dofs on interfaces to ghost cells
3693 // that we don't own
3694 // - in case of hp-support, unify dofs
3695 std::vector<::types::global_dof_index> renumbering(
3696 n_initial_local_dofs, enumeration_dof_index);
3697
3698 // first, we invalidate degrees of freedom that belong to processors
3699 // of a lower rank, from which we will receive the final (and lower)
3700 // degrees of freedom later.
3701 Implementation::
3702 invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering(
3703 renumbering, subdomain_id, *dof_handler);
3704
3705 // then, we identify DoF duplicates if the DoFHandler has hp-
3706 // capabilities
3707 std::vector<std::map<types::global_dof_index, types::global_dof_index>>
3708 all_constrained_indices(dim);
3709 Implementation::compute_dof_identities(all_constrained_indices,
3710 *dof_handler);
3711
3712 // --------- Phase 3: re-enumerate the valid degrees of freedom
3713 // consecutively. thus, we finally receive the
3714 // correct number of locally owned DoFs after
3715 // this step.
3716 //
3717 // the order in which we handle Phases 2 and 3 is important,
3718 // since we want to clarify ownership of degrees of freedom before
3719 // we actually unify and enumerate their indices. otherwise, we could
3720 // end up having a degree of freedom to which only invalid indices will
3721 // be assigned.
3722 types::global_dof_index n_identity_constrained_indices = 0;
3723 for (const auto &constrained_indices : all_constrained_indices)
3724 for (const auto &index : constrained_indices)
3725 if (renumbering[index.first] != numbers::invalid_dof_index)
3726 ++n_identity_constrained_indices;
3727
3728 const types::global_dof_index n_locally_owned_dofs =
3729 std::count(renumbering.begin(),
3730 renumbering.end(),
3731 enumeration_dof_index) -
3732 n_identity_constrained_indices;
3733
3734 // --------- Phase 4: shift indices so that each processor has a unique
3735 // range of indices
3736 const auto [my_shift, n_global_dofs] =
3737 Utilities::MPI::partial_and_total_sum(
3738 n_locally_owned_dofs, triangulation->get_mpi_communicator());
3739
3740
3741 // make dof indices globally consecutive
3742 Implementation::enumerate_dof_indices_for_renumbering(
3743 renumbering, all_constrained_indices, my_shift);
3744
3745 // now re-enumerate all dofs to this shifted and condensed
3746 // numbering form. we renumber some dofs as invalid, so
3747 // choose the nocheck-version.
3748 Implementation::renumber_dofs(renumbering,
3749 IndexSet(0),
3750 *dof_handler,
3751 /*check_validity=*/false);
3752
3753 NumberCache number_cache;
3754 number_cache.n_global_dofs = n_global_dofs;
3755 number_cache.n_locally_owned_dofs = n_locally_owned_dofs;
3756 number_cache.locally_owned_dofs = IndexSet(n_global_dofs);
3757 number_cache.locally_owned_dofs.add_range(my_shift,
3758 my_shift +
3759 n_locally_owned_dofs);
3760 number_cache.locally_owned_dofs.compress();
3761
3762 // this ends the phase where we enumerate degrees of freedom on
3763 // each processor. what is missing is communicating DoF indices
3764 // on ghost cells
3765
3766 // --------- Phase 5: for all locally owned cells that are ghost
3767 // cells somewhere else, send our own DoF indices
3768 // to the appropriate set of other processors
3769 {
3770 // mark all cells that either have to send data (locally
3771 // owned cells that are adjacent to ghost neighbors in some
3772 // way) or receive data (all ghost cells) via the user flags
3773 std::vector<bool> cell_marked(triangulation->n_active_cells());
3774 for (const auto &cell : dof_handler->active_cell_iterators())
3775 if (cell->is_ghost())
3776 cell_marked[cell->active_cell_index()] = true;
3777
3778 // Send and receive cells. After this, only the local cells
3779 // are marked, that received new data. This has to be
3780 // communicated in a second communication step.
3781 //
3782 // as explained in the 'distributed' paper, this has to be
3783 // done twice
3784 communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
3785
3786 // If the DoFHandler has hp-capabilities enabled, then we may have
3787 // received valid indices of degrees of freedom that are dominated
3788 // by a FE object adjacent to a ghost interface. thus, we overwrite
3789 // the remaining invalid indices with the valid ones in this step.
3790 Implementation::merge_invalid_dof_indices_on_ghost_interfaces(
3791 *dof_handler);
3792
3793 // --------- Phase 6: all locally owned cells have their correct
3794 // DoF indices set. however, some ghost cells
3795 // may still have invalid ones. thus, exchange
3796 // one more time.
3797 communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
3798
3799 // at this point, we must have taken care of the data transfer
3800 // on all cells we had previously marked. verify this
3801 if constexpr (running_in_debug_mode())
3802 {
3803 for (const auto &cell : dof_handler->active_cell_iterators())
3804 Assert(cell_marked[cell->active_cell_index()] == false,
3805 ExcInternalError());
3806 }
3807 }
3808
3809 if constexpr (running_in_debug_mode())
3810 {
3811 // check that we are really done
3812 {
3813 std::vector<::types::global_dof_index> local_dof_indices;
3814
3815 for (const auto &cell : dof_handler->active_cell_iterators())
3816 if (!cell->is_artificial())
3817 {
3818 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
3819 cell->get_dof_indices(local_dof_indices);
3820 if (local_dof_indices.end() !=
3821 std::find(local_dof_indices.begin(),
3822 local_dof_indices.end(),
3823 numbers::invalid_dof_index))
3824 {
3825 if (cell->is_ghost())
3826 {
3827 Assert(false,
3828 ExcMessage(
3829 "A ghost cell ended up with incomplete "
3830 "DoF index information. This should not "
3831 "have happened!"));
3832 }
3833 else
3834 {
3835 Assert(
3836 false,
3837 ExcMessage(
3838 "A locally owned cell ended up with incomplete "
3839 "DoF index information. This should not "
3840 "have happened!"));
3841 }
3842 }
3843 }
3844 }
3845 } // DEBUG
3846 return number_cache;
3847#endif // DEAL_II_WITH_MPI
3848 }
3849
3850
3851
3852 template <int dim, int spacedim>
3853 std::vector<NumberCache>
3854 ParallelDistributed<dim, spacedim>::distribute_mg_dofs() const
3855 {
3856#ifndef DEAL_II_WITH_MPI
3857 DEAL_II_NOT_IMPLEMENTED();
3858 return std::vector<NumberCache>();
3859#else
3860
3861 ::parallel::DistributedTriangulationBase<dim, spacedim>
3862 *triangulation =
3863 (dynamic_cast<
3864 ::parallel::DistributedTriangulationBase<dim, spacedim> *>(
3865 const_cast<::Triangulation<dim, spacedim> *>(
3866 &dof_handler->get_triangulation())));
3867 Assert(triangulation != nullptr, ExcInternalError());
3868
3869 AssertThrow((triangulation->is_multilevel_hierarchy_constructed()),
3870 ExcMessage(
3871 "Multigrid DoFs can only be distributed on a parallel "
3872 "Triangulation if the flag construct_multigrid_hierarchy "
3873 "is set in the constructor."));
3874
3875 // loop over all levels that exist globally (across all
3876 // processors), even if the current processor does not in fact
3877 // have any cells on that level or if the local part of the
3878 // Triangulation has fewer levels. we need to do this because
3879 // we need to communicate across all processors on all levels
3880 const unsigned int n_levels = triangulation->n_global_levels();
3881 std::vector<NumberCache> number_caches;
3882 number_caches.reserve(n_levels);
3883 for (unsigned int level = 0; level < n_levels; ++level)
3884 {
3885 NumberCache level_number_cache;
3886
3887 //* 1. distribute on own subdomain
3888 const unsigned int n_initial_local_dofs =
3889 Implementation::distribute_dofs_on_level(
3890 triangulation->locally_owned_subdomain(), *dof_handler, level);
3891
3892 //* 2. iterate over ghostcells and kill dofs that are not
3893 // owned by us, which we mark by invalid_dof_index
3894 std::vector<::types::global_dof_index> renumbering(
3895 n_initial_local_dofs, enumeration_dof_index);
3896
3897 if (level < triangulation->n_levels())
3898 {
3899 std::vector<::types::global_dof_index> local_dof_indices;
3900
3901 for (const auto &cell :
3902 dof_handler->cell_iterators_on_level(level))
3903 if (cell->level_subdomain_id() !=
3904 numbers::artificial_subdomain_id &&
3905 (cell->level_subdomain_id() <
3906 triangulation->locally_owned_subdomain()))
3907 {
3908 // we found a neighboring ghost cell whose
3909 // subdomain is "stronger" than our own
3910 // subdomain
3911
3912 // delete all dofs that live there and that we
3913 // have previously assigned a number to
3914 // (i.e. the ones on the interface)
3915 local_dof_indices.resize(
3916 cell->get_fe().n_dofs_per_cell());
3917 cell->get_mg_dof_indices(local_dof_indices);
3918 for (unsigned int i = 0;
3919 i < cell->get_fe().n_dofs_per_cell();
3920 ++i)
3921 if (local_dof_indices[i] != numbers::invalid_dof_index)
3922 renumbering[local_dof_indices[i]] =
3923 numbers::invalid_dof_index;
3924 }
3925 }
3926
3927 level_number_cache.n_locally_owned_dofs =
3928 std::count(renumbering.begin(),
3929 renumbering.end(),
3930 enumeration_dof_index);
3931
3932 //* 3. communicate local dofcount and shift ids to make
3933 // them unique
3934 const auto [my_shift, n_global_dofs] =
3935 Utilities::MPI::partial_and_total_sum(
3936 level_number_cache.n_locally_owned_dofs,
3937 triangulation->get_mpi_communicator());
3938 level_number_cache.n_global_dofs = n_global_dofs;
3939
3940 // assign appropriate indices
3941 types::global_dof_index next_free_index = my_shift;
3942 for (types::global_dof_index &index : renumbering)
3943 if (index == enumeration_dof_index)
3944 index = next_free_index++;
3945
3946 // now re-enumerate all dofs to this shifted and condensed
3947 // numbering form. we renumber some dofs as invalid, so
3948 // choose the nocheck-version of the function
3949 //
3950 // of course there is nothing for us to renumber if the
3951 // level we are currently dealing with doesn't even exist
3952 // within the current triangulation, so skip renumbering
3953 // in that case
3954 if (level < triangulation->n_levels())
3955 Implementation::renumber_mg_dofs(
3956 renumbering, IndexSet(0), *dof_handler, level, false);
3957
3958 // now a little bit of housekeeping
3959 level_number_cache.locally_owned_dofs =
3960 IndexSet(level_number_cache.n_global_dofs);
3961 level_number_cache.locally_owned_dofs.add_range(
3962 next_free_index - level_number_cache.n_locally_owned_dofs,
3963 next_free_index);
3964 level_number_cache.locally_owned_dofs.compress();
3965
3966 number_caches.emplace_back(level_number_cache);
3967 }
3968
3969
3970 //* communicate ghost DoFs
3971 // We mark all ghost cells by setting the user_flag and then request
3972 // these cells from the corresponding owners. As this information
3973 // can be incomplete,
3974 {
3975 std::vector<std::vector<bool>> cell_marked(triangulation->n_levels());
3976 for (unsigned int l = 0; l < triangulation->n_levels(); ++l)
3977 cell_marked[l].resize(triangulation->n_raw_cells(l));
3978 for (const auto &cell : dof_handler->cell_iterators())
3979 if (cell->is_ghost_on_level())
3980 cell_marked[cell->level()][cell->index()] = true;
3981
3982 // Phase 1. Request all marked cells from corresponding owners. If we
3983 // managed to get every DoF, remove the user_flag, otherwise we
3984 // will request them again in the step below.
3985 communicate_mg_ghost_cells(*dof_handler, cell_marked);
3986
3987 // Phase 2, only request the cells that were not completed
3988 // in Phase 1.
3989 communicate_mg_ghost_cells(*dof_handler, cell_marked);
3990
3991 if constexpr (running_in_debug_mode())
3992 {
3993 // make sure we have finished all cells:
3994 for (const auto &cell : dof_handler->cell_iterators())
3995 Assert(cell_marked[cell->level()][cell->index()] == false,
3996 ExcInternalError());
3997 }
3998 }
3999
4000
4001
4002 if constexpr (running_in_debug_mode())
4003 {
4004 // check that we are really done
4005 {
4006 std::vector<::types::global_dof_index> local_dof_indices;
4007 for (const auto &cell : dof_handler->cell_iterators())
4008 if (cell->level_subdomain_id() !=
4009 ::numbers::artificial_subdomain_id)
4010 {
4011 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
4012 cell->get_mg_dof_indices(local_dof_indices);
4013 if (local_dof_indices.end() !=
4014 std::find(local_dof_indices.begin(),
4015 local_dof_indices.end(),
4016 numbers::invalid_dof_index))
4017 {
4018 Assert(false,
4019 ExcMessage("not all DoFs got distributed!"));
4020 }
4021 }
4022 }
4023 } // DEBUG
4024
4025 return number_caches;
4026
4027#endif // DEAL_II_WITH_MPI
4028 }
4029
4030
4031 template <int dim, int spacedim>
4032 NumberCache
4033 ParallelDistributed<dim, spacedim>::renumber_dofs(
4034 const std::vector<::types::global_dof_index> &new_numbers) const
4035 {
4036 (void)new_numbers;
4037
4038 Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
4039 ExcInternalError());
4040
4041#ifndef DEAL_II_WITH_MPI
4042 DEAL_II_NOT_IMPLEMENTED();
4043 return NumberCache();
4044#else
4045
4046 ::parallel::DistributedTriangulationBase<dim, spacedim>
4047 *triangulation =
4048 (dynamic_cast<
4049 ::parallel::DistributedTriangulationBase<dim, spacedim> *>(
4050 const_cast<::Triangulation<dim, spacedim> *>(
4051 &dof_handler->get_triangulation())));
4052 Assert(triangulation != nullptr, ExcInternalError());
4053
4054
4055 // We start by checking whether only the numbering within the MPI
4056 // ranks changed, in which case we do not need to find a new index
4057 // set.
4058 const IndexSet &owned_dofs = dof_handler->locally_owned_dofs();
4059 const bool locally_owned_set_changes =
4060 std::any_of(new_numbers.cbegin(),
4061 new_numbers.cend(),
4062 [&owned_dofs](const types::global_dof_index i) {
4063 return owned_dofs.is_element(i) == false;
4064 });
4065
4066 IndexSet my_locally_owned_new_dof_indices = owned_dofs;
4067 if (locally_owned_set_changes && owned_dofs.n_elements() > 0)
4068 {
4069 std::vector<::types::global_dof_index> new_numbers_sorted =
4070 new_numbers;
4071 std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4072
4073 my_locally_owned_new_dof_indices = IndexSet(dof_handler->n_dofs());
4074 my_locally_owned_new_dof_indices.add_indices(
4075 new_numbers_sorted.begin(), new_numbers_sorted.end());
4076 my_locally_owned_new_dof_indices.compress();
4077
4078 Assert(my_locally_owned_new_dof_indices.n_elements() ==
4079 new_numbers.size(),
4080 ExcInternalError());
4081 }
4082
4083 // delete all knowledge of DoF indices that are not locally
4084 // owned. we do so by getting DoF indices on cells, checking
4085 // whether they are locally owned, if not, setting them to
4086 // an invalid value, and then setting them again on the current
4087 // cell
4088 //
4089 // DoFs we (i) know about, and (ii) don't own locally must be
4090 // located either on ghost cells, or on the interface between a
4091 // locally owned cell and a ghost cell. In any case, it is
4092 // sufficient to kill them only from the ghost side cell, so loop
4093 // only over ghost cells
4094 for (auto cell : dof_handler->active_cell_iterators())
4095 if (cell->is_ghost())
4096 {
4097 DoFAccessorImplementation::Implementation::process_dof_indices(
4098 *cell,
4099 std::make_tuple(),
4100 cell->active_fe_index(),
4101 DoFAccessorImplementation::Implementation::
4102 DoFIndexProcessor<dim, spacedim>(),
4103 [&owned_dofs](auto &stored_index, auto) {
4104 // delete a DoF index if it has not already been
4105 // deleted (e.g., by visiting a neighboring cell, if
4106 // it is on the boundary), and if we don't own it
4107 if (stored_index != numbers::invalid_dof_index &&
4108 (!owned_dofs.is_element(stored_index)))
4109 stored_index = numbers::invalid_dof_index;
4110 },
4111 false);
4112 }
4113
4114
4115 // renumber. Skip when there is nothing to do because we own no DoF.
4116 if (owned_dofs.n_elements() > 0)
4117 Implementation::renumber_dofs(new_numbers,
4118 owned_dofs,
4119 *dof_handler,
4120 /*check_validity=*/false);
4121
4122 // Communicate newly assigned DoF indices to other processors
4123 // and get the same information for our own ghost cells.
4124 //
4125 // This is the same as phase 5+6 in the distribute_dofs() algorithm,
4126 // taking into account that we have to unify a few DoFs in between
4127 // then communication phases if we do hp-numbering
4128 {
4129 // mark all ghost cells for transfer
4130 std::vector<bool> cell_marked(triangulation->n_active_cells());
4131 for (const auto &cell : dof_handler->active_cell_iterators())
4132 if (cell->is_ghost())
4133 cell_marked[cell->active_cell_index()] = true;
4134
4135 // Send and receive cells. After this, only the local cells
4136 // are marked, that received new data. This has to be
4137 // communicated in a second communication step.
4138 //
4139 // as explained in the 'distributed' paper, this has to be
4140 // done twice
4141 communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
4142
4143 // if the DoFHandler has hp-capabilities then we may have
4144 // received valid indices of degrees of freedom that are
4145 // dominated by a FE object adjacent to a ghost interface.
4146 // thus, we overwrite the remaining invalid indices with the
4147 // valid ones in this step.
4148 Implementation::merge_invalid_dof_indices_on_ghost_interfaces(
4149 *dof_handler);
4150
4151 communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
4152 }
4153
4154 NumberCache number_cache;
4155 number_cache.locally_owned_dofs =
4156 std::move(my_locally_owned_new_dof_indices);
4157 number_cache.n_global_dofs = dof_handler->n_dofs();
4158 number_cache.n_locally_owned_dofs =
4159 number_cache.locally_owned_dofs.n_elements();
4160 return number_cache;
4161#endif
4162 }
4163
4164
4165
4166 template <int dim, int spacedim>
4167 NumberCache
4168 ParallelDistributed<dim, spacedim>::renumber_mg_dofs(
4169 const unsigned int level,
4170 const std::vector<types::global_dof_index> &new_numbers) const
4171 {
4172#ifndef DEAL_II_WITH_MPI
4173
4174 (void)level;
4175 (void)new_numbers;
4176
4177 DEAL_II_NOT_IMPLEMENTED();
4178 return NumberCache();
4179#else
4180
4181 ::parallel::DistributedTriangulationBase<dim, spacedim>
4182 *triangulation =
4183 (dynamic_cast<
4184 ::parallel::DistributedTriangulationBase<dim, spacedim> *>(
4185 const_cast<::Triangulation<dim, spacedim> *>(
4186 &dof_handler->get_triangulation())));
4187 Assert(triangulation != nullptr, ExcInternalError());
4188
4189 // This code is very close to the respective code in renumber_dofs,
4190 // with the difference that we work on different entities with
4191 // different objects.
4192 const IndexSet &owned_dofs = dof_handler->locally_owned_mg_dofs(level);
4193 AssertDimension(new_numbers.size(), owned_dofs.n_elements());
4194
4195 const bool locally_owned_set_changes =
4196 std::any_of(new_numbers.cbegin(),
4197 new_numbers.cend(),
4198 [&owned_dofs](const types::global_dof_index i) {
4199 return owned_dofs.is_element(i) == false;
4200 });
4201
4202 IndexSet my_locally_owned_new_dof_indices = owned_dofs;
4203 if (locally_owned_set_changes && owned_dofs.n_elements() > 0)
4204 {
4205 std::vector<::types::global_dof_index> new_numbers_sorted =
4206 new_numbers;
4207 std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4208
4209 my_locally_owned_new_dof_indices =
4210 IndexSet(dof_handler->n_dofs(level));
4211 my_locally_owned_new_dof_indices.add_indices(
4212 new_numbers_sorted.begin(), new_numbers_sorted.end());
4213 my_locally_owned_new_dof_indices.compress();
4214
4215 Assert(my_locally_owned_new_dof_indices.n_elements() ==
4216 new_numbers.size(),
4217 ExcInternalError());
4218 }
4219
4220 // delete all knowledge of DoF indices that are not locally
4221 // owned
4222 for (auto cell : dof_handler->cell_iterators_on_level(level))
4223 if (cell->is_ghost_on_level())
4224 {
4225 DoFAccessorImplementation::Implementation::process_dof_indices(
4226 *cell,
4227 std::make_tuple(),
4228 0,
4229 DoFAccessorImplementation::Implementation::
4230 MGDoFIndexProcessor<dim, spacedim>(cell->level()),
4231 [&owned_dofs](auto &stored_index, auto) {
4232 if ((stored_index != numbers::invalid_dof_index) &&
4233 (!owned_dofs.is_element(stored_index)))
4234 stored_index = numbers::invalid_dof_index;
4235 },
4236 true);
4237 }
4238
4239 // renumber. Skip when there is nothing to do because we own no DoF.
4240 if (level < triangulation->n_levels() && owned_dofs.n_elements() > 0)
4241 Implementation::renumber_mg_dofs(
4242 new_numbers, owned_dofs, *dof_handler, level, false);
4243
4244 // communicate newly assigned DoF indices with other processors
4245 {
4246 std::vector<std::vector<bool>> cell_marked(triangulation->n_levels());
4247 for (unsigned int l = 0; l < triangulation->n_levels(); ++l)
4248 cell_marked[l].resize(triangulation->n_raw_cells(l));
4249 for (const auto &cell : dof_handler->cell_iterators_on_level(level))
4250 if (cell->is_ghost_on_level())
4251 cell_marked[cell->level()][cell->index()] = true;
4252
4253 communicate_mg_ghost_cells(*dof_handler, cell_marked);
4254
4255 communicate_mg_ghost_cells(*dof_handler, cell_marked);
4256 }
4257
4258 NumberCache number_cache;
4259 number_cache.locally_owned_dofs =
4260 std::move(my_locally_owned_new_dof_indices);
4261 number_cache.n_global_dofs = dof_handler->n_dofs(level);
4262 number_cache.n_locally_owned_dofs =
4263 number_cache.locally_owned_dofs.n_elements();
4264 return number_cache;
4265#endif
4266 }
4267 } // namespace Policy
4268 } // namespace DoFHandlerImplementation
4269} // namespace internal
4270
4271
4272
4273/*-------------- Explicit Instantiations -------------------------------*/
4274#include "dofs/dof_handler_policy.inst"
4275
4276
4277DEAL_II_NAMESPACE_CLOSE
DoFHandler::mg_levels
std::vector< std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel< dim > > > mg_levels
Definition dof_handler.h:1563
DoFHandler::fe_collection
hp::FECollection< dim, spacedim > fe_collection
Definition dof_handler.h:1483
DoFHandler::get_fe_collection
const hp::FECollection< dim, spacedim > & get_fe_collection() const
DoFHandler::get_fe
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
DoFHandler::locally_owned_mg_dofs
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
DoFHandler::mg_vertex_dofs
std::vector< MGVertexDoFs > mg_vertex_dofs
Definition dof_handler.h:1556
DoFHandler::object_dof_indices
std::vector< std::array< std::vector< types::global_dof_index >, dim+1 > > object_dof_indices
Definition dof_handler.h:1514
DoFHandler::get_triangulation
const Triangulation< dim, spacedim > & get_triangulation() const
DoFHandler::locally_owned_dofs
const IndexSet & locally_owned_dofs() const
DoFHandler::mg_faces
std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces< dim > > mg_faces
Definition dof_handler.h:1569
DoFHandler::hp_capability_enabled
bool hp_capability_enabled
Definition dof_handler.h:1470
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
DoFHandler::n_locally_owned_dofs
types::global_dof_index n_locally_owned_dofs() const
FiniteElementData::n_dofs_per_vertex
unsigned int n_dofs_per_vertex() const
FiniteElementData::n_dofs_per_line
unsigned int n_dofs_per_line() const
FiniteElementData::max_dofs_per_quad
unsigned int max_dofs_per_quad() const
FiniteElementData::n_dofs_per_quad
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
IndexSet::size
size_type size() const
Definition index_set.h:1759
IndexSet::index_within_set
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1977
IndexSet::n_elements
size_type n_elements() const
Definition index_set.h:1917
IndexSet::is_element
bool is_element(const size_type index) const
Definition index_set.h:1877
IndexSet::add_range
void add_range(const size_type begin, const size_type end)
Definition index_set.h:1786
IndexSet::nth_index_in_set
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1958
IndexSet::compress
void compress() const
Definition index_set.h:1767
IndexSet::add_indices
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1814
Triangulation::begin
cell_iterator begin(const unsigned int level=0) const
Triangulation::n_raw_lines
unsigned int n_raw_lines() const
Triangulation::locally_owned_subdomain
virtual types::subdomain_id locally_owned_subdomain() const
Triangulation::n_active_cells
unsigned int n_active_cells() const
Triangulation::n_levels
unsigned int n_levels() const
Triangulation::end
cell_iterator end() const
Triangulation::n_raw_cells
unsigned int n_raw_cells(const unsigned int level) const
Triangulation::vertex_used
bool vertex_used(const unsigned int index) const
Triangulation::n_global_levels
virtual unsigned int n_global_levels() const
Triangulation::n_raw_quads
unsigned int n_raw_quads() const
Triangulation::get_used_vertices
const std::vector< bool > & get_used_vertices() const
Triangulation::n_vertices
unsigned int n_vertices() const
hp::Collection::size
unsigned int size() const
Definition collection.h:314
hp::FECollection::find_dominating_fe
unsigned int find_dominating_fe(const std::set< unsigned int > &fes, const unsigned int codim=0) const
Definition fe_collection.cc:197
hp::FECollection::max_dofs_per_cell
unsigned int max_dofs_per_cell() const
Definition fe_collection.h:1040
internal::DoFHandlerImplementation::Policy::ParallelDistributed::renumber_dofs
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
Definition dof_handler_policy.cc:4033
internal::DoFHandlerImplementation::Policy::ParallelDistributed::distribute_mg_dofs
virtual std::vector< NumberCache > distribute_mg_dofs() const override
Definition dof_handler_policy.cc:3854
internal::DoFHandlerImplementation::Policy::ParallelDistributed::renumber_mg_dofs
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
Definition dof_handler_policy.cc:4168
internal::DoFHandlerImplementation::Policy::ParallelShared::distribute_mg_dofs
virtual std::vector< NumberCache > distribute_mg_dofs() const override
Definition dof_handler_policy.cc:3062
internal::DoFHandlerImplementation::Policy::ParallelShared::renumber_mg_dofs
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
Definition dof_handler_policy.cc:3423
internal::DoFHandlerImplementation::Policy::ParallelShared::ParallelShared
ParallelShared(DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:2772
internal::DoFHandlerImplementation::Policy::ParallelShared::renumber_dofs
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
Definition dof_handler_policy.cc:3256
internal::DoFHandlerImplementation::Policy::Sequential::distribute_mg_dofs
virtual std::vector< NumberCache > distribute_mg_dofs() const override
Definition dof_handler_policy.cc:2705
internal::DoFHandlerImplementation::Policy::Sequential::Sequential
Sequential(DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:2677
internal::DoFHandlerImplementation::Policy::Sequential::renumber_dofs
virtual NumberCache renumber_dofs(const std::vector< types::global_dof_index > &new_numbers) const override
Definition dof_handler_policy.cc:2729
internal::DoFHandlerImplementation::Policy::Sequential::renumber_mg_dofs
virtual NumberCache renumber_mg_dofs(const unsigned int level, const std::vector< types::global_dof_index > &new_numbers) const override
Definition dof_handler_policy.cc:2756
parallel::TriangulationBase::locally_owned_subdomain
types::subdomain_id locally_owned_subdomain() const override
Definition tria_base.cc:315
parallel::TriangulationBase::n_global_levels
virtual unsigned int n_global_levels() const override
Definition tria_base.cc:139
parallel::TriangulationBase::get_mpi_communicator
virtual MPI_Comm get_mpi_communicator() const override
Definition tria_base.cc:158
parallel::TriangulationBase::is_multilevel_hierarchy_constructed
virtual bool is_multilevel_hierarchy_constructed() const =0
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
running_in_debug_mode
constexpr bool running_in_debug_mode()
Definition config.h:76
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
dof_accessor.h
dof_handler.h
dof_handler_policy.h
DEAL_II_NOT_IMPLEMENTED
#define DEAL_II_NOT_IMPLEMENTED()
Definition exception_macros.h:705
geometry_info.h
second
Point< 2 > second
Definition grid_out.cc:4640
first
Point< 2 > first
Definition grid_out.cc:4639
level
unsigned int level
Definition grid_out.cc:4642
grid_tools.h
DoFHandler::active_cell_iterators
IteratorRange< active_cell_iterator > active_cell_iterators() const
DoFHandler::cell_iterators_on_level
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
DoFHandler::cell_iterators
IteratorRange< cell_iterator > cell_iterators() const
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertThrowMPI
#define AssertThrowMPI(error_code)
Definition exception_macros.h:963
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
Threads::new_task
Task< RT > new_task(const std::function< RT()> &function)
Definition thread_management.h:1173
memory_consumption.h
data
std::vector< index_type > data
Definition mpi.cc:734
n_procs
const unsigned int n_procs
Definition mpi.cc:923
DoFTools::locally_owned_dofs_per_subdomain
std::vector< IndexSet > locally_owned_dofs_per_subdomain(const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_tools.cc:1784
GridTools::exchange_cell_data_to_level_ghosts
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
GridTools::exchange_cell_data_to_ghosts
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})
PlasticityLab::const
*  *  *  *  std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters   const
Definition ALE_Finite_Strain_Plasticity.h:9082
ReferenceCells::Quadrilateral
constexpr ReferenceCell< 2 > Quadrilateral
Definition reference_cell.h:1350
Utilities::MPI::partial_and_total_sum
std::pair< T, T > partial_and_total_sum(const T &value, const MPI_Comm comm)
Utilities::MPI::logical_and
bool logical_and(const bool t, const MPI_Comm mpi_communicator)
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:103
Utilities::MPI::all_gather
std::vector< T > all_gather(const MPI_Comm comm, const T &object_to_send)
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:118
Utilities::pack
std::size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition utilities.h:1352
Utilities::unpack
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
Definition utilities.h:1509
numbers::invalid_dof_index
constexpr types::global_dof_index invalid_dof_index
Definition types.h:259
numbers::invalid_unsigned_int
constexpr unsigned int invalid_unsigned_int
Definition types.h:228
numbers::artificial_subdomain_id
constexpr types::subdomain_id artificial_subdomain_id
Definition types.h:406
numbers::invalid_subdomain_id
constexpr types::subdomain_id invalid_subdomain_id
Definition types.h:385
numbers::invalid_fe_index
constexpr types::fe_index invalid_fe_index
Definition types.h:250
types::fe_index
unsigned short int fe_index
Definition types.h:70
shared_tria.h
internal::DoFAccessorImplementation::Implementation::process_dof_indices
static void process_dof_indices(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > &accessor, const DoFIndicesType &const_dof_indices, const types::fe_index fe_index_, const DoFOperation &dof_operation, const DoFProcessor &dof_processor, const bool count_level_dofs)
Definition dof_accessor.h:3132
internal::DoFAccessorImplementation::Implementation::get_dof_indices
static void get_dof_indices(const ::DoFAccessor< structdim, dim, spacedim, level_dof_access > &accessor, std::vector< types::global_dof_index > &dof_indices, const types::fe_index fe_index)
Definition dof_accessor.h:3417
internal::DoFAccessorImplementation::Implementation::get_active_fe_indices
static std::set< types::fe_index > get_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const std::integral_constant< int, structdim > &t)
Definition dof_accessor.h:2842
internal::DoFAccessorImplementation::Implementation::n_active_fe_indices
static unsigned int n_active_fe_indices(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const std::integral_constant< int, structdim > &)
Definition dof_accessor.h:2751
internal::DoFAccessorImplementation::Implementation::set_dof_index
static void set_dof_index(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, structdim > &dd, const types::global_dof_index global_index)
Definition dof_accessor.h:2672
internal::DoFAccessorImplementation::Implementation::get_mg_dof_index
static types::global_dof_index & get_mg_dof_index(const DoFHandler< dim, spacedim > &dof_handler, const std::unique_ptr< internal::DoFHandlerImplementation::DoFLevel< dim > > &mg_level, const std::unique_ptr< internal::DoFHandlerImplementation::DoFFaces< dim > > &, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, dim >)
Definition dof_accessor.h:3520
internal::DoFAccessorImplementation::Implementation::nth_active_fe_index
static types::fe_index nth_active_fe_index(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const unsigned int local_index, const std::integral_constant< int, structdim > &)
Definition dof_accessor.h:2790
internal::DoFAccessorImplementation::Implementation::get_dof_index
static types::global_dof_index get_dof_index(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int obj_level, const unsigned int obj_index, const types::fe_index fe_index, const unsigned int local_index, const std::integral_constant< int, structdim > &dd)
Definition dof_accessor.h:2703
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 3, spacedim > &dof_handler)
Definition dof_handler_policy.cc:2191
internal::DoFHandlerImplementation::Policy::Implementation::compute_quad_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:737
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_quad_dofs_on_ghost_interfaces
static void merge_invalid_quad_dofs_on_ghost_interfaces(DoFHandler< 3, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1495
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_mg_dofs
static void renumber_face_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
Definition dof_handler_policy.cc:2518
internal::DoFHandlerImplementation::Policy::Implementation::renumber_dofs
static void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, const DoFHandler< dim, space_dim > &dof_handler, const bool check_validity)
Definition dof_handler_policy.cc:2369
internal::DoFHandlerImplementation::Policy::Implementation::compute_line_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const DoFHandler< 1, spacedim > &dof_handler)
Definition dof_handler_policy.cc:381
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_line_dofs_on_ghost_interfaces
static void merge_invalid_line_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1224
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_dof_indices_on_ghost_interfaces
static void merge_invalid_dof_indices_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1635
internal::DoFHandlerImplementation::Policy::Implementation::renumber_vertex_dofs
static void renumber_vertex_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const bool check_validity)
Definition dof_handler_policy.cc:1853
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_vertex_dofs_on_ghost_interfaces
static void merge_invalid_vertex_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1052
internal::DoFHandlerImplementation::Policy::Implementation::renumber_cell_mg_dofs
static void renumber_cell_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
Definition dof_handler_policy.cc:2468
internal::DoFHandlerImplementation::Policy::Implementation::compute_quad_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_quad_dof_identities(const DoFHandler< 3, spacedim > &dof_handler)
Definition dof_handler_policy.cc:756
internal::DoFHandlerImplementation::Policy::Implementation::renumber_mg_dofs
static void renumber_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level, const bool check_validity)
Definition dof_handler_policy.cc:2629
internal::DoFHandlerImplementation::Policy::Implementation::renumber_vertex_mg_dofs
static void renumber_vertex_mg_dofs(const std::vector<::types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
Definition dof_handler_policy.cc:2417
internal::DoFHandlerImplementation::Policy::Implementation::compute_vertex_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_vertex_dof_identities(const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:219
internal::DoFHandlerImplementation::Policy::Implementation::enumerate_dof_indices_for_renumbering
static types::global_dof_index enumerate_dof_indices_for_renumbering(std::vector< types::global_dof_index > &new_dof_indices, const std::vector< std::map< types::global_dof_index, types::global_dof_index > > &all_constrained_indices, const types::global_dof_index start_dof_index)
Definition dof_handler_policy.cc:957
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_line_dofs_on_ghost_interfaces
static void merge_invalid_line_dofs_on_ghost_interfaces(DoFHandler< 1, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1213
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< 2, spacedim > &dof_handler)
Definition dof_handler_policy.cc:2096
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
static void renumber_face_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1836
internal::DoFHandlerImplementation::Policy::Implementation::compute_dof_identities
static void compute_dof_identities(std::vector< std::map< types::global_dof_index, types::global_dof_index > > &all_constrained_indices, const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:896
internal::DoFHandlerImplementation::Policy::Implementation::unify_dof_indices
static types::global_dof_index unify_dof_indices(const DoFHandler< dim, spacedim > &dof_handler, const types::global_dof_index n_dofs_before_identification, const bool check_validity)
Definition dof_handler_policy.cc:1019
internal::DoFHandlerImplementation::Policy::Implementation::renumber_cell_dofs
static void renumber_cell_dofs(const std::vector< types::global_dof_index > &new_numbers, const IndexSet &indices_we_care_about, DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:2010
internal::DoFHandlerImplementation::Policy::Implementation::compute_line_dof_identities
static std::map< types::global_dof_index, types::global_dof_index > compute_line_dof_identities(const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:393
internal::DoFHandlerImplementation::Policy::Implementation::invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering
static void invalidate_dof_indices_on_weaker_ghost_cells_for_renumbering(std::vector< types::global_dof_index > &renumbering, const types::subdomain_id subdomain_id, const DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1739
internal::DoFHandlerImplementation::Policy::Implementation::distribute_dofs
static types::global_dof_index distribute_dofs(const types::subdomain_id subdomain_id, DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1676
internal::DoFHandlerImplementation::Policy::Implementation::merge_invalid_quad_dofs_on_ghost_interfaces
static void merge_invalid_quad_dofs_on_ghost_interfaces(DoFHandler< dim, spacedim > &dof_handler)
Definition dof_handler_policy.cc:1478
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_mg_dofs
static void renumber_face_mg_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &, const unsigned int, const bool)
Definition dof_handler_policy.cc:2504
internal::DoFHandlerImplementation::Policy::Implementation::distribute_dofs_on_level
static types::global_dof_index distribute_dofs_on_level(const types::subdomain_id level_subdomain_id, DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
Definition dof_handler_policy.cc:1774
internal::DoFHandlerImplementation::Policy::Implementation::renumber_face_dofs
static void renumber_face_dofs(const std::vector< types::global_dof_index > &, const IndexSet &, DoFHandler< 1, spacedim > &)
Definition dof_handler_policy.cc:2083
thread_management.h
tria_iterator.h
work_stream.h
Generated by doxygen 1.9.8