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