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