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