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