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