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