Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3087-ga35b476a78 2025-04-19 20:30:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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(
3315 new_numbers_copy.data(),
3316 new_numbers_copy.size(),
3317 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
3318 gathered_new_numbers.data(),
3319 rcounts.data(),
3320 displacements.data(),
3321 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
3322 tr->get_mpi_communicator());
3323 AssertThrowMPI(ierr);
3324 }
3325
3326 // put new numbers according to the current
3327 // locally_owned_dofs_per_processor IndexSets
3328 types::global_dof_index shift = 0;
3329 // flag_1 and flag_2 are
3330 // used to control that there is a
3331 // one-to-one relation between old and new DoFs.
3332 std::vector<unsigned int> flag_1(this->dof_handler->n_dofs(), 0);
3333 std::vector<unsigned int> flag_2(this->dof_handler->n_dofs(), 0);
3334 std::vector<IndexSet> locally_owned_dofs_per_processor =
3336 tr->get_mpi_communicator(),
3337 this->dof_handler->locally_owned_dofs());
3338 for (unsigned int i = 0; i < n_cpu; ++i)
3339 {
3340 const IndexSet &iset = locally_owned_dofs_per_processor[i];
3341 for (types::global_dof_index ind = 0; ind < iset.n_elements();
3342 ind++)
3343 {
3344 const types::global_dof_index target =
3345 iset.nth_index_in_set(ind);
3347 gathered_new_numbers[shift + ind];
3348 Assert(target < this->dof_handler->n_dofs(),
3350 Assert(value < this->dof_handler->n_dofs(),
3352 global_gathered_numbers[target] = value;
3353 flag_1[target]++;
3354 flag_2[value]++;
3355 }
3356 shift += iset.n_elements();
3357 }
3358
3359 Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1,
3361 Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1,
3363 Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1,
3365 Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1,
3367 }
3368
3369 // let the sequential algorithm do its magic; ignore the
3370 // return type, but reconstruct the number cache based on
3371 // which DoFs each process owns
3372 Implementation::renumber_dofs(global_gathered_numbers,
3373 IndexSet(0),
3374 *this->dof_handler,
3375 /*check_validity=*/true);
3376
3377 const NumberCache number_cache(
3379 this->dof_handler->get_triangulation().locally_owned_subdomain());
3380
3381 return number_cache;
3382#endif
3383 }
3384
3385
3386
3387 template <int dim, int spacedim>
3390 const unsigned int /*level*/,
3391 const std::vector<types::global_dof_index> & /*new_numbers*/) const
3392 {
3393 // multigrid is not currently implemented for shared triangulations
3395
3396 return {};
3397 }
3398
3399
3400
3401 /* --------------------- class ParallelDistributed ---------------- */
3402
3403#ifdef DEAL_II_WITH_MPI
3404
3405 namespace
3406 {
3407 template <int dim, int spacedim>
3408 void
3409 communicate_mg_ghost_cells(DoFHandler<dim, spacedim> &dof_handler,
3410 std::vector<std::vector<bool>> &cell_marked)
3411 {
3412 const auto pack = [](const auto &cell) {
3413 // why would somebody request a cell that is not ours?
3414 Assert(cell->is_locally_owned_on_level(), ExcInternalError());
3415
3416 std::vector<::types::global_dof_index> data(
3417 cell->get_fe().n_dofs_per_cell());
3418 cell->get_mg_dof_indices(data);
3419
3420 return data;
3421 };
3422
3423 const auto unpack = [&cell_marked](const auto &cell,
3424 const auto &dofs) {
3425 Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3427
3428 Assert(cell->level_subdomain_id() !=
3431
3432 bool complete = true;
3434 *cell,
3435 dofs,
3436 0,
3438 MGDoFIndexProcessor<dim, spacedim>(cell->level()),
3439
3440 // Intel ICC 18 and earlier for some reason believe that
3441 // numbers::invalid_dof_index is not a valid object
3442 // inside the lambda function. Fix this by creating a
3443 // local variable initialized by the global one.
3444 //
3445 // Intel ICC 19 and earlier have trouble with our Assert
3446 // macros inside the lambda function. We disable the macro
3447 // for these compilers.
3448 [&complete, invalid_dof_index = numbers::invalid_dof_index](
3449 auto &stored_index, auto received_index) {
3450 if (*received_index != invalid_dof_index)
3451 {
3452# if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
3453 Assert((stored_index == invalid_dof_index) ||
3454 (stored_index == *received_index),
3456# endif
3457 stored_index = *received_index;
3458 }
3459 else
3460 complete = false;
3461 },
3462 true);
3463
3464 if (!complete)
3465 {
3466 // We should have the cell already marked
3467 Assert(cell_marked[cell->level()][cell->index()],
3469 }
3470 else
3471 cell_marked[cell->level()][cell->index()] = false;
3472 };
3473
3474 const auto filter = [&cell_marked](const auto &cell) {
3475 return cell_marked[cell->level()][cell->index()];
3476 };
3477
3479 std::vector<types::global_dof_index>,
3480 DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter);
3481 }
3482
3483
3484
3503 template <int dim, int spacedim>
3504 void
3505 communicate_dof_indices_on_marked_cells(
3506 const DoFHandler<dim, spacedim> &dof_handler,
3507 std::vector<bool> &cell_marked)
3508 {
3509# ifndef DEAL_II_WITH_MPI
3510 (void)dof_handler;
3512# else
3513
3514 // define functions that pack data on cells that are ghost cells
3515 // somewhere else, and unpack data on cells where we get information
3516 // from elsewhere
3517 const auto pack = [](const auto &cell) {
3518 Assert(cell->is_locally_owned(), ExcInternalError());
3519
3520 std::vector<::types::global_dof_index> data(
3521 cell->get_fe().n_dofs_per_cell());
3522
3523 // bypass the cache which is not filled yet
3525 get_dof_indices(*cell, data, cell->active_fe_index());
3526
3527 return data;
3528 };
3529
3530 const auto unpack = [&cell_marked](const auto &cell,
3531 const auto &dofs) {
3532 Assert(cell->get_fe().n_dofs_per_cell() == dofs.size(),
3534
3535 Assert(cell->is_ghost(), ExcInternalError());
3536
3537 // Use a combined read/set function on the entities of the dof
3538 // indices to speed things up against get_dof_indices +
3539 // set_dof_indices
3540 bool complete = true;
3542 *cell,
3543 dofs,
3544 cell->active_fe_index(),
3545 DoFAccessorImplementation::Implementation::
3546 DoFIndexProcessor<dim, spacedim>(),
3547
3548 // Intel ICC 18 and earlier for some reason believe that
3549 // numbers::invalid_dof_index is not a valid object
3550 // inside the lambda function. Fix this by creating a
3551 // local variable initialized by the global one.
3552 //
3553 // Intel ICC 19 and earlier have trouble with our Assert
3554 // macros inside the lambda function. We disable the macro
3555 // for these compilers.
3557 auto &stored_index, const auto received_index) {
3558 if (*received_index != invalid_dof_index)
3559 {
3560# if !defined(__INTEL_COMPILER) || __INTEL_COMPILER >= 1900
3561 Assert((stored_index == invalid_dof_index) ||
3562 (stored_index == *received_index),
3564# endif
3565 stored_index = *received_index;
3566 }
3567 else
3568 complete = false;
3569 },
3570 false);
3571
3572 if (!complete)
3573 {
3574 // We should have the cell already marked
3575 Assert(cell_marked[cell->active_cell_index()],
3577 }
3578 else
3579 cell_marked[cell->active_cell_index()] = false;
3580 };
3581
3582 const auto filter = [&cell_marked](const auto &cell) {
3583 return cell_marked[cell->active_cell_index()];
3584 };
3585
3587 std::vector<types::global_dof_index>,
3588 DoFHandler<dim, spacedim>>(dof_handler, pack, unpack, filter);
3589# endif
3590 }
3591
3592
3593
3594 } // namespace
3595
3596#endif // DEAL_II_WITH_MPI
3597
3598
3599
3600 template <int dim, int spacedim>
3602 DoFHandler<dim, spacedim> &dof_handler)
3603 : dof_handler(&dof_handler)
3604 {}
3605
3606
3607
3608 template <int dim, int spacedim>
3611 {
3612#ifndef DEAL_II_WITH_MPI
3614 return NumberCache();
3615#else
3616
3618 *triangulation =
3619 (dynamic_cast<
3621 const_cast<::Triangulation<dim, spacedim> *>(
3622 &dof_handler->get_triangulation())));
3623 Assert(triangulation != nullptr, ExcInternalError());
3624
3625 const types::subdomain_id subdomain_id =
3626 triangulation->locally_owned_subdomain();
3627
3628
3629 /*
3630 The following algorithm has a number of stages that are all
3631 documented in the paper that describes the parallel::distributed
3632 functionality:
3633
3634 1/ locally enumerate dofs on locally owned cells
3635 2/ eliminate dof duplicates on all cells.
3636 un-numerate those that are on interfaces with ghost
3637 cells and that we don't own based on the tie-breaking
3638 criterion. unify dofs afterwards.
3639 3/ unify dofs and re-enumerate the remaining valid ones.
3640 the end result is that we only enumerate locally owned
3641 DoFs
3642 4/ shift indices so that each processor has a unique
3643 range of indices
3644 5/ for all locally owned cells that are ghost
3645 cells somewhere else, send our own DoF indices
3646 to the appropriate set of other processors.
3647 overwrite invalid DoF indices on ghost interfaces
3648 with the corresponding valid ones that we now know.
3649 6/ send DoF indices again to get the correct indices
3650 on ghost cells that we may not have known earlier
3651 */
3652
3653 // --------- Phase 1: enumerate dofs on locally owned cells
3654 const types::global_dof_index n_initial_local_dofs =
3655 Implementation::distribute_dofs(subdomain_id, *dof_handler);
3656
3657 // --------- Phase 2: eliminate dof duplicates on all cells:
3658 // - un-numerate dofs on interfaces to ghost cells
3659 // that we don't own
3660 // - in case of hp-support, unify dofs
3661 std::vector<::types::global_dof_index> renumbering(
3662 n_initial_local_dofs, enumeration_dof_index);
3663
3664 // first, we invalidate degrees of freedom that belong to processors
3665 // of a lower rank, from which we will receive the final (and lower)
3666 // degrees of freedom later.
3669 renumbering, subdomain_id, *dof_handler);
3670
3671 // then, we identify DoF duplicates if the DoFHandler has hp-
3672 // capabilities
3673 std::vector<std::map<types::global_dof_index, types::global_dof_index>>
3674 all_constrained_indices(dim);
3675 Implementation::compute_dof_identities(all_constrained_indices,
3676 *dof_handler);
3677
3678 // --------- Phase 3: re-enumerate the valid degrees of freedom
3679 // consecutively. thus, we finally receive the
3680 // correct number of locally owned DoFs after
3681 // this step.
3682 //
3683 // the order in which we handle Phases 2 and 3 is important,
3684 // since we want to clarify ownership of degrees of freedom before
3685 // we actually unify and enumerate their indices. otherwise, we could
3686 // end up having a degree of freedom to which only invalid indices will
3687 // be assigned.
3688 types::global_dof_index n_identity_constrained_indices = 0;
3689 for (const auto &constrained_indices : all_constrained_indices)
3690 for (const auto &index : constrained_indices)
3691 if (renumbering[index.first] != numbers::invalid_dof_index)
3692 ++n_identity_constrained_indices;
3693
3694 const types::global_dof_index n_locally_owned_dofs =
3695 std::count(renumbering.begin(),
3696 renumbering.end(),
3697 enumeration_dof_index) -
3698 n_identity_constrained_indices;
3699
3700 // --------- Phase 4: shift indices so that each processor has a unique
3701 // range of indices
3702 const auto [my_shift, n_global_dofs] =
3704 n_locally_owned_dofs, triangulation->get_mpi_communicator());
3705
3706
3707 // make dof indices globally consecutive
3709 renumbering, all_constrained_indices, my_shift);
3710
3711 // now re-enumerate all dofs to this shifted and condensed
3712 // numbering form. we renumber some dofs as invalid, so
3713 // choose the nocheck-version.
3715 IndexSet(0),
3716 *dof_handler,
3717 /*check_validity=*/false);
3718
3719 NumberCache number_cache;
3720 number_cache.n_global_dofs = n_global_dofs;
3721 number_cache.n_locally_owned_dofs = n_locally_owned_dofs;
3722 number_cache.locally_owned_dofs = IndexSet(n_global_dofs);
3723 number_cache.locally_owned_dofs.add_range(my_shift,
3724 my_shift +
3725 n_locally_owned_dofs);
3726 number_cache.locally_owned_dofs.compress();
3727
3728 // this ends the phase where we enumerate degrees of freedom on
3729 // each processor. what is missing is communicating DoF indices
3730 // on ghost cells
3731
3732 // --------- Phase 5: for all locally owned cells that are ghost
3733 // cells somewhere else, send our own DoF indices
3734 // to the appropriate set of other processors
3735 {
3736 // mark all cells that either have to send data (locally
3737 // owned cells that are adjacent to ghost neighbors in some
3738 // way) or receive data (all ghost cells) via the user flags
3739 std::vector<bool> cell_marked(triangulation->n_active_cells());
3740 for (const auto &cell : dof_handler->active_cell_iterators())
3741 if (cell->is_ghost())
3742 cell_marked[cell->active_cell_index()] = true;
3743
3744 // Send and receive cells. After this, only the local cells
3745 // are marked, that received new data. This has to be
3746 // communicated in a second communication step.
3747 //
3748 // as explained in the 'distributed' paper, this has to be
3749 // done twice
3750 communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
3751
3752 // If the DoFHandler has hp-capabilities enabled, then we may have
3753 // received valid indices of degrees of freedom that are dominated
3754 // by a FE object adjacent to a ghost interface. thus, we overwrite
3755 // the remaining invalid indices with the valid ones in this step.
3757 *dof_handler);
3758
3759 // --------- Phase 6: all locally owned cells have their correct
3760 // DoF indices set. however, some ghost cells
3761 // may still have invalid ones. thus, exchange
3762 // one more time.
3763 communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
3764
3765 // at this point, we must have taken care of the data transfer
3766 // on all cells we had previously marked. verify this
3767 if constexpr (running_in_debug_mode())
3768 {
3769 for (const auto &cell : dof_handler->active_cell_iterators())
3770 Assert(cell_marked[cell->active_cell_index()] == false,
3772 }
3773 }
3774
3775 if constexpr (running_in_debug_mode())
3776 {
3777 // check that we are really done
3778 {
3779 std::vector<::types::global_dof_index> local_dof_indices;
3780
3781 for (const auto &cell : dof_handler->active_cell_iterators())
3782 if (!cell->is_artificial())
3783 {
3784 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
3785 cell->get_dof_indices(local_dof_indices);
3786 if (local_dof_indices.end() !=
3787 std::find(local_dof_indices.begin(),
3788 local_dof_indices.end(),
3790 {
3791 if (cell->is_ghost())
3792 {
3793 Assert(false,
3794 ExcMessage(
3795 "A ghost cell ended up with incomplete "
3796 "DoF index information. This should not "
3797 "have happened!"));
3798 }
3799 else
3800 {
3801 Assert(
3802 false,
3803 ExcMessage(
3804 "A locally owned cell ended up with incomplete "
3805 "DoF index information. This should not "
3806 "have happened!"));
3807 }
3808 }
3809 }
3810 }
3811 } // DEBUG
3812 return number_cache;
3813#endif // DEAL_II_WITH_MPI
3814 }
3815
3816
3817
3818 template <int dim, int spacedim>
3819 std::vector<NumberCache>
3821 {
3822#ifndef DEAL_II_WITH_MPI
3824 return std::vector<NumberCache>();
3825#else
3826
3828 *triangulation =
3829 (dynamic_cast<
3831 const_cast<::Triangulation<dim, spacedim> *>(
3832 &dof_handler->get_triangulation())));
3833 Assert(triangulation != nullptr, ExcInternalError());
3834
3835 AssertThrow((triangulation->is_multilevel_hierarchy_constructed()),
3836 ExcMessage(
3837 "Multigrid DoFs can only be distributed on a parallel "
3838 "Triangulation if the flag construct_multigrid_hierarchy "
3839 "is set in the constructor."));
3840
3841 // loop over all levels that exist globally (across all
3842 // processors), even if the current processor does not in fact
3843 // have any cells on that level or if the local part of the
3844 // Triangulation has fewer levels. we need to do this because
3845 // we need to communicate across all processors on all levels
3846 const unsigned int n_levels = triangulation->n_global_levels();
3847 std::vector<NumberCache> number_caches;
3848 number_caches.reserve(n_levels);
3849 for (unsigned int level = 0; level < n_levels; ++level)
3850 {
3851 NumberCache level_number_cache;
3852
3853 //* 1. distribute on own subdomain
3854 const unsigned int n_initial_local_dofs =
3856 triangulation->locally_owned_subdomain(), *dof_handler, level);
3857
3858 //* 2. iterate over ghostcells and kill dofs that are not
3859 // owned by us, which we mark by invalid_dof_index
3860 std::vector<::types::global_dof_index> renumbering(
3861 n_initial_local_dofs, enumeration_dof_index);
3862
3863 if (level < triangulation->n_levels())
3864 {
3865 std::vector<::types::global_dof_index> local_dof_indices;
3866
3867 for (const auto &cell :
3868 dof_handler->cell_iterators_on_level(level))
3869 if (cell->level_subdomain_id() !=
3871 (cell->level_subdomain_id() <
3872 triangulation->locally_owned_subdomain()))
3873 {
3874 // we found a neighboring ghost cell whose
3875 // subdomain is "stronger" than our own
3876 // subdomain
3877
3878 // delete all dofs that live there and that we
3879 // have previously assigned a number to
3880 // (i.e. the ones on the interface)
3881 local_dof_indices.resize(
3882 cell->get_fe().n_dofs_per_cell());
3883 cell->get_mg_dof_indices(local_dof_indices);
3884 for (unsigned int i = 0;
3885 i < cell->get_fe().n_dofs_per_cell();
3886 ++i)
3887 if (local_dof_indices[i] != numbers::invalid_dof_index)
3888 renumbering[local_dof_indices[i]] =
3890 }
3891 }
3892
3893 level_number_cache.n_locally_owned_dofs =
3894 std::count(renumbering.begin(),
3895 renumbering.end(),
3896 enumeration_dof_index);
3897
3898 //* 3. communicate local dofcount and shift ids to make
3899 // them unique
3900 const auto [my_shift, n_global_dofs] =
3902 level_number_cache.n_locally_owned_dofs,
3903 triangulation->get_mpi_communicator());
3904 level_number_cache.n_global_dofs = n_global_dofs;
3905
3906 // assign appropriate indices
3907 types::global_dof_index next_free_index = my_shift;
3908 for (types::global_dof_index &index : renumbering)
3909 if (index == enumeration_dof_index)
3910 index = next_free_index++;
3911
3912 // now re-enumerate all dofs to this shifted and condensed
3913 // numbering form. we renumber some dofs as invalid, so
3914 // choose the nocheck-version of the function
3915 //
3916 // of course there is nothing for us to renumber if the
3917 // level we are currently dealing with doesn't even exist
3918 // within the current triangulation, so skip renumbering
3919 // in that case
3920 if (level < triangulation->n_levels())
3922 renumbering, IndexSet(0), *dof_handler, level, false);
3923
3924 // now a little bit of housekeeping
3925 level_number_cache.locally_owned_dofs =
3926 IndexSet(level_number_cache.n_global_dofs);
3927 level_number_cache.locally_owned_dofs.add_range(
3928 next_free_index - level_number_cache.n_locally_owned_dofs,
3929 next_free_index);
3930 level_number_cache.locally_owned_dofs.compress();
3931
3932 number_caches.emplace_back(level_number_cache);
3933 }
3934
3935
3936 //* communicate ghost DoFs
3937 // We mark all ghost cells by setting the user_flag and then request
3938 // these cells from the corresponding owners. As this information
3939 // can be incomplete,
3940 {
3941 std::vector<std::vector<bool>> cell_marked(triangulation->n_levels());
3942 for (unsigned int l = 0; l < triangulation->n_levels(); ++l)
3943 cell_marked[l].resize(triangulation->n_raw_cells(l));
3944 for (const auto &cell : dof_handler->cell_iterators())
3945 if (cell->is_ghost_on_level())
3946 cell_marked[cell->level()][cell->index()] = true;
3947
3948 // Phase 1. Request all marked cells from corresponding owners. If we
3949 // managed to get every DoF, remove the user_flag, otherwise we
3950 // will request them again in the step below.
3951 communicate_mg_ghost_cells(*dof_handler, cell_marked);
3952
3953 // Phase 2, only request the cells that were not completed
3954 // in Phase 1.
3955 communicate_mg_ghost_cells(*dof_handler, cell_marked);
3956
3957 if constexpr (running_in_debug_mode())
3958 {
3959 // make sure we have finished all cells:
3960 for (const auto &cell : dof_handler->cell_iterators())
3961 Assert(cell_marked[cell->level()][cell->index()] == false,
3963 }
3964 }
3965
3966
3967
3968 if constexpr (running_in_debug_mode())
3969 {
3970 // check that we are really done
3971 {
3972 std::vector<::types::global_dof_index> local_dof_indices;
3973 for (const auto &cell : dof_handler->cell_iterators())
3974 if (cell->level_subdomain_id() !=
3976 {
3977 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
3978 cell->get_mg_dof_indices(local_dof_indices);
3979 if (local_dof_indices.end() !=
3980 std::find(local_dof_indices.begin(),
3981 local_dof_indices.end(),
3983 {
3984 Assert(false,
3985 ExcMessage("not all DoFs got distributed!"));
3986 }
3987 }
3988 }
3989 } // DEBUG
3990
3991 return number_caches;
3992
3993#endif // DEAL_II_WITH_MPI
3994 }
3995
3996
3997 template <int dim, int spacedim>
4000 const std::vector<::types::global_dof_index> &new_numbers) const
4001 {
4002 (void)new_numbers;
4003
4004 Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
4006
4007#ifndef DEAL_II_WITH_MPI
4009 return NumberCache();
4010#else
4011
4013 *triangulation =
4014 (dynamic_cast<
4016 const_cast<::Triangulation<dim, spacedim> *>(
4017 &dof_handler->get_triangulation())));
4018 Assert(triangulation != nullptr, ExcInternalError());
4019
4020
4021 // We start by checking whether only the numbering within the MPI
4022 // ranks changed, in which case we do not need to find a new index
4023 // set.
4024 const IndexSet &owned_dofs = dof_handler->locally_owned_dofs();
4025 const bool locally_owned_set_changes =
4026 std::any_of(new_numbers.cbegin(),
4027 new_numbers.cend(),
4028 [&owned_dofs](const types::global_dof_index i) {
4029 return owned_dofs.is_element(i) == false;
4030 });
4031
4032 IndexSet my_locally_owned_new_dof_indices = owned_dofs;
4033 if (locally_owned_set_changes && owned_dofs.n_elements() > 0)
4034 {
4035 std::vector<::types::global_dof_index> new_numbers_sorted =
4036 new_numbers;
4037 std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4038
4039 my_locally_owned_new_dof_indices = IndexSet(dof_handler->n_dofs());
4040 my_locally_owned_new_dof_indices.add_indices(
4041 new_numbers_sorted.begin(), new_numbers_sorted.end());
4042 my_locally_owned_new_dof_indices.compress();
4043
4044 Assert(my_locally_owned_new_dof_indices.n_elements() ==
4045 new_numbers.size(),
4047 }
4048
4049 // delete all knowledge of DoF indices that are not locally
4050 // owned. we do so by getting DoF indices on cells, checking
4051 // whether they are locally owned, if not, setting them to
4052 // an invalid value, and then setting them again on the current
4053 // cell
4054 //
4055 // DoFs we (i) know about, and (ii) don't own locally must be
4056 // located either on ghost cells, or on the interface between a
4057 // locally owned cell and a ghost cell. In any case, it is
4058 // sufficient to kill them only from the ghost side cell, so loop
4059 // only over ghost cells
4060 for (auto cell : dof_handler->active_cell_iterators())
4061 if (cell->is_ghost())
4062 {
4064 *cell,
4065 std::make_tuple(),
4066 cell->active_fe_index(),
4068 DoFIndexProcessor<dim, spacedim>(),
4069 [&owned_dofs](auto &stored_index, auto) {
4070 // delete a DoF index if it has not already been
4071 // deleted (e.g., by visiting a neighboring cell, if
4072 // it is on the boundary), and if we don't own it
4073 if (stored_index != numbers::invalid_dof_index &&
4074 (!owned_dofs.is_element(stored_index)))
4075 stored_index = numbers::invalid_dof_index;
4076 },
4077 false);
4078 }
4079
4080
4081 // renumber. Skip when there is nothing to do because we own no DoF.
4082 if (owned_dofs.n_elements() > 0)
4084 owned_dofs,
4085 *dof_handler,
4086 /*check_validity=*/false);
4087
4088 // Communicate newly assigned DoF indices to other processors
4089 // and get the same information for our own ghost cells.
4090 //
4091 // This is the same as phase 5+6 in the distribute_dofs() algorithm,
4092 // taking into account that we have to unify a few DoFs in between
4093 // then communication phases if we do hp-numbering
4094 {
4095 // mark all ghost cells for transfer
4096 std::vector<bool> cell_marked(triangulation->n_active_cells());
4097 for (const auto &cell : dof_handler->active_cell_iterators())
4098 if (cell->is_ghost())
4099 cell_marked[cell->active_cell_index()] = true;
4100
4101 // Send and receive cells. After this, only the local cells
4102 // are marked, that received new data. This has to be
4103 // communicated in a second communication step.
4104 //
4105 // as explained in the 'distributed' paper, this has to be
4106 // done twice
4107 communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
4108
4109 // if the DoFHandler has hp-capabilities then we may have
4110 // received valid indices of degrees of freedom that are
4111 // dominated by a FE object adjacent to a ghost interface.
4112 // thus, we overwrite the remaining invalid indices with the
4113 // valid ones in this step.
4115 *dof_handler);
4116
4117 communicate_dof_indices_on_marked_cells(*dof_handler, cell_marked);
4118 }
4119
4120 NumberCache number_cache;
4121 number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
4122 number_cache.n_global_dofs = dof_handler->n_dofs();
4123 number_cache.n_locally_owned_dofs =
4124 number_cache.locally_owned_dofs.n_elements();
4125 return number_cache;
4126#endif
4127 }
4128
4129
4130
4131 template <int dim, int spacedim>
4134 const unsigned int level,
4135 const std::vector<types::global_dof_index> &new_numbers) const
4136 {
4137#ifndef DEAL_II_WITH_MPI
4138
4139 (void)level;
4140 (void)new_numbers;
4141
4143 return NumberCache();
4144#else
4145
4147 *triangulation =
4148 (dynamic_cast<
4150 const_cast<::Triangulation<dim, spacedim> *>(
4151 &dof_handler->get_triangulation())));
4152 Assert(triangulation != nullptr, ExcInternalError());
4153
4154 // This code is very close to the respective code in renumber_dofs,
4155 // with the difference that we work on different entities with
4156 // different objects.
4157 const IndexSet &owned_dofs = dof_handler->locally_owned_mg_dofs(level);
4158 AssertDimension(new_numbers.size(), owned_dofs.n_elements());
4159
4160 const bool locally_owned_set_changes =
4161 std::any_of(new_numbers.cbegin(),
4162 new_numbers.cend(),
4163 [&owned_dofs](const types::global_dof_index i) {
4164 return owned_dofs.is_element(i) == false;
4165 });
4166
4167 IndexSet my_locally_owned_new_dof_indices = owned_dofs;
4168 if (locally_owned_set_changes && owned_dofs.n_elements() > 0)
4169 {
4170 std::vector<::types::global_dof_index> new_numbers_sorted =
4171 new_numbers;
4172 std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());
4173
4174 my_locally_owned_new_dof_indices =
4175 IndexSet(dof_handler->n_dofs(level));
4176 my_locally_owned_new_dof_indices.add_indices(
4177 new_numbers_sorted.begin(), new_numbers_sorted.end());
4178 my_locally_owned_new_dof_indices.compress();
4179
4180 Assert(my_locally_owned_new_dof_indices.n_elements() ==
4181 new_numbers.size(),
4183 }
4184
4185 // delete all knowledge of DoF indices that are not locally
4186 // owned
4187 for (auto cell : dof_handler->cell_iterators_on_level(level))
4188 if (cell->is_ghost_on_level())
4189 {
4191 *cell,
4192 std::make_tuple(),
4193 0,
4195 MGDoFIndexProcessor<dim, spacedim>(cell->level()),
4196 [&owned_dofs](auto &stored_index, auto) {
4197 if ((stored_index != numbers::invalid_dof_index) &&
4198 (!owned_dofs.is_element(stored_index)))
4199 stored_index = numbers::invalid_dof_index;
4200 },
4201 true);
4202 }
4203
4204 // renumber. Skip when there is nothing to do because we own no DoF.
4205 if (level < triangulation->n_levels() && owned_dofs.n_elements() > 0)
4207 new_numbers, owned_dofs, *dof_handler, level, false);
4208
4209 // communicate newly assigned DoF indices with other processors
4210 {
4211 std::vector<std::vector<bool>> cell_marked(triangulation->n_levels());
4212 for (unsigned int l = 0; l < triangulation->n_levels(); ++l)
4213 cell_marked[l].resize(triangulation->n_raw_cells(l));
4214 for (const auto &cell : dof_handler->cell_iterators_on_level(level))
4215 if (cell->is_ghost_on_level())
4216 cell_marked[cell->level()][cell->index()] = true;
4217
4218 communicate_mg_ghost_cells(*dof_handler, cell_marked);
4219
4220 communicate_mg_ghost_cells(*dof_handler, cell_marked);
4221 }
4222
4223 NumberCache number_cache;
4224 number_cache.locally_owned_dofs = my_locally_owned_new_dof_indices;
4225 number_cache.n_global_dofs = dof_handler->n_dofs(level);
4226 number_cache.n_locally_owned_dofs =
4227 number_cache.locally_owned_dofs.n_elements();
4228 return number_cache;
4229#endif
4230 }
4231 } // namespace Policy
4232 } // namespace DoFHandlerImplementation
4233} // namespace internal
4234
4235
4236
4237/*-------------- Explicit Instantiations -------------------------------*/
4238#include "dofs/dof_handler_policy.inst"
4239
4240
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:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#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 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:269
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
constexpr types::subdomain_id artificial_subdomain_id
Definition types.h:402
constexpr types::subdomain_id invalid_subdomain_id
Definition types.h:381
constexpr types::fe_index invalid_fe_index
Definition types.h:260
unsigned short int fe_index
Definition types.h:72
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 > &)