Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3110-g10dd77059b 2025-04-22 10:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
dof_renumbering.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17#include <deal.II/base/types.h>
19
21
26
27#include <deal.II/fe/fe.h>
30
31#include <deal.II/grid/tria.h>
33
36
41
43
46
47#define BOOST_BIND_GLOBAL_PLACEHOLDERS
48#include <boost/config.hpp>
49#include <boost/graph/adjacency_list.hpp>
50#include <boost/graph/bandwidth.hpp>
51#include <boost/graph/cuthill_mckee_ordering.hpp>
52#include <boost/graph/king_ordering.hpp>
53#include <boost/graph/minimum_degree_ordering.hpp>
54#include <boost/graph/properties.hpp>
55#include <boost/random.hpp>
56#include <boost/random/uniform_int_distribution.hpp>
57
58#undef BOOST_BIND_GLOBAL_PLACEHOLDERS
59
60#include <algorithm>
61#include <cmath>
62#include <functional>
63#include <map>
64#include <vector>
65
66
68
69namespace DoFRenumbering
70{
71 namespace boost
72 {
73 namespace boosttypes
74 {
75 using namespace ::boost;
76
77 using Graph = adjacency_list<vecS,
78 vecS,
79 undirectedS,
80 property<vertex_color_t,
81 default_color_type,
82 property<vertex_degree_t, int>>>;
83 using Vertex = graph_traits<Graph>::vertex_descriptor;
84 using size_type = graph_traits<Graph>::vertices_size_type;
85
86 using Pair = std::pair<size_type, size_type>;
87 } // namespace boosttypes
88
89
90 namespace internal
91 {
92 template <int dim, int spacedim>
93 void
95 const bool use_constraints,
96 boosttypes::Graph &graph,
97 boosttypes::property_map<boosttypes::Graph,
98 boosttypes::vertex_degree_t>::type
99 &graph_degree)
100 {
101 {
102 // create intermediate sparsity pattern (faster than directly
103 // submitting indices)
104 AffineConstraints<double> constraints;
105 if (use_constraints)
106 DoFTools::make_hanging_node_constraints(dof_handler, constraints);
107 constraints.close();
108 DynamicSparsityPattern dsp(dof_handler.n_dofs(),
109 dof_handler.n_dofs());
110 DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
111
112 // submit the entries to the boost graph
113 for (unsigned int row = 0; row < dsp.n_rows(); ++row)
114 for (unsigned int col = 0; col < dsp.row_length(row); ++col)
115 add_edge(row, dsp.column_number(row, col), graph);
116 }
117
118 boosttypes::graph_traits<boosttypes::Graph>::vertex_iterator ui, ui_end;
119
120 graph_degree = get(::boost::vertex_degree, graph);
121 for (::boost::tie(ui, ui_end) = vertices(graph); ui != ui_end; ++ui)
122 graph_degree[*ui] = degree(*ui, graph);
123 }
124 } // namespace internal
125
126
127 template <int dim, int spacedim>
128 void
130 const bool reversed_numbering,
131 const bool use_constraints)
132 {
133 std::vector<types::global_dof_index> renumbering(
134 dof_handler.n_dofs(), numbers::invalid_dof_index);
135 compute_Cuthill_McKee(renumbering,
136 dof_handler,
137 reversed_numbering,
138 use_constraints);
139
140 // actually perform renumbering;
141 // this is dimension specific and
142 // thus needs an own function
143 dof_handler.renumber_dofs(renumbering);
144 }
145
146
147 template <int dim, int spacedim>
148 void
149 compute_Cuthill_McKee(std::vector<types::global_dof_index> &new_dof_indices,
150 const DoFHandler<dim, spacedim> &dof_handler,
151 const bool reversed_numbering,
152 const bool use_constraints)
153 {
154 boosttypes::Graph graph(dof_handler.n_dofs());
155 boosttypes::property_map<boosttypes::Graph,
156 boosttypes::vertex_degree_t>::type graph_degree;
157
158 internal::create_graph(dof_handler, use_constraints, graph, graph_degree);
159
160 boosttypes::property_map<boosttypes::Graph,
161 boosttypes::vertex_index_t>::type index_map =
162 get(::boost::vertex_index, graph);
163
164
165 std::vector<boosttypes::Vertex> inv_perm(num_vertices(graph));
166
167 if (reversed_numbering == false)
168 ::boost::cuthill_mckee_ordering(graph,
169 inv_perm.rbegin(),
170 get(::boost::vertex_color, graph),
171 make_degree_map(graph));
172 else
173 ::boost::cuthill_mckee_ordering(graph,
174 inv_perm.begin(),
175 get(::boost::vertex_color, graph),
176 make_degree_map(graph));
177
178 for (boosttypes::size_type c = 0; c != inv_perm.size(); ++c)
179 new_dof_indices[index_map[inv_perm[c]]] = c;
180
181 Assert(std::find(new_dof_indices.begin(),
182 new_dof_indices.end(),
183 numbers::invalid_dof_index) == new_dof_indices.end(),
185 }
186
187
188
189 template <int dim, int spacedim>
190 void
192 const bool reversed_numbering,
193 const bool use_constraints)
194 {
195 std::vector<types::global_dof_index> renumbering(
196 dof_handler.n_dofs(), numbers::invalid_dof_index);
197 compute_king_ordering(renumbering,
198 dof_handler,
199 reversed_numbering,
200 use_constraints);
201
202 // actually perform renumbering;
203 // this is dimension specific and
204 // thus needs an own function
205 dof_handler.renumber_dofs(renumbering);
206 }
207
208
209 template <int dim, int spacedim>
210 void
211 compute_king_ordering(std::vector<types::global_dof_index> &new_dof_indices,
212 const DoFHandler<dim, spacedim> &dof_handler,
213 const bool reversed_numbering,
214 const bool use_constraints)
215 {
216 boosttypes::Graph graph(dof_handler.n_dofs());
217 boosttypes::property_map<boosttypes::Graph,
218 boosttypes::vertex_degree_t>::type graph_degree;
219
220 internal::create_graph(dof_handler, use_constraints, graph, graph_degree);
221
222 boosttypes::property_map<boosttypes::Graph,
223 boosttypes::vertex_index_t>::type index_map =
224 get(::boost::vertex_index, graph);
225
226
227 std::vector<boosttypes::Vertex> inv_perm(num_vertices(graph));
228
229 if (reversed_numbering == false)
230 ::boost::king_ordering(graph, inv_perm.rbegin());
231 else
232 ::boost::king_ordering(graph, inv_perm.begin());
233
234 for (boosttypes::size_type c = 0; c != inv_perm.size(); ++c)
235 new_dof_indices[index_map[inv_perm[c]]] = c;
236
237 Assert(std::find(new_dof_indices.begin(),
238 new_dof_indices.end(),
239 numbers::invalid_dof_index) == new_dof_indices.end(),
241 }
242
243
244
245 template <int dim, int spacedim>
246 void
248 const bool reversed_numbering,
249 const bool use_constraints)
250 {
251 std::vector<types::global_dof_index> renumbering(
252 dof_handler.n_dofs(), numbers::invalid_dof_index);
253 compute_minimum_degree(renumbering,
254 dof_handler,
255 reversed_numbering,
256 use_constraints);
257
258 // actually perform renumbering;
259 // this is dimension specific and
260 // thus needs an own function
261 dof_handler.renumber_dofs(renumbering);
262 }
263
264
265 template <int dim, int spacedim>
266 void
268 std::vector<types::global_dof_index> &new_dof_indices,
269 const DoFHandler<dim, spacedim> &dof_handler,
270 const bool reversed_numbering,
271 const bool use_constraints)
272 {
273 (void)use_constraints;
274 Assert(use_constraints == false, ExcNotImplemented());
275
276 // the following code is pretty
277 // much a verbatim copy of the
278 // sample code for the
279 // minimum_degree_ordering manual
280 // page from the BOOST Graph
281 // Library
282 using namespace ::boost;
283
284 int delta = 0;
285
286 // must be BGL directed graph now
287 using Graph = adjacency_list<vecS, vecS, directedS>;
288
289 int n = dof_handler.n_dofs();
290
291 Graph G(n);
292
293 std::vector<::types::global_dof_index> dofs_on_this_cell;
294
295 for (const auto &cell : dof_handler.active_cell_iterators())
296 {
297 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
298
299 dofs_on_this_cell.resize(dofs_per_cell);
300
301 cell->get_active_or_mg_dof_indices(dofs_on_this_cell);
302 for (unsigned int i = 0; i < dofs_per_cell; ++i)
303 for (unsigned int j = 0; j < dofs_per_cell; ++j)
304 if (dofs_on_this_cell[i] > dofs_on_this_cell[j])
305 {
306 add_edge(dofs_on_this_cell[i], dofs_on_this_cell[j], G);
307 add_edge(dofs_on_this_cell[j], dofs_on_this_cell[i], G);
308 }
309 }
310
311
312 using Vector = std::vector<int>;
313
314
315 Vector inverse_perm(n, 0);
316
317 Vector perm(n, 0);
318
319
320 Vector supernode_sizes(n, 1);
321 // init has to be 1
322
323 ::boost::property_map<Graph, vertex_index_t>::type id =
324 get(vertex_index, G);
325
326
327 Vector degree(n, 0);
328
329
330 minimum_degree_ordering(
331 G,
332 make_iterator_property_map(degree.begin(), id, degree[0]),
333 inverse_perm.data(),
334 perm.data(),
335 make_iterator_property_map(supernode_sizes.begin(),
336 id,
337 supernode_sizes[0]),
338 delta,
339 id);
340
341
342 for (int i = 0; i < n; ++i)
343 {
344 Assert(std::find(perm.begin(), perm.end(), i) != perm.end(),
346 Assert(std::find(inverse_perm.begin(), inverse_perm.end(), i) !=
347 inverse_perm.end(),
349 Assert(inverse_perm[perm[i]] == i, ExcInternalError());
350 }
351
352 if (reversed_numbering == true)
353 std::copy(perm.begin(), perm.end(), new_dof_indices.begin());
354 else
355 std::copy(inverse_perm.begin(),
356 inverse_perm.end(),
357 new_dof_indices.begin());
358 }
359
360 } // namespace boost
361
362
363
364 template <int dim, int spacedim>
365 void
367 const bool reversed_numbering,
368 const bool use_constraints,
369 const std::vector<types::global_dof_index> &starting_indices)
370 {
371 std::vector<types::global_dof_index> renumbering(
372 dof_handler.locally_owned_dofs().n_elements(),
374 compute_Cuthill_McKee(renumbering,
375 dof_handler,
376 reversed_numbering,
377 use_constraints,
378 starting_indices);
379
380 // actually perform renumbering;
381 // this is dimension specific and
382 // thus needs an own function
383 dof_handler.renumber_dofs(renumbering);
384 }
385
386
387
388 template <int dim, int spacedim>
389 void
391 std::vector<types::global_dof_index> &new_indices,
392 const DoFHandler<dim, spacedim> &dof_handler,
393 const bool reversed_numbering,
394 const bool use_constraints,
395 const std::vector<types::global_dof_index> &starting_indices,
396 const unsigned int level)
397 {
398 const bool reorder_level_dofs =
399 (level == numbers::invalid_unsigned_int) ? false : true;
400
401 // see if there is anything to do at all or whether we can skip the work on
402 // this processor
403 if (dof_handler.locally_owned_dofs().n_elements() == 0)
404 {
405 Assert(new_indices.empty(), ExcInternalError());
406 return;
407 }
408
409 // make the connection graph
410 //
411 // note that if constraints are not requested, then the 'constraints'
412 // object will be empty and using it has no effect
413 if (reorder_level_dofs == true)
416
417 const IndexSet locally_relevant_dofs =
418 (reorder_level_dofs == false ?
421 const IndexSet &locally_owned_dofs =
422 (reorder_level_dofs == false ? dof_handler.locally_owned_dofs() :
423 dof_handler.locally_owned_mg_dofs(level));
424
425 AffineConstraints<double> constraints;
426 if (use_constraints)
427 {
428 // reordering with constraints is not yet implemented on a level basis
429 Assert(reorder_level_dofs == false, ExcNotImplemented());
430
431 constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
432 DoFTools::make_hanging_node_constraints(dof_handler, constraints);
433 }
434 constraints.close();
435
436 // see if we can get away with the sequential algorithm
437 if (locally_owned_dofs.n_elements() == locally_owned_dofs.size())
438 {
439 AssertDimension(new_indices.size(), locally_owned_dofs.n_elements());
440
441 DynamicSparsityPattern dsp(locally_owned_dofs.size(),
442 locally_owned_dofs.size());
443 if (reorder_level_dofs == false)
444 {
445 DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
446 }
447 else
448 {
449 MGTools::make_sparsity_pattern(dof_handler, dsp, level);
450 }
451
453 new_indices,
454 starting_indices);
455 if (reversed_numbering)
456 new_indices = Utilities::reverse_permutation(new_indices);
457 }
458 else
459 {
460 // we are in the parallel case where we need to work in the
461 // local index space, i.e., the locally owned part of the
462 // sparsity pattern.
463 //
464 // first figure out whether the user only gave us starting
465 // indices that are locally owned, or that are only locally
466 // relevant. in the process, also check that all indices
467 // really belong to at least the locally relevant ones
468 const IndexSet locally_active_dofs =
469 (reorder_level_dofs == false ?
472
473 bool needs_locally_active = false;
474 for (const auto starting_index : starting_indices)
475 {
476 if ((needs_locally_active ==
477 /* previously already set to */ true) ||
478 (locally_owned_dofs.is_element(starting_index) == false))
479 {
480 Assert(
481 locally_active_dofs.is_element(starting_index),
483 "You specified global degree of freedom " +
484 std::to_string(starting_index) +
485 " as a starting index, but this index is not among the "
486 "locally active ones on this processor, as required "
487 "for this function."));
488 needs_locally_active = true;
489 }
490 }
491
492 const IndexSet index_set_to_use =
493 (needs_locally_active ? locally_active_dofs : locally_owned_dofs);
494
495 // if this process doesn't own any DoFs (on this level), there is
496 // nothing to do
497 if (index_set_to_use.n_elements() == 0)
498 return;
499
500 // then create first the global sparsity pattern, and then the local
501 // sparsity pattern from the global one by transferring its indices to
502 // processor-local (locally owned or locally active) index space
503 DynamicSparsityPattern dsp(index_set_to_use.size(),
504 index_set_to_use.size(),
505 index_set_to_use);
506 if (reorder_level_dofs == false)
507 {
508 DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
509 }
510 else
511 {
512 MGTools::make_sparsity_pattern(dof_handler, dsp, level);
513 }
514
515 DynamicSparsityPattern local_sparsity(index_set_to_use.n_elements(),
516 index_set_to_use.n_elements());
517 std::vector<types::global_dof_index> row_entries;
518 for (unsigned int i = 0; i < index_set_to_use.n_elements(); ++i)
519 {
520 const types::global_dof_index row =
521 index_set_to_use.nth_index_in_set(i);
522 const unsigned int row_length = dsp.row_length(row);
523 row_entries.clear();
524 for (unsigned int j = 0; j < row_length; ++j)
525 {
526 const unsigned int col = dsp.column_number(row, j);
527 if (col != row && index_set_to_use.is_element(col))
528 row_entries.push_back(index_set_to_use.index_within_set(col));
529 }
530 local_sparsity.add_entries(i,
531 row_entries.begin(),
532 row_entries.end(),
533 true);
534 }
535
536 // translate starting indices from global to local indices
537 std::vector<types::global_dof_index> local_starting_indices(
538 starting_indices.size());
539 for (unsigned int i = 0; i < starting_indices.size(); ++i)
540 local_starting_indices[i] =
541 index_set_to_use.index_within_set(starting_indices[i]);
542
543 // then do the renumbering on the locally owned portion
544 AssertDimension(new_indices.size(), locally_owned_dofs.n_elements());
545 std::vector<types::global_dof_index> my_new_indices(
546 index_set_to_use.n_elements());
548 my_new_indices,
549 local_starting_indices);
550 if (reversed_numbering)
551 my_new_indices = Utilities::reverse_permutation(my_new_indices);
552
553 // now that we have a re-enumeration of all DoFs, we need to throw
554 // out the ones that are not locally owned in case we have worked
555 // with the locally active ones. that's because the renumbering
556 // functions only want new indices for the locally owned DoFs (other
557 // processors are responsible for renumbering the ones that are
558 // on cell interfaces)
559 if (needs_locally_active == true)
560 {
561 // first step: figure out which DoF indices to eliminate
562 IndexSet active_but_not_owned_dofs = locally_active_dofs;
563 active_but_not_owned_dofs.subtract_set(locally_owned_dofs);
564
565 std::set<types::global_dof_index> erase_these_indices;
566 for (const auto p : active_but_not_owned_dofs)
567 {
568 const auto index = index_set_to_use.index_within_set(p);
569 Assert(index < index_set_to_use.n_elements(),
571 erase_these_indices.insert(my_new_indices[index]);
572 my_new_indices[index] = numbers::invalid_dof_index;
573 }
574 Assert(erase_these_indices.size() ==
575 active_but_not_owned_dofs.n_elements(),
577 Assert(static_cast<unsigned int>(
578 std::count(my_new_indices.begin(),
579 my_new_indices.end(),
581 active_but_not_owned_dofs.n_elements(),
583
584 // then compute a renumbering of the remaining ones
585 std::vector<types::global_dof_index> translate_indices(
586 my_new_indices.size());
587 {
588 std::set<types::global_dof_index>::const_iterator
589 next_erased_index = erase_these_indices.begin();
590 types::global_dof_index next_new_index = 0;
591 for (unsigned int i = 0; i < translate_indices.size(); ++i)
592 if ((next_erased_index != erase_these_indices.end()) &&
593 (*next_erased_index == i))
594 {
595 translate_indices[i] = numbers::invalid_dof_index;
596 ++next_erased_index;
597 }
598 else
599 {
600 translate_indices[i] = next_new_index;
601 ++next_new_index;
602 }
603 Assert(next_new_index == locally_owned_dofs.n_elements(),
605 }
606
607 // and then do the renumbering of the result of the
608 // Cuthill-McKee algorithm above, right into the output array
609 new_indices.clear();
610 new_indices.reserve(locally_owned_dofs.n_elements());
611 for (const auto &p : my_new_indices)
613 {
614 Assert(translate_indices[p] != numbers::invalid_dof_index,
616 new_indices.push_back(translate_indices[p]);
617 }
618 Assert(new_indices.size() == locally_owned_dofs.n_elements(),
620 }
621 else
622 new_indices = std::move(my_new_indices);
623
624 // convert indices back to global index space. in both of the branches
625 // above, we ended up with new_indices only containing the local
626 // indices of the locally-owned DoFs. so that's where we get the
627 // indices
628 for (types::global_dof_index &new_index : new_indices)
629 new_index = locally_owned_dofs.nth_index_in_set(new_index);
630 }
631 }
632
633
634
635 template <int dim, int spacedim>
636 void
638 const unsigned int level,
639 const bool reversed_numbering,
640 const std::vector<types::global_dof_index> &starting_indices)
641 {
644
645 std::vector<types::global_dof_index> new_indices(
648
649 compute_Cuthill_McKee(new_indices,
650 dof_handler,
651 reversed_numbering,
652 false,
653 starting_indices,
654 level);
655
656 // actually perform renumbering;
657 // this is dimension specific and
658 // thus needs an own function
659 dof_handler.renumber_dofs(level, new_indices);
660 }
661
662
663
664 template <int dim, int spacedim>
665 void
667 const std::vector<unsigned int> &component_order_arg)
668 {
669 std::vector<types::global_dof_index> renumbering(
671
672 const types::global_dof_index result =
673 compute_component_wise<dim, spacedim>(renumbering,
674 dof_handler.begin_active(),
675 dof_handler.end(),
676 component_order_arg,
677 false);
678 (void)result;
679
680 // If we don't have a renumbering (i.e., when there is 1 component) then
681 // return
682 if (Utilities::MPI::max(renumbering.size(),
683 dof_handler.get_mpi_communicator()) == 0)
684 return;
685
686 // verify that the last numbered
687 // degree of freedom is either
688 // equal to the number of degrees
689 // of freedom in total (the
690 // sequential case) or in the
691 // distributed case at least
692 // makes sense
693 Assert((result == dof_handler.n_locally_owned_dofs()) ||
694 ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
695 (result <= dof_handler.n_dofs())),
697
698 dof_handler.renumber_dofs(renumbering);
699 }
700
701
702
703 template <int dim, int spacedim>
704 void
706 const unsigned int level,
707 const std::vector<unsigned int> &component_order_arg)
708 {
711
712 std::vector<types::global_dof_index> renumbering(
715
717 dof_handler.begin(level);
719 dof_handler.end(level);
720
721 const types::global_dof_index result =
722 compute_component_wise<dim, spacedim>(
723 renumbering, start, end, component_order_arg, true);
724 (void)result;
725
726 // If we don't have a renumbering (i.e., when there is 1 component) then
727 // return
728 if (Utilities::MPI::max(renumbering.size(),
729 dof_handler.get_mpi_communicator()) == 0)
730 return;
731
732 // verify that the last numbered
733 // degree of freedom is either
734 // equal to the number of degrees
735 // of freedom in total (the
736 // sequential case) or in the
737 // distributed case at least
738 // makes sense
739 Assert((result == dof_handler.locally_owned_mg_dofs(level).n_elements()) ||
740 ((dof_handler.locally_owned_mg_dofs(level).n_elements() <
741 dof_handler.n_dofs(level)) &&
742 (result <= dof_handler.n_dofs(level))),
744
745 dof_handler.renumber_dofs(level, renumbering);
746 }
747
748
749
750 template <int dim, int spacedim, typename CellIterator>
752 compute_component_wise(std::vector<types::global_dof_index> &new_indices,
753 const CellIterator &start,
755 const std::vector<unsigned int> &component_order_arg,
756 const bool is_level_operation)
757 {
758 const hp::FECollection<dim, spacedim> &fe_collection =
759 start->get_dof_handler().get_fe_collection();
760
761 // do nothing if the FE has only
762 // one component
763 if (fe_collection.n_components() == 1)
764 {
765 new_indices.resize(0);
766 return 0;
767 }
768
769 // Get a reference to the set of dofs. Note that we assume that all cells
770 // are assumed to be on the same level, otherwise the operation doesn't make
771 // much sense (we will assert this below).
772 const IndexSet &locally_owned_dofs =
773 is_level_operation ?
774 start->get_dof_handler().locally_owned_mg_dofs(start->level()) :
775 start->get_dof_handler().locally_owned_dofs();
776
777 // Copy last argument into a
778 // writable vector.
779 std::vector<unsigned int> component_order(component_order_arg);
780 // If the last argument was an
781 // empty vector, set up things to
782 // store components in the order
783 // found in the system.
784 if (component_order.empty())
785 for (unsigned int i = 0; i < fe_collection.n_components(); ++i)
786 component_order.push_back(i);
787
788 Assert(component_order.size() == fe_collection.n_components(),
789 ExcDimensionMismatch(component_order.size(),
790 fe_collection.n_components()));
791
792 for (const unsigned int component : component_order)
793 {
794 (void)component;
795 AssertIndexRange(component, fe_collection.n_components());
796 }
797
798 // vector to hold the dof indices on
799 // the cell we visit at a time
800 std::vector<types::global_dof_index> local_dof_indices;
801
802 // prebuilt list to which component
803 // a given dof on a cell
804 // should go. note that we get into
805 // trouble here if the shape
806 // function is not primitive, since
807 // then there is no single vector
808 // component to which it
809 // belongs. in this case, assign it
810 // to the first vector component to
811 // which it belongs
812 std::vector<std::vector<unsigned int>> component_list(fe_collection.size());
813 for (unsigned int f = 0; f < fe_collection.size(); ++f)
814 {
815 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
816 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
817 component_list[f].resize(dofs_per_cell);
818 for (unsigned int i = 0; i < dofs_per_cell; ++i)
819 if (fe.is_primitive(i))
820 component_list[f][i] =
821 component_order[fe.system_to_component_index(i).first];
822 else
823 {
824 const unsigned int comp =
826
827 // then associate this degree
828 // of freedom with this
829 // component
830 component_list[f][i] = component_order[comp];
831 }
832 }
833
834 // set up a map where for each
835 // component the respective degrees
836 // of freedom are collected.
837 //
838 // note that this map is sorted by
839 // component but that within each
840 // component it is NOT sorted by
841 // dof index. note also that some
842 // dof indices are entered
843 // multiply, so we will have to
844 // take care of that
845 std::vector<std::vector<types::global_dof_index>> component_to_dof_map(
846 fe_collection.n_components());
847 for (CellIterator cell = start; cell != end; ++cell)
848 {
849 if (is_level_operation)
850 {
851 // we are dealing with mg dofs, skip foreign level cells:
852 if (!cell->is_locally_owned_on_level())
853 continue;
854 }
855 else
856 {
857 // we are dealing with active dofs, skip the loop if not locally
858 // owned:
859 if (!cell->is_locally_owned())
860 continue;
861 }
862
863 if (is_level_operation)
864 Assert(
865 cell->level() == start->level(),
867 "Multigrid renumbering in compute_component_wise() needs to be applied to a single level!"));
868
869 // on each cell: get dof indices
870 // and insert them into the global
871 // list using their component
872 const types::fe_index fe_index = cell->active_fe_index();
873 const unsigned int dofs_per_cell =
874 fe_collection[fe_index].n_dofs_per_cell();
875 local_dof_indices.resize(dofs_per_cell);
876 cell->get_active_or_mg_dof_indices(local_dof_indices);
877
878 for (unsigned int i = 0; i < dofs_per_cell; ++i)
879 if (locally_owned_dofs.is_element(local_dof_indices[i]))
880 component_to_dof_map[component_list[fe_index][i]].push_back(
881 local_dof_indices[i]);
882 }
883
884 // now we've got all indices sorted
885 // into buckets labeled by their
886 // target component number. we've
887 // only got to traverse this list
888 // and assign the new indices
889 //
890 // however, we first want to sort
891 // the indices entered into the
892 // buckets to preserve the order
893 // within each component and during
894 // this also remove duplicate
895 // entries
896 //
897 // note that we no longer have to
898 // care about non-primitive shape
899 // functions since the buckets
900 // corresponding to the second and
901 // following vector components of a
902 // non-primitive FE will simply be
903 // empty, everything being shoved
904 // into the first one. The same
905 // holds if several components were
906 // joined into a single target.
907 for (unsigned int component = 0; component < fe_collection.n_components();
908 ++component)
909 {
910 std::sort(component_to_dof_map[component].begin(),
911 component_to_dof_map[component].end());
912 component_to_dof_map[component].erase(
913 std::unique(component_to_dof_map[component].begin(),
914 component_to_dof_map[component].end()),
915 component_to_dof_map[component].end());
916 }
917
918 // calculate the number of locally owned
919 // DoFs per bucket
920 const unsigned int n_buckets = fe_collection.n_components();
921 std::vector<types::global_dof_index> shifts(n_buckets);
922
924 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
925 &start->get_dof_handler().get_triangulation())))
926 {
927#ifdef DEAL_II_WITH_MPI
928 std::vector<types::global_dof_index> local_dof_count(n_buckets);
929
930 for (unsigned int c = 0; c < n_buckets; ++c)
931 local_dof_count[c] = component_to_dof_map[c].size();
932
933 std::vector<types::global_dof_index> prefix_dof_count(n_buckets);
934 const int ierr = MPI_Exscan(
935 local_dof_count.data(),
936 prefix_dof_count.data(),
937 n_buckets,
938 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
939 MPI_SUM,
940 tria->get_mpi_communicator());
941 AssertThrowMPI(ierr);
942
943 std::vector<types::global_dof_index> global_dof_count(n_buckets);
944 Utilities::MPI::sum(local_dof_count,
945 tria->get_mpi_communicator(),
946 global_dof_count);
947
948 // calculate shifts
949 types::global_dof_index cumulated = 0;
950 for (unsigned int c = 0; c < n_buckets; ++c)
951 {
952 shifts[c] = prefix_dof_count[c] + cumulated;
953 cumulated += global_dof_count[c];
954 }
955#else
956 (void)tria;
958#endif
959 }
960 else
961 {
962 shifts[0] = 0;
963 for (unsigned int c = 1; c < fe_collection.n_components(); ++c)
964 shifts[c] = shifts[c - 1] + component_to_dof_map[c - 1].size();
965 }
966
967
968
969 // now concatenate all the
970 // components in the order the user
971 // desired to see
972 types::global_dof_index next_free_index = 0;
973 for (unsigned int component = 0; component < fe_collection.n_components();
974 ++component)
975 {
976 next_free_index = shifts[component];
977
978 for (const types::global_dof_index dof_index :
979 component_to_dof_map[component])
980 {
981 Assert(locally_owned_dofs.index_within_set(dof_index) <
982 new_indices.size(),
984 new_indices[locally_owned_dofs.index_within_set(dof_index)] =
985 next_free_index;
986
987 ++next_free_index;
988 }
989 }
990
991 return next_free_index;
992 }
993
994
995
996 template <int dim, int spacedim>
997 void
999 {
1000 std::vector<types::global_dof_index> renumbering(
1002
1004 dim,
1005 spacedim,
1008 renumbering, dof_handler.begin_active(), dof_handler.end(), false);
1009 if (result == 0)
1010 return;
1011
1012 // verify that the last numbered
1013 // degree of freedom is either
1014 // equal to the number of degrees
1015 // of freedom in total (the
1016 // sequential case) or in the
1017 // distributed case at least
1018 // makes sense
1019 Assert((result == dof_handler.n_locally_owned_dofs()) ||
1020 ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
1021 (result <= dof_handler.n_dofs())),
1023
1024 dof_handler.renumber_dofs(renumbering);
1025 }
1026
1027
1028
1029 template <int dim, int spacedim>
1030 void
1031 block_wise(DoFHandler<dim, spacedim> &dof_handler, const unsigned int level)
1032 {
1035
1036 std::vector<types::global_dof_index> renumbering(
1037 dof_handler.locally_owned_mg_dofs(level).n_elements(),
1039
1041 dof_handler.begin(level);
1043 dof_handler.end(level);
1044
1046 dim,
1047 spacedim,
1050 start,
1051 end,
1052 true);
1053 (void)result;
1054
1055 Assert(result == 0 || result == dof_handler.n_dofs(level),
1057
1058 if (Utilities::MPI::max(renumbering.size(),
1059 dof_handler.get_mpi_communicator()) > 0)
1060 dof_handler.renumber_dofs(level, renumbering);
1061 }
1062
1063
1064
1065 template <int dim, int spacedim, class IteratorType, class EndIteratorType>
1067 compute_block_wise(std::vector<types::global_dof_index> &new_indices,
1068 const IteratorType &start,
1069 const EndIteratorType &end,
1070 const bool is_level_operation)
1071 {
1072 const hp::FECollection<dim, spacedim> &fe_collection =
1073 start->get_dof_handler().get_fe_collection();
1074
1075 // do nothing if the FE has only
1076 // one component
1077 if (fe_collection.n_blocks() == 1)
1078 {
1079 new_indices.resize(0);
1080 return 0;
1081 }
1082
1083 // Get a reference to the set of dofs. Note that we assume that all cells
1084 // are assumed to be on the same level, otherwise the operation doesn't make
1085 // much sense (we will assert this below).
1086 const IndexSet &locally_owned_dofs =
1087 is_level_operation ?
1088 start->get_dof_handler().locally_owned_mg_dofs(start->level()) :
1089 start->get_dof_handler().locally_owned_dofs();
1090
1091 // vector to hold the dof indices on
1092 // the cell we visit at a time
1093 std::vector<types::global_dof_index> local_dof_indices;
1094
1095 // prebuilt list to which block
1096 // a given dof on a cell
1097 // should go.
1098 std::vector<std::vector<types::global_dof_index>> block_list(
1099 fe_collection.size());
1100 for (unsigned int f = 0; f < fe_collection.size(); ++f)
1101 {
1102 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
1103 block_list[f].resize(fe.n_dofs_per_cell());
1104 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
1105 block_list[f][i] = fe.system_to_block_index(i).first;
1106 }
1107
1108 // set up a map where for each
1109 // block the respective degrees
1110 // of freedom are collected.
1111 //
1112 // note that this map is sorted by
1113 // block but that within each
1114 // block it is NOT sorted by
1115 // dof index. note also that some
1116 // dof indices are entered
1117 // multiply, so we will have to
1118 // take care of that
1119 std::vector<std::vector<types::global_dof_index>> block_to_dof_map(
1120 fe_collection.n_blocks());
1121 for (IteratorType cell = start; cell != end; ++cell)
1122 {
1123 if (is_level_operation)
1124 {
1125 // we are dealing with mg dofs, skip foreign level cells:
1126 if (!cell->is_locally_owned_on_level())
1127 continue;
1128 }
1129 else
1130 {
1131 // we are dealing with active dofs, skip the loop if not locally
1132 // owned:
1133 if (!cell->is_locally_owned())
1134 continue;
1135 }
1136
1137 if (is_level_operation)
1138 Assert(
1139 cell->level() == start->level(),
1140 ExcMessage(
1141 "Multigrid renumbering in compute_block_wise() needs to be applied to a single level!"));
1142
1143 // on each cell: get dof indices
1144 // and insert them into the global
1145 // list using their component
1146 const types::fe_index fe_index = cell->active_fe_index();
1147 const unsigned int dofs_per_cell =
1148 fe_collection[fe_index].n_dofs_per_cell();
1149 local_dof_indices.resize(dofs_per_cell);
1150 cell->get_active_or_mg_dof_indices(local_dof_indices);
1151
1152 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1153 if (locally_owned_dofs.is_element(local_dof_indices[i]))
1154 block_to_dof_map[block_list[fe_index][i]].push_back(
1155 local_dof_indices[i]);
1156 }
1157
1158 // now we've got all indices sorted
1159 // into buckets labeled by their
1160 // target block number. we've
1161 // only got to traverse this list
1162 // and assign the new indices
1163 //
1164 // however, we first want to sort
1165 // the indices entered into the
1166 // buckets to preserve the order
1167 // within each component and during
1168 // this also remove duplicate
1169 // entries
1170 for (unsigned int block = 0; block < fe_collection.n_blocks(); ++block)
1171 {
1172 std::sort(block_to_dof_map[block].begin(),
1173 block_to_dof_map[block].end());
1174 block_to_dof_map[block].erase(
1175 std::unique(block_to_dof_map[block].begin(),
1176 block_to_dof_map[block].end()),
1177 block_to_dof_map[block].end());
1178 }
1179
1180 // calculate the number of locally owned
1181 // DoFs per bucket
1182 const unsigned int n_buckets = fe_collection.n_blocks();
1183 std::vector<types::global_dof_index> shifts(n_buckets);
1184
1186 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1187 &start->get_dof_handler().get_triangulation())))
1188 {
1189#ifdef DEAL_II_WITH_MPI
1190 std::vector<types::global_dof_index> local_dof_count(n_buckets);
1191
1192 for (unsigned int c = 0; c < n_buckets; ++c)
1193 local_dof_count[c] = block_to_dof_map[c].size();
1194
1195 std::vector<types::global_dof_index> prefix_dof_count(n_buckets);
1196 const int ierr = MPI_Exscan(
1197 local_dof_count.data(),
1198 prefix_dof_count.data(),
1199 n_buckets,
1200 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
1201 MPI_SUM,
1202 tria->get_mpi_communicator());
1203 AssertThrowMPI(ierr);
1204
1205 std::vector<types::global_dof_index> global_dof_count(n_buckets);
1206 Utilities::MPI::sum(local_dof_count,
1207 tria->get_mpi_communicator(),
1208 global_dof_count);
1209
1210 // calculate shifts
1211 types::global_dof_index cumulated = 0;
1212 for (unsigned int c = 0; c < n_buckets; ++c)
1213 {
1214 shifts[c] = prefix_dof_count[c] + cumulated;
1215 cumulated += global_dof_count[c];
1216 }
1217#else
1218 (void)tria;
1220#endif
1221 }
1222 else
1223 {
1224 shifts[0] = 0;
1225 for (unsigned int c = 1; c < fe_collection.n_blocks(); ++c)
1226 shifts[c] = shifts[c - 1] + block_to_dof_map[c - 1].size();
1227 }
1228
1229
1230
1231 // now concatenate all the
1232 // components in the order the user
1233 // desired to see
1234 types::global_dof_index next_free_index = 0;
1235 for (unsigned int block = 0; block < fe_collection.n_blocks(); ++block)
1236 {
1237 const typename std::vector<types::global_dof_index>::const_iterator
1238 begin_of_component = block_to_dof_map[block].begin(),
1239 end_of_component = block_to_dof_map[block].end();
1240
1241 next_free_index = shifts[block];
1242
1243 for (typename std::vector<types::global_dof_index>::const_iterator
1244 dof_index = begin_of_component;
1245 dof_index != end_of_component;
1246 ++dof_index)
1247 {
1248 Assert(locally_owned_dofs.index_within_set(*dof_index) <
1249 new_indices.size(),
1251 new_indices[locally_owned_dofs.index_within_set(*dof_index)] =
1252 next_free_index++;
1253 }
1254 }
1255
1256 return next_free_index;
1257 }
1258
1259
1260
1261 namespace
1262 {
1263 // Helper function for DoFRenumbering::hierarchical(). This function
1264 // recurses into the given cell or, if that should be an active (terminal)
1265 // cell, renumbers DoF indices on it. The function starts renumbering with
1266 // 'next_free_dof_index' and returns the first still unused DoF index at the
1267 // end of its operation.
1268 template <int dim, typename CellIteratorType>
1270 compute_hierarchical_recursive(
1271 const types::global_dof_index next_free_dof_offset,
1272 const types::global_dof_index my_starting_index,
1273 const CellIteratorType &cell,
1274 const IndexSet &locally_owned_dof_indices,
1275 std::vector<types::global_dof_index> &new_indices)
1276 {
1277 types::global_dof_index current_next_free_dof_offset =
1278 next_free_dof_offset;
1279
1280 if (cell->has_children())
1281 {
1282 // recursion
1283 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1284 ++c)
1285 current_next_free_dof_offset =
1286 compute_hierarchical_recursive<dim>(current_next_free_dof_offset,
1287 my_starting_index,
1288 cell->child(c),
1289 locally_owned_dof_indices,
1290 new_indices);
1291 }
1292 else
1293 {
1294 // this is a terminal cell. we need to renumber its DoF indices. there
1295 // are now three cases to decide:
1296 // - this is a sequential triangulation: we can just go ahead and
1297 // number
1298 // the DoFs in the order in which we encounter cells. in this case,
1299 // all cells are actually locally owned
1300 // - if this is a parallel::distributed::Triangulation, then we only
1301 // need to work on the locally owned cells since they contain
1302 // all locally owned DoFs.
1303 // - if this is a parallel::shared::Triangulation, then the same
1304 // applies
1305 //
1306 // in all cases, each processor starts new indices so that we get
1307 // a consecutive numbering on each processor, and disjoint ownership
1308 // of the global range of DoF indices
1309 if (cell->is_locally_owned())
1310 {
1311 // first get the existing DoF indices
1312 const unsigned int dofs_per_cell =
1313 cell->get_fe().n_dofs_per_cell();
1314 std::vector<types::global_dof_index> local_dof_indices(
1315 dofs_per_cell);
1316 cell->get_dof_indices(local_dof_indices);
1317
1318 // then loop over the existing DoF indices on this cell
1319 // and see whether it has already been re-numbered (it
1320 // may have been on a face or vertex to a neighboring
1321 // cell that we have encountered before). if not,
1322 // give it a new number and store that number in the
1323 // output array (new_indices)
1324 //
1325 // if this is a parallel triangulation and a DoF index is
1326 // not locally owned, then don't touch it. since
1327 // we don't actually *change* DoF indices (just record new
1328 // numbers in an array), we don't need to worry about
1329 // the decision whether a DoF is locally owned or not changing
1330 // as we progress in renumbering DoFs -- all adjacent cells
1331 // will always agree that a DoF is locally owned or not.
1332 // that said, the first cell to encounter a locally owned DoF
1333 // gets to number it, so the order in which we traverse cells
1334 // matters
1335 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1336 if (locally_owned_dof_indices.is_element(local_dof_indices[i]))
1337 {
1338 // this is a locally owned DoF, assign new number if not
1339 // assigned a number yet
1340 const unsigned int idx =
1341 locally_owned_dof_indices.index_within_set(
1342 local_dof_indices[i]);
1343 if (new_indices[idx] == numbers::invalid_dof_index)
1344 {
1345 new_indices[idx] =
1346 my_starting_index + current_next_free_dof_offset;
1347 ++current_next_free_dof_offset;
1348 }
1349 }
1350 }
1351 }
1352
1353 return current_next_free_dof_offset;
1354 }
1355 } // namespace
1356
1357
1358
1359 template <int dim, int spacedim>
1360 void
1362 {
1363 std::vector<types::global_dof_index> renumbering(
1365
1366 types::global_dof_index next_free_dof_offset = 0;
1367 const IndexSet locally_owned = dof_handler.locally_owned_dofs();
1368
1369 // in the function we call recursively, we want to number DoFs so
1370 // that global cell zero starts with DoF zero, regardless of how
1371 // DoFs were previously numbered. to this end, we need to figure
1372 // out which DoF index the current processor should start with.
1373 //
1374 // if this is a sequential triangulation, then obviously the starting
1375 // index is zero. otherwise, make sure we get contiguous, successive
1376 // ranges on each processor. note that since the number of locally owned
1377 // DoFs is an invariant under renumbering, we can easily compute this
1378 // starting index by just accumulating over the number of locally owned
1379 // DoFs for all previous processes
1380 types::global_dof_index my_starting_index = 0;
1381
1383 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1384 &dof_handler.get_triangulation()))
1385 {
1386#ifdef DEAL_II_WITH_MPI
1388 dof_handler.locally_owned_dofs().n_elements();
1389 const int ierr = MPI_Exscan(
1391 &my_starting_index,
1392 1,
1393 Utilities::MPI::mpi_type_id_for_type<types::global_dof_index>,
1394 MPI_SUM,
1395 tria->get_mpi_communicator());
1396 AssertThrowMPI(ierr);
1397#endif
1398 }
1399
1402 *>(&dof_handler.get_triangulation()))
1403 {
1404#ifdef DEAL_II_WITH_P4EST
1405 // this is a distributed Triangulation. we need to traverse the coarse
1406 // cells in the order p4est does to match the z-order actually used
1407 // by p4est. this requires using the renumbering of coarse cells
1408 // we do before we hand things off to p4est
1409 for (unsigned int c = 0; c < tria->n_cells(0); ++c)
1410 {
1411 const unsigned int coarse_cell_index =
1412 tria->get_p4est_tree_to_coarse_cell_permutation()[c];
1413
1415 this_cell(tria, 0, coarse_cell_index, &dof_handler);
1416
1417 next_free_dof_offset =
1418 compute_hierarchical_recursive<dim>(next_free_dof_offset,
1419 my_starting_index,
1420 this_cell,
1421 locally_owned,
1422 renumbering);
1423 }
1424#else
1426#endif
1427 }
1428 else
1429 {
1430 // this is not a distributed Triangulation, so we can traverse coarse
1431 // cells in the normal order
1432 for (typename DoFHandler<dim, spacedim>::cell_iterator cell =
1433 dof_handler.begin(0);
1434 cell != dof_handler.end(0);
1435 ++cell)
1436 next_free_dof_offset =
1437 compute_hierarchical_recursive<dim>(next_free_dof_offset,
1438 my_starting_index,
1439 cell,
1440 locally_owned,
1441 renumbering);
1442 }
1443
1444 // verify that the last numbered degree of freedom is either
1445 // equal to the number of degrees of freedom in total (the
1446 // sequential case) or in the distributed case at least
1447 // makes sense
1448 Assert((next_free_dof_offset == dof_handler.n_locally_owned_dofs()) ||
1449 ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
1450 (next_free_dof_offset <= dof_handler.n_dofs())),
1452
1453 // make sure that all local DoFs got new numbers assigned
1454 Assert(std::find(renumbering.begin(),
1455 renumbering.end(),
1456 numbers::invalid_dof_index) == renumbering.end(),
1458
1459 dof_handler.renumber_dofs(renumbering);
1460 }
1461
1462
1463
1464 template <int dim, int spacedim>
1465 void
1467 const std::vector<bool> &selected_dofs)
1468 {
1469 std::vector<types::global_dof_index> renumbering(
1470 dof_handler.n_dofs(), numbers::invalid_dof_index);
1471 compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs);
1472
1473 dof_handler.renumber_dofs(renumbering);
1474 }
1475
1476
1477
1478 template <int dim, int spacedim>
1479 void
1481 const std::vector<bool> &selected_dofs,
1482 const unsigned int level)
1483 {
1486
1487 std::vector<types::global_dof_index> renumbering(
1490 dof_handler,
1491 selected_dofs,
1492 level);
1493
1494 dof_handler.renumber_dofs(level, renumbering);
1495 }
1496
1497
1498
1499 template <int dim, int spacedim>
1500 void
1502 std::vector<types::global_dof_index> &new_indices,
1503 const DoFHandler<dim, spacedim> &dof_handler,
1504 const std::vector<bool> &selected_dofs)
1505 {
1506 const types::global_dof_index n_dofs = dof_handler.n_dofs();
1507 Assert(selected_dofs.size() == n_dofs,
1508 ExcDimensionMismatch(selected_dofs.size(), n_dofs));
1509
1510 // re-sort the dofs according to
1511 // their selection state
1512 Assert(new_indices.size() == n_dofs,
1513 ExcDimensionMismatch(new_indices.size(), n_dofs));
1514
1515 const types::global_dof_index n_selected_dofs =
1516 std::count(selected_dofs.begin(), selected_dofs.end(), false);
1517
1518 types::global_dof_index next_unselected = 0;
1519 types::global_dof_index next_selected = n_selected_dofs;
1520 for (types::global_dof_index i = 0; i < n_dofs; ++i)
1521 if (selected_dofs[i] == false)
1522 {
1523 new_indices[i] = next_unselected;
1524 ++next_unselected;
1525 }
1526 else
1527 {
1528 new_indices[i] = next_selected;
1529 ++next_selected;
1530 }
1531 Assert(next_unselected == n_selected_dofs, ExcInternalError());
1532 Assert(next_selected == n_dofs, ExcInternalError());
1533 }
1534
1535
1536
1537 template <int dim, int spacedim>
1538 void
1540 std::vector<types::global_dof_index> &new_indices,
1541 const DoFHandler<dim, spacedim> &dof_handler,
1542 const std::vector<bool> &selected_dofs,
1543 const unsigned int level)
1544 {
1547
1548 const unsigned int n_dofs = dof_handler.n_dofs(level);
1549 Assert(selected_dofs.size() == n_dofs,
1550 ExcDimensionMismatch(selected_dofs.size(), n_dofs));
1551
1552 // re-sort the dofs according to
1553 // their selection state
1554 Assert(new_indices.size() == n_dofs,
1555 ExcDimensionMismatch(new_indices.size(), n_dofs));
1556
1557 const unsigned int n_selected_dofs =
1558 std::count(selected_dofs.begin(), selected_dofs.end(), false);
1559
1560 unsigned int next_unselected = 0;
1561 unsigned int next_selected = n_selected_dofs;
1562 for (unsigned int i = 0; i < n_dofs; ++i)
1563 if (selected_dofs[i] == false)
1564 {
1565 new_indices[i] = next_unselected;
1566 ++next_unselected;
1567 }
1568 else
1569 {
1570 new_indices[i] = next_selected;
1571 ++next_selected;
1572 }
1573 Assert(next_unselected == n_selected_dofs, ExcInternalError());
1574 Assert(next_selected == n_dofs, ExcInternalError());
1575 }
1576
1577
1578
1579 template <int dim, int spacedim>
1580 void
1583 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
1584 &cells)
1585 {
1586 std::vector<types::global_dof_index> renumbering(
1587 dof.n_locally_owned_dofs());
1588 std::vector<types::global_dof_index> reverse(dof.n_locally_owned_dofs());
1589 compute_cell_wise(renumbering, reverse, dof, cells);
1590
1591 dof.renumber_dofs(renumbering);
1592 }
1593
1594
1595 template <int dim, int spacedim>
1596 void
1598 std::vector<types::global_dof_index> &new_indices,
1599 std::vector<types::global_dof_index> &reverse,
1600 const DoFHandler<dim, spacedim> &dof,
1601 const typename std::vector<
1603 {
1605 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1606 &dof.get_triangulation()))
1607 {
1608 AssertDimension(cells.size(), p->n_locally_owned_active_cells());
1609 }
1610 else
1611 {
1612 AssertDimension(cells.size(), dof.get_triangulation().n_active_cells());
1613 }
1614
1615 const auto n_owned_dofs = dof.n_locally_owned_dofs();
1616
1617 // Actually, we compute the inverse of the reordering vector, called reverse
1618 // here. Later, its inverse is computed into new_indices, which is the
1619 // return argument.
1620
1621 AssertDimension(new_indices.size(), n_owned_dofs);
1622 AssertDimension(reverse.size(), n_owned_dofs);
1623
1624 // For continuous elements, we must make sure, that each dof is reordered
1625 // only once.
1626 std::vector<bool> already_sorted(n_owned_dofs, false);
1627 std::vector<types::global_dof_index> cell_dofs;
1628
1629 const auto &owned_dofs = dof.locally_owned_dofs();
1630
1631 unsigned int index = 0;
1632
1633 for (const auto &cell : cells)
1634 {
1635 // Determine the number of dofs on this cell and reinit the
1636 // vector storing these numbers.
1637 const unsigned int n_cell_dofs = cell->get_fe().n_dofs_per_cell();
1638 cell_dofs.resize(n_cell_dofs);
1639
1640 cell->get_active_or_mg_dof_indices(cell_dofs);
1641
1642 // Sort here to make sure that degrees of freedom inside a single cell
1643 // are in the same order after renumbering.
1644 std::sort(cell_dofs.begin(), cell_dofs.end());
1645
1646 for (const auto dof : cell_dofs)
1647 {
1648 const auto local_dof = owned_dofs.index_within_set(dof);
1649 if (local_dof != numbers::invalid_dof_index &&
1650 !already_sorted[local_dof])
1651 {
1652 already_sorted[local_dof] = true;
1653 reverse[index++] = local_dof;
1654 }
1655 }
1656 }
1657 Assert(index == n_owned_dofs,
1658 ExcMessage(
1659 "Traversing over the given set of cells did not cover all "
1660 "degrees of freedom in the DoFHandler. Does the set of cells "
1661 "not include all active cells?"));
1662
1663 for (types::global_dof_index i = 0; i < reverse.size(); ++i)
1664 new_indices[reverse[i]] = owned_dofs.nth_index_in_set(i);
1665 }
1666
1667
1668
1669 template <int dim, int spacedim>
1670 void
1672 const unsigned int level,
1673 const typename std::vector<
1675 {
1678
1679 std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1680 std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1681
1682 compute_cell_wise(renumbering, reverse, dof, level, cells);
1683 dof.renumber_dofs(level, renumbering);
1684 }
1685
1686
1687
1688 template <int dim, int spacedim>
1689 void
1691 std::vector<types::global_dof_index> &new_order,
1692 std::vector<types::global_dof_index> &reverse,
1693 const DoFHandler<dim, spacedim> &dof,
1694 const unsigned int level,
1695 const typename std::vector<
1697 {
1698 Assert(cells.size() == dof.get_triangulation().n_cells(level),
1699 ExcDimensionMismatch(cells.size(),
1701 Assert(new_order.size() == dof.n_dofs(level),
1702 ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
1703 Assert(reverse.size() == dof.n_dofs(level),
1704 ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
1705
1706 unsigned int n_global_dofs = dof.n_dofs(level);
1707 unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
1708
1709 std::vector<bool> already_sorted(n_global_dofs, false);
1710 std::vector<types::global_dof_index> cell_dofs(n_cell_dofs);
1711
1712 unsigned int global_index = 0;
1713
1714 for (const auto &cell : cells)
1715 {
1716 Assert(cell->level() == static_cast<int>(level), ExcInternalError());
1717
1718 cell->get_active_or_mg_dof_indices(cell_dofs);
1719 std::sort(cell_dofs.begin(), cell_dofs.end());
1720
1721 for (unsigned int i = 0; i < n_cell_dofs; ++i)
1722 {
1723 if (!already_sorted[cell_dofs[i]])
1724 {
1725 already_sorted[cell_dofs[i]] = true;
1726 reverse[global_index++] = cell_dofs[i];
1727 }
1728 }
1729 }
1730 Assert(global_index == n_global_dofs,
1731 ExcMessage(
1732 "Traversing over the given set of cells did not cover all "
1733 "degrees of freedom in the DoFHandler. Does the set of cells "
1734 "not include all cells of the specified level?"));
1735
1736 for (types::global_dof_index i = 0; i < new_order.size(); ++i)
1737 new_order[reverse[i]] = i;
1738 }
1739
1740
1741
1742 template <int dim, int spacedim>
1743 void
1745 const Tensor<1, spacedim> &direction,
1746 const bool dof_wise_renumbering)
1747 {
1748 std::vector<types::global_dof_index> renumbering(dof.n_dofs());
1749 std::vector<types::global_dof_index> reverse(dof.n_dofs());
1751 renumbering, reverse, dof, direction, dof_wise_renumbering);
1752
1753 dof.renumber_dofs(renumbering);
1754 }
1755
1756
1757
1758 template <int dim, int spacedim>
1759 void
1760 compute_downstream(std::vector<types::global_dof_index> &new_indices,
1761 std::vector<types::global_dof_index> &reverse,
1762 const DoFHandler<dim, spacedim> &dof,
1763 const Tensor<1, spacedim> &direction,
1764 const bool dof_wise_renumbering)
1765 {
1767 &dof.get_triangulation()) == nullptr),
1769
1770 if (dof_wise_renumbering == false)
1771 {
1772 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
1773 ordered_cells;
1774 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
1775 const CompareDownstream<
1777 spacedim>
1778 comparator(direction);
1779
1780 for (const auto &cell : dof.active_cell_iterators())
1781 ordered_cells.push_back(cell);
1782
1783 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
1784
1785 compute_cell_wise(new_indices, reverse, dof, ordered_cells);
1786 }
1787 else
1788 {
1789 // similar code as for
1790 // DoFTools::map_dofs_to_support_points, but
1791 // need to do this for general DoFHandler<dim, spacedim> classes and
1792 // want to be able to sort the result
1793 // (otherwise, could use something like
1794 // DoFTools::map_support_points_to_dofs)
1795 const unsigned int n_dofs = dof.n_dofs();
1796 std::vector<std::pair<Point<spacedim>, unsigned int>>
1797 support_point_list(n_dofs);
1798
1799 const hp::FECollection<dim> &fe_collection = dof.get_fe_collection();
1800 Assert(fe_collection[0].has_support_points(),
1802 hp::QCollection<dim> quadrature_collection;
1803 for (unsigned int comp = 0; comp < fe_collection.size(); ++comp)
1804 {
1805 Assert(fe_collection[comp].has_support_points(),
1807 quadrature_collection.push_back(
1808 Quadrature<dim>(fe_collection[comp].get_unit_support_points()));
1809 }
1810 hp::FEValues<dim, spacedim> hp_fe_values(fe_collection,
1811 quadrature_collection,
1813
1814 std::vector<bool> already_touched(n_dofs, false);
1815
1816 std::vector<types::global_dof_index> local_dof_indices;
1817
1818 for (const auto &cell : dof.active_cell_iterators())
1819 {
1820 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1821 local_dof_indices.resize(dofs_per_cell);
1822 hp_fe_values.reinit(cell);
1823 const FEValues<dim> &fe_values =
1824 hp_fe_values.get_present_fe_values();
1825 cell->get_active_or_mg_dof_indices(local_dof_indices);
1826 const std::vector<Point<spacedim>> &points =
1827 fe_values.get_quadrature_points();
1828 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1829 if (!already_touched[local_dof_indices[i]])
1830 {
1831 support_point_list[local_dof_indices[i]].first = points[i];
1832 support_point_list[local_dof_indices[i]].second =
1833 local_dof_indices[i];
1834 already_touched[local_dof_indices[i]] = true;
1835 }
1836 }
1837
1838 ComparePointwiseDownstream<spacedim> comparator(direction);
1839 std::sort(support_point_list.begin(),
1840 support_point_list.end(),
1841 comparator);
1842 for (types::global_dof_index i = 0; i < n_dofs; ++i)
1843 new_indices[support_point_list[i].second] = i;
1844 }
1845 }
1846
1847
1848
1849 template <int dim, int spacedim>
1850 void
1852 const unsigned int level,
1853 const Tensor<1, spacedim> &direction,
1854 const bool dof_wise_renumbering)
1855 {
1856 std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1857 std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1859 renumbering, reverse, dof, level, direction, dof_wise_renumbering);
1860
1861 dof.renumber_dofs(level, renumbering);
1862 }
1863
1864
1865
1866 template <int dim, int spacedim>
1867 void
1868 compute_downstream(std::vector<types::global_dof_index> &new_indices,
1869 std::vector<types::global_dof_index> &reverse,
1870 const DoFHandler<dim, spacedim> &dof,
1871 const unsigned int level,
1872 const Tensor<1, spacedim> &direction,
1873 const bool dof_wise_renumbering)
1874 {
1875 if (dof_wise_renumbering == false)
1876 {
1877 std::vector<typename DoFHandler<dim, spacedim>::level_cell_iterator>
1878 ordered_cells;
1879 ordered_cells.reserve(dof.get_triangulation().n_cells(level));
1880 const CompareDownstream<
1882 spacedim>
1883 comparator(direction);
1884
1886 dof.begin(level);
1888 dof.end(level);
1889
1890 while (p != end)
1891 {
1892 ordered_cells.push_back(p);
1893 ++p;
1894 }
1895 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
1896
1897 compute_cell_wise(new_indices, reverse, dof, level, ordered_cells);
1898 }
1899 else
1900 {
1903 const unsigned int n_dofs = dof.n_dofs(level);
1904 std::vector<std::pair<Point<spacedim>, unsigned int>>
1905 support_point_list(n_dofs);
1906
1907 const Quadrature<dim> q_dummy(dof.get_fe().get_unit_support_points());
1908 FEValues<dim, spacedim> fe_values(dof.get_fe(),
1909 q_dummy,
1911
1912 std::vector<bool> already_touched(dof.n_dofs(), false);
1913
1914 const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
1915 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1917 dof.begin(level);
1919 dof.end(level);
1920 for (; begin != end; ++begin)
1921 {
1923 &begin_tria = begin;
1924 begin->get_active_or_mg_dof_indices(local_dof_indices);
1925 fe_values.reinit(begin_tria);
1926 const std::vector<Point<spacedim>> &points =
1927 fe_values.get_quadrature_points();
1928 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1929 if (!already_touched[local_dof_indices[i]])
1930 {
1931 support_point_list[local_dof_indices[i]].first = points[i];
1932 support_point_list[local_dof_indices[i]].second =
1933 local_dof_indices[i];
1934 already_touched[local_dof_indices[i]] = true;
1935 }
1936 }
1937
1938 ComparePointwiseDownstream<spacedim> comparator(direction);
1939 std::sort(support_point_list.begin(),
1940 support_point_list.end(),
1941 comparator);
1942 for (types::global_dof_index i = 0; i < n_dofs; ++i)
1943 new_indices[support_point_list[i].second] = i;
1944 }
1945 }
1946
1947
1948
1952 namespace internal
1953 {
1954 template <int dim>
1956 {
1965
1970 : center(center)
1971 , counter(counter)
1972 {}
1973
1977 template <class DHCellIterator>
1978 bool
1979 operator()(const DHCellIterator &c1, const DHCellIterator &c2) const
1980 {
1981 // dispatch to
1982 // dimension-dependent functions
1983 return compare(c1, c2, std::integral_constant<int, dim>());
1984 }
1985
1986 private:
1990 template <class DHCellIterator, int xdim>
1991 bool
1992 compare(const DHCellIterator &c1,
1993 const DHCellIterator &c2,
1994 std::integral_constant<int, xdim>) const
1995 {
1996 const Tensor<1, dim> v1 = c1->center() - center;
1997 const Tensor<1, dim> v2 = c2->center() - center;
1998 const double s1 = std::atan2(v1[0], v1[1]);
1999 const double s2 = std::atan2(v2[0], v2[1]);
2000 return (counter ? (s1 > s2) : (s2 > s1));
2001 }
2002
2003
2008 template <class DHCellIterator>
2009 bool
2010 compare(const DHCellIterator &,
2011 const DHCellIterator &,
2012 std::integral_constant<int, 1>) const
2013 {
2014 Assert(dim >= 2,
2015 ExcMessage("This operation only makes sense for dim>=2."));
2016 return false;
2017 }
2018 };
2019 } // namespace internal
2020
2021
2022
2023 template <int dim, int spacedim>
2024 void
2026 const Point<spacedim> &center,
2027 const bool counter)
2028 {
2029 std::vector<types::global_dof_index> renumbering(dof.n_dofs());
2030 compute_clockwise_dg(renumbering, dof, center, counter);
2031
2032 dof.renumber_dofs(renumbering);
2033 }
2034
2035
2036
2037 template <int dim, int spacedim>
2038 void
2039 compute_clockwise_dg(std::vector<types::global_dof_index> &new_indices,
2040 const DoFHandler<dim, spacedim> &dof,
2041 const Point<spacedim> &center,
2042 const bool counter)
2043 {
2044 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2045 ordered_cells;
2046 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
2047 internal::ClockCells<spacedim> comparator(center, counter);
2048
2049 for (const auto &cell : dof.active_cell_iterators())
2050 ordered_cells.push_back(cell);
2051
2052 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
2053
2054 std::vector<types::global_dof_index> reverse(new_indices.size());
2055 compute_cell_wise(new_indices, reverse, dof, ordered_cells);
2056 }
2057
2058
2059
2060 template <int dim, int spacedim>
2061 void
2063 const unsigned int level,
2064 const Point<spacedim> &center,
2065 const bool counter)
2066 {
2067 std::vector<typename DoFHandler<dim, spacedim>::level_cell_iterator>
2068 ordered_cells;
2069 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
2070 internal::ClockCells<spacedim> comparator(center, counter);
2071
2073 dof.begin(level);
2075 dof.end(level);
2076
2077 while (p != end)
2078 {
2079 ordered_cells.push_back(p);
2080 ++p;
2081 }
2082 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
2083
2084 cell_wise(dof, level, ordered_cells);
2085 }
2086
2087
2088
2089 template <int dim, int spacedim>
2090 void
2092 {
2093 std::vector<types::global_dof_index> renumbering(
2094 dof_handler.n_dofs(), numbers::invalid_dof_index);
2095 compute_random(renumbering, dof_handler);
2096
2097 dof_handler.renumber_dofs(renumbering);
2098 }
2099
2100
2101
2102 template <int dim, int spacedim>
2103 void
2104 random(DoFHandler<dim, spacedim> &dof_handler, const unsigned int level)
2105 {
2108
2109 std::vector<types::global_dof_index> renumbering(
2110 dof_handler.locally_owned_mg_dofs(level).n_elements(),
2112
2113 compute_random(renumbering, dof_handler, level);
2114
2115 dof_handler.renumber_dofs(level, renumbering);
2116 }
2117
2118
2119
2120 template <int dim, int spacedim>
2121 void
2122 compute_random(std::vector<types::global_dof_index> &new_indices,
2123 const DoFHandler<dim, spacedim> &dof_handler)
2124 {
2125 const types::global_dof_index n_dofs = dof_handler.n_dofs();
2126 Assert(new_indices.size() == n_dofs,
2127 ExcDimensionMismatch(new_indices.size(), n_dofs));
2128
2129 std::iota(new_indices.begin(),
2130 new_indices.end(),
2132
2133 // shuffle the elements; the following is essentially std::shuffle (which
2134 // is new in C++11) but with a boost URNG
2135 // we could use std::mt19937 here but doing so results in compiler-dependent
2136 // output
2137 ::boost::mt19937 random_number_generator;
2138 for (unsigned int i = 1; i < n_dofs; ++i)
2139 {
2140 // get a random number between 0 and i (inclusive)
2141 const unsigned int j =
2142 ::boost::random::uniform_int_distribution<>(0, i)(
2143 random_number_generator);
2144
2145 // if possible, swap the elements
2146 if (i != j)
2147 std::swap(new_indices[i], new_indices[j]);
2148 }
2149 }
2150
2151
2152
2153 template <int dim, int spacedim>
2154 void
2155 compute_random(std::vector<types::global_dof_index> &new_indices,
2156 const DoFHandler<dim, spacedim> &dof_handler,
2157 const unsigned int level)
2158 {
2159 const types::global_dof_index n_dofs = dof_handler.n_dofs(level);
2160 Assert(new_indices.size() == n_dofs,
2161 ExcDimensionMismatch(new_indices.size(), n_dofs));
2162
2163 std::iota(new_indices.begin(),
2164 new_indices.end(),
2166
2167 // shuffle the elements; the following is essentially std::shuffle (which
2168 // is new in C++11) but with a boost URNG
2169 // we could use std::mt19937 here but doing so results in
2170 // compiler-dependent output
2171 ::boost::mt19937 random_number_generator;
2172 for (unsigned int i = 1; i < n_dofs; ++i)
2173 {
2174 // get a random number between 0 and i (inclusive)
2175 const unsigned int j =
2176 ::boost::random::uniform_int_distribution<>(0, i)(
2177 random_number_generator);
2178
2179 // if possible, swap the elements
2180 if (i != j)
2181 std::swap(new_indices[i], new_indices[j]);
2182 }
2183 }
2184
2185
2186
2187 template <int dim, int spacedim>
2188 void
2190 {
2191 Assert(
2192 (!dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2193 &dof_handler.get_triangulation())),
2194 ExcMessage(
2195 "Parallel triangulations are already enumerated according to their MPI process id."));
2196
2197 std::vector<types::global_dof_index> renumbering(
2198 dof_handler.n_dofs(), numbers::invalid_dof_index);
2199 compute_subdomain_wise(renumbering, dof_handler);
2200
2201 dof_handler.renumber_dofs(renumbering);
2202 }
2203
2204
2205
2206 template <int dim, int spacedim>
2207 void
2208 compute_subdomain_wise(std::vector<types::global_dof_index> &new_dof_indices,
2209 const DoFHandler<dim, spacedim> &dof_handler)
2210 {
2211 const types::global_dof_index n_dofs = dof_handler.n_dofs();
2212 Assert(new_dof_indices.size() == n_dofs,
2213 ExcDimensionMismatch(new_dof_indices.size(), n_dofs));
2214
2215 // first get the association of each dof
2216 // with a subdomain and determine the total
2217 // number of subdomain ids used
2218 std::vector<types::subdomain_id> subdomain_association(n_dofs);
2219 DoFTools::get_subdomain_association(dof_handler, subdomain_association);
2220 const unsigned int n_subdomains =
2221 *std::max_element(subdomain_association.begin(),
2222 subdomain_association.end()) +
2223 1;
2224
2225 // then renumber the subdomains by first
2226 // looking at those belonging to subdomain
2227 // 0, then those of subdomain 1, etc. note
2228 // that the algorithm is stable, i.e. if
2229 // two dofs i,j have i<j and belong to the
2230 // same subdomain, then they will be in
2231 // this order also after reordering
2232 std::fill(new_dof_indices.begin(),
2233 new_dof_indices.end(),
2235 types::global_dof_index next_free_index = 0;
2236 for (types::subdomain_id subdomain = 0; subdomain < n_subdomains;
2237 ++subdomain)
2238 for (types::global_dof_index i = 0; i < n_dofs; ++i)
2239 if (subdomain_association[i] == subdomain)
2240 {
2241 Assert(new_dof_indices[i] == numbers::invalid_dof_index,
2243 new_dof_indices[i] = next_free_index;
2244 ++next_free_index;
2245 }
2246
2247 // we should have numbered all dofs
2248 Assert(next_free_index == n_dofs, ExcInternalError());
2249 Assert(std::find(new_dof_indices.begin(),
2250 new_dof_indices.end(),
2251 numbers::invalid_dof_index) == new_dof_indices.end(),
2253 }
2254
2255
2256
2257 template <int dim, int spacedim>
2258 void
2260 {
2261 std::vector<types::global_dof_index> renumbering(
2263 compute_support_point_wise(renumbering, dof_handler);
2264
2265 // If there is only one component then there is nothing to do, so check
2266 // first:
2267 if (Utilities::MPI::max(renumbering.size(),
2268 dof_handler.get_mpi_communicator()) > 0)
2269 dof_handler.renumber_dofs(renumbering);
2270 }
2271
2272
2273
2274 template <int dim, int spacedim>
2275 void
2277 std::vector<types::global_dof_index> &new_dof_indices,
2278 const DoFHandler<dim, spacedim> &dof_handler)
2279 {
2280 const types::global_dof_index n_dofs = dof_handler.n_locally_owned_dofs();
2281 Assert(new_dof_indices.size() == n_dofs,
2282 ExcDimensionMismatch(new_dof_indices.size(), n_dofs));
2283
2284 // This renumbering occurs in three steps:
2285 // 1. Compute the component-wise renumbering so that all DoFs of component
2286 // i are less than the DoFs of component i + 1.
2287 // 2. Compute a second renumbering component_to_nodal in which the
2288 // renumbering is now, for two components, [u0, v0, u1, v1, ...]: i.e.,
2289 // DoFs are first sorted by component and then by support point.
2290 // 3. Compose the two renumberings to obtain the final result.
2291
2292 // Step 1:
2293 std::vector<types::global_dof_index> component_renumbering(
2295 compute_component_wise<dim, spacedim>(component_renumbering,
2296 dof_handler.begin_active(),
2297 dof_handler.end(),
2298 std::vector<unsigned int>(),
2299 false);
2300
2301 if constexpr (running_in_debug_mode())
2302 {
2303 {
2304 const std::vector<types::global_dof_index> dofs_per_component =
2305 DoFTools::count_dofs_per_fe_component(dof_handler, true);
2306 for (const auto &dpc : dofs_per_component)
2307 Assert(dofs_per_component[0] == dpc, ExcNotImplemented());
2308 }
2309 }
2310 const unsigned int n_components =
2311 dof_handler.get_fe_collection().n_components();
2312 Assert(dof_handler.n_dofs() % n_components == 0, ExcInternalError());
2313 const types::global_dof_index dofs_per_component =
2314 dof_handler.n_dofs() / n_components;
2315 const types::global_dof_index local_dofs_per_component =
2316 dof_handler.n_locally_owned_dofs() / n_components;
2317
2318 // At this point we have no more communication to do - simplify things by
2319 // returning early if possible
2320 if (component_renumbering.empty())
2321 {
2322 new_dof_indices.resize(0);
2323 return;
2324 }
2325 std::fill(new_dof_indices.begin(),
2326 new_dof_indices.end(),
2328 // This index set equals what dof_handler.locally_owned_dofs() would be if
2329 // we executed the componentwise renumbering.
2330 IndexSet component_renumbered_dofs(dof_handler.n_dofs());
2331 // DoFs in each component are now consecutive, which IndexSet::add_indices()
2332 // can exploit by avoiding calls to sort. Make use of that by adding DoFs
2333 // one component at a time:
2334 std::vector<types::global_dof_index> component_dofs(
2335 local_dofs_per_component);
2336 for (unsigned int component = 0; component < n_components; ++component)
2337 {
2338 for (std::size_t i = 0; i < local_dofs_per_component; ++i)
2339 component_dofs[i] =
2340 component_renumbering[n_components * i + component];
2341 component_renumbered_dofs.add_indices(component_dofs.begin(),
2342 component_dofs.end());
2343 }
2344 component_renumbered_dofs.compress();
2345 if constexpr (running_in_debug_mode())
2346 {
2347 {
2348 IndexSet component_renumbered_dofs2(dof_handler.n_dofs());
2349 component_renumbered_dofs2.add_indices(component_renumbering.begin(),
2350 component_renumbering.end());
2351 Assert(component_renumbered_dofs2 == component_renumbered_dofs,
2353 }
2354 }
2355 for (const FiniteElement<dim, spacedim> &fe :
2356 dof_handler.get_fe_collection())
2357 {
2358 AssertThrow(fe.dofs_per_cell == 0 || fe.has_support_points(),
2360 for (unsigned int i = 0; i < fe.n_base_elements(); ++i)
2362 fe.base_element(0).get_unit_support_points() ==
2363 fe.base_element(i).get_unit_support_points(),
2364 ExcMessage(
2365 "All base elements should have the same support points."));
2366 }
2367
2368 std::vector<types::global_dof_index> component_to_nodal(
2370
2371 // Step 2:
2372 std::vector<types::global_dof_index> cell_dofs;
2373 std::vector<types::global_dof_index> component_renumbered_cell_dofs;
2374 const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
2375 // Reuse the original index space for the new renumbering: it is the right
2376 // size and is contiguous on the current processor
2377 auto next_dof_it = locally_owned_dofs.begin();
2378 for (const auto &cell : dof_handler.active_cell_iterators())
2379 if (cell->is_locally_owned())
2380 {
2381 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
2382 cell_dofs.resize(fe.dofs_per_cell);
2383 component_renumbered_cell_dofs.resize(fe.dofs_per_cell);
2384 cell->get_dof_indices(cell_dofs);
2385 // Apply the component renumbering while skipping any ghost dofs. This
2386 // algorithm assumes that all locally owned DoFs before the component
2387 // renumbering are still locally owned afterwards (just with a new
2388 // index).
2389 for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
2390 {
2391 if (locally_owned_dofs.is_element(cell_dofs[i]))
2392 {
2393 const auto local_index =
2394 locally_owned_dofs.index_within_set(cell_dofs[i]);
2395 component_renumbered_cell_dofs[i] =
2396 component_renumbering[local_index];
2397 }
2398 else
2399 {
2400 component_renumbered_cell_dofs[i] =
2402 }
2403 }
2404
2405 for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
2406 {
2407 if (fe.system_to_component_index(i).first == 0 &&
2408 component_renumbered_dofs.is_element(
2409 component_renumbered_cell_dofs[i]))
2410 {
2411 for (unsigned int component = 0;
2412 component < fe.n_components();
2413 ++component)
2414 {
2415 // Since we are indexing in an odd way here it is much
2416 // simpler to compute the permutation separately and
2417 // combine it at the end instead of doing both at once
2418 const auto local_index =
2419 component_renumbered_dofs.index_within_set(
2420 component_renumbered_cell_dofs[i] +
2421 dofs_per_component * component);
2422
2423 if (component_to_nodal[local_index] ==
2425 {
2426 component_to_nodal[local_index] = *next_dof_it;
2427 ++next_dof_it;
2428 }
2429 }
2430 }
2431 }
2432 }
2433
2434 // Step 3:
2435 for (std::size_t i = 0; i < dof_handler.n_locally_owned_dofs(); ++i)
2436 {
2437 const auto local_index =
2438 component_renumbered_dofs.index_within_set(component_renumbering[i]);
2439 new_dof_indices[i] = component_to_nodal[local_index];
2440 }
2441 }
2442
2443
2444
2445 template <int dim,
2446 int spacedim,
2447 typename Number,
2448 typename VectorizedArrayType>
2449 void
2451 DoFHandler<dim, spacedim> &dof_handler,
2453 {
2454 const std::vector<types::global_dof_index> new_global_numbers =
2455 compute_matrix_free_data_locality(dof_handler, matrix_free);
2456 if (matrix_free.get_mg_level() == ::numbers::invalid_unsigned_int)
2457 dof_handler.renumber_dofs(new_global_numbers);
2458 else
2459 dof_handler.renumber_dofs(matrix_free.get_mg_level(), new_global_numbers);
2460 }
2461
2462
2463
2464 template <int dim, int spacedim, typename Number, typename AdditionalDataType>
2465 void
2467 const AffineConstraints<Number> &constraints,
2468 const AdditionalDataType &matrix_free_data)
2469 {
2470 const std::vector<types::global_dof_index> new_global_numbers =
2472 constraints,
2473 matrix_free_data);
2474 if (matrix_free_data.mg_level == ::numbers::invalid_unsigned_int)
2475 dof_handler.renumber_dofs(new_global_numbers);
2476 else
2477 dof_handler.renumber_dofs(matrix_free_data.mg_level, new_global_numbers);
2478 }
2479
2480
2481
2482 template <int dim, int spacedim, typename Number, typename AdditionalDataType>
2483 std::vector<types::global_dof_index>
2485 const DoFHandler<dim, spacedim> &dof_handler,
2486 const AffineConstraints<Number> &constraints,
2487 const AdditionalDataType &matrix_free_data)
2488 {
2489 AdditionalDataType my_mf_data = matrix_free_data;
2490 my_mf_data.initialize_mapping = false;
2491 my_mf_data.tasks_parallel_scheme = AdditionalDataType::none;
2492
2493 typename AdditionalDataType::MatrixFreeType separate_matrix_free;
2494 separate_matrix_free.reinit(::MappingQ1<dim>(),
2495 dof_handler,
2496 constraints,
2497 ::QGauss<1>(2),
2498 my_mf_data);
2499 return compute_matrix_free_data_locality(dof_handler, separate_matrix_free);
2500 }
2501
2502
2503
2504 // Implementation details for matrix-free renumbering
2505 namespace
2506 {
2507 // Compute a vector of lists with number of unknowns of the same category
2508 // in terms of influence from other MPI processes, starting with unknowns
2509 // touched only by the local process and finally a new set of indices
2510 // going to different MPI neighbors. Later passes of the algorithm will
2511 // re-order unknowns within each of these sets.
2512 std::vector<std::vector<unsigned int>>
2513 group_dofs_by_rank_access(
2514 const ::Utilities::MPI::Partitioner &partitioner)
2515 {
2516 // count the number of times a locally owned DoF is accessed by the
2517 // remote ghost data
2518 std::vector<unsigned int> touch_count(partitioner.locally_owned_size());
2519 for (const auto &p : partitioner.import_indices())
2520 for (unsigned int i = p.first; i < p.second; ++i)
2521 touch_count[i]++;
2522
2523 // category 0: DoFs never touched by ghosts
2524 std::vector<std::vector<unsigned int>> result(1);
2525 for (unsigned int i = 0; i < touch_count.size(); ++i)
2526 if (touch_count[i] == 0)
2527 result.back().push_back(i);
2528
2529 // DoFs with 1 appearance can be simply categorized by their (single)
2530 // MPI rank, whereas we need to go an extra round for the remaining DoFs
2531 // by collecting the owning processes by hand
2532 std::map<unsigned int, std::vector<unsigned int>>
2533 multiple_ranks_access_dof;
2534 const std::vector<std::pair<unsigned int, unsigned int>> &import_targets =
2535 partitioner.import_targets();
2536 auto it = partitioner.import_indices().begin();
2537 for (const std::pair<unsigned int, unsigned int> &proc : import_targets)
2538 {
2539 result.emplace_back();
2540 unsigned int count_dofs = 0;
2541 while (count_dofs < proc.second)
2542 {
2543 for (unsigned int i = it->first; i < it->second;
2544 ++i, ++count_dofs)
2545 {
2546 if (touch_count[i] == 1)
2547 result.back().push_back(i);
2548 else
2549 multiple_ranks_access_dof[i].push_back(proc.first);
2550 }
2551 ++it;
2552 }
2553 }
2554 Assert(it == partitioner.import_indices().end(), ExcInternalError());
2555
2556 // Now go from the computed map on DoFs to a map on the processes for
2557 // DoFs with multiple owners, and append the DoFs found this way to our
2558 // global list
2559 std::map<std::vector<unsigned int>,
2560 std::vector<unsigned int>,
2561 std::function<bool(const std::vector<unsigned int> &,
2562 const std::vector<unsigned int> &)>>
2563 dofs_by_rank{[](const std::vector<unsigned int> &a,
2564 const std::vector<unsigned int> &b) {
2565 if (a.size() < b.size())
2566 return true;
2567 if (a.size() == b.size())
2568 {
2569 for (unsigned int i = 0; i < a.size(); ++i)
2570 if (a[i] < b[i])
2571 return true;
2572 else if (a[i] > b[i])
2573 return false;
2574 }
2575 return false;
2576 }};
2577 for (const auto &entry : multiple_ranks_access_dof)
2578 dofs_by_rank[entry.second].push_back(entry.first);
2579
2580 for (const auto &procs : dofs_by_rank)
2581 result.push_back(procs.second);
2582
2583 return result;
2584 }
2585
2586
2587
2588 // Compute two vectors, the first indicating the best numbers for a
2589 // MatrixFree::cell_loop and the second the count of how often a DoF is
2590 // touched by different cell groups, in order to later figure out DoFs
2591 // with far reach and those with only local influence.
2592 template <int dim, typename Number, typename VectorizedArrayType>
2593 std::pair<std::vector<unsigned int>, std::vector<unsigned char>>
2594 compute_mf_numbering(
2596 const unsigned int component)
2597 {
2598 const IndexSet &owned_dofs = matrix_free.get_dof_info(component)
2599 .vector_partitioner->locally_owned_range();
2600 const unsigned int n_comp =
2601 matrix_free.get_dof_handler(component).get_fe().n_components();
2602 Assert(
2603 matrix_free.get_dof_handler(component).get_fe().n_base_elements() == 1,
2605 const bool is_fe_q = dynamic_cast<const FE_Q_Base<dim> *>(
2606 &matrix_free.get_dof_handler(component).get_fe().base_element(0));
2607
2608 const unsigned int fe_degree =
2609 matrix_free.get_dof_handler(component).get_fe().degree;
2610 const unsigned int nn = fe_degree - 1;
2611
2612 // Data structure used further down for the identification of various
2613 // entities in the hierarchical numbering of unknowns. The first number
2614 // indicates the offset from which a given object starts its range of
2615 // numbers in the hierarchical DoF numbering of FE_Q, and the second the
2616 // number of unknowns per component on that particular component. The
2617 // numbers are group by the 3^dim possible objects, listed in
2618 // lexicographic order.
2619 std::array<std::pair<unsigned int, unsigned int>,
2620 ::Utilities::pow(3, dim)>
2621 dofs_on_objects;
2622 if (dim == 1)
2623 {
2624 dofs_on_objects[0] = std::make_pair(0U, 1U);
2625 dofs_on_objects[1] = std::make_pair(2 * n_comp, nn);
2626 dofs_on_objects[2] = std::make_pair(n_comp, 1U);
2627 }
2628 else if (dim == 2)
2629 {
2630 dofs_on_objects[0] = std::make_pair(0U, 1U);
2631 dofs_on_objects[1] = std::make_pair(n_comp * (4 + 2 * nn), nn);
2632 dofs_on_objects[2] = std::make_pair(n_comp, 1U);
2633 dofs_on_objects[3] = std::make_pair(n_comp * 4, nn);
2634 dofs_on_objects[4] = std::make_pair(n_comp * (4 + 4 * nn), nn * nn);
2635 dofs_on_objects[5] = std::make_pair(n_comp * (4 + 1 * nn), nn);
2636 dofs_on_objects[6] = std::make_pair(2 * n_comp, 1U);
2637 dofs_on_objects[7] = std::make_pair(n_comp * (4 + 3 * nn), nn);
2638 dofs_on_objects[8] = std::make_pair(3 * n_comp, 1U);
2639 }
2640 else if (dim == 3)
2641 {
2642 dofs_on_objects[0] = std::make_pair(0U, 1U);
2643 dofs_on_objects[1] = std::make_pair(n_comp * (8 + 2 * nn), nn);
2644 dofs_on_objects[2] = std::make_pair(n_comp, 1U);
2645 dofs_on_objects[3] = std::make_pair(n_comp * 8, nn);
2646 dofs_on_objects[4] =
2647 std::make_pair(n_comp * (8 + 12 * nn + 4 * nn * nn), nn * nn);
2648 dofs_on_objects[5] = std::make_pair(n_comp * (8 + 1 * nn), nn);
2649 dofs_on_objects[6] = std::make_pair(n_comp * 2, 1U);
2650 dofs_on_objects[7] = std::make_pair(n_comp * (8 + 3 * nn), nn);
2651 dofs_on_objects[8] = std::make_pair(n_comp * 3, 1U);
2652 dofs_on_objects[9] = std::make_pair(n_comp * (8 + 8 * nn), nn);
2653 dofs_on_objects[10] =
2654 std::make_pair(n_comp * (8 + 12 * nn + 2 * nn * nn), nn * nn);
2655 dofs_on_objects[11] = std::make_pair(n_comp * (8 + 9 * nn), nn);
2656 dofs_on_objects[12] = std::make_pair(n_comp * (8 + 12 * nn), nn * nn);
2657 dofs_on_objects[13] =
2658 std::make_pair(n_comp * (8 + 12 * nn + 6 * nn * nn), nn * nn * nn);
2659 dofs_on_objects[14] =
2660 std::make_pair(n_comp * (8 + 12 * nn + 1 * nn * nn), nn * nn);
2661 dofs_on_objects[15] = std::make_pair(n_comp * (8 + 10 * nn), nn);
2662 dofs_on_objects[16] =
2663 std::make_pair(n_comp * (8 + 12 * nn + 3 * nn * nn), nn * nn);
2664 dofs_on_objects[17] = std::make_pair(n_comp * (8 + 11 * nn), nn);
2665 dofs_on_objects[18] = std::make_pair(n_comp * 4, 1U);
2666 dofs_on_objects[19] = std::make_pair(n_comp * (8 + 6 * nn), nn);
2667 dofs_on_objects[20] = std::make_pair(n_comp * 5, 1U);
2668 dofs_on_objects[21] = std::make_pair(n_comp * (8 + 4 * nn), nn);
2669 dofs_on_objects[22] =
2670 std::make_pair(n_comp * (8 + 12 * nn + 5 * nn * nn), nn * nn);
2671 dofs_on_objects[23] = std::make_pair(n_comp * (8 + 5 * nn), nn);
2672 dofs_on_objects[24] = std::make_pair(n_comp * 6, 1U);
2673 dofs_on_objects[25] = std::make_pair(n_comp * (8 + 7 * nn), nn);
2674 dofs_on_objects[26] = std::make_pair(n_comp * 7, 1U);
2675 }
2676
2677 const auto renumber_func = [](const types::global_dof_index dof_index,
2678 const IndexSet &owned_dofs,
2679 std::vector<unsigned int> &result,
2680 unsigned int &counter_dof_numbers) {
2681 const types::global_dof_index local_dof_index =
2682 owned_dofs.index_within_set(dof_index);
2683 if (local_dof_index != numbers::invalid_dof_index)
2684 {
2685 AssertIndexRange(local_dof_index, result.size());
2686 if (result[local_dof_index] == numbers::invalid_unsigned_int)
2687 result[local_dof_index] = counter_dof_numbers++;
2688 }
2689 };
2690
2691 unsigned int counter_dof_numbers = 0;
2692 std::vector<unsigned int> dofs_extracted;
2693 std::vector<types::global_dof_index> dof_indices(
2694 matrix_free.get_dof_handler(component).get_fe().dofs_per_cell);
2695
2696 // We now define a lambda function that does two things: (a) determine
2697 // DoF numbers in a way that fit with the order a MatrixFree loop
2698 // travels through the cells (variable 'dof_numbers_mf_order'), and (b)
2699 // determine which unknowns are only touched from within a single range
2700 // of cells and which ones span multiple ranges (variable
2701 // 'touch_count'). Note that this process is done by calling into
2702 // MatrixFree::cell_loop, which gives the right level of granularity
2703 // (when executed in serial) for the scheduled vector operations. Note
2704 // that we pick the unconstrained indices in the hierarchical order for
2705 // part (a) as this makes it easy to identify the DoFs on the individual
2706 // entities, whereas we pick the numbers with constraints eliminated for
2707 // part (b). For the latter, we keep track of each DoF's interaction
2708 // with different ranges of cell batches, i.e., call-backs into the
2709 // operation_on_cell_range() function.
2710 const unsigned int n_owned_dofs = owned_dofs.n_elements();
2711 std::vector<unsigned int> dof_numbers_mf_order(
2712 n_owned_dofs, ::numbers::invalid_unsigned_int);
2713 std::vector<unsigned int> last_touch_by_cell_batch_range(
2714 n_owned_dofs, ::numbers::invalid_unsigned_int);
2715 std::vector<unsigned char> touch_count(n_owned_dofs);
2716
2717 const auto operation_on_cell_range =
2719 unsigned int &,
2720 const unsigned int &,
2721 const std::pair<unsigned int, unsigned int> &cell_range) {
2722 for (unsigned int cell = cell_range.first; cell < cell_range.second;
2723 ++cell)
2724 {
2725 // part (a): assign beneficial numbers
2726 for (unsigned int v = 0;
2727 v < data.n_active_entries_per_cell_batch(cell);
2728 ++v)
2729 {
2730 // get the indices for the dofs in cell_batch
2731 if (data.get_mg_level() == numbers::invalid_unsigned_int)
2732 data.get_cell_iterator(cell, v, component)
2733 ->get_dof_indices(dof_indices);
2734 else
2735 data.get_cell_iterator(cell, v, component)
2736 ->get_mg_dof_indices(dof_indices);
2737
2738 if (is_fe_q)
2739 for (unsigned int a = 0; a < dofs_on_objects.size(); ++a)
2740 {
2741 const auto &r = dofs_on_objects[a];
2742 if (a == 10 || a == 16)
2743 // switch order x-z for y faces in 3d to lexicographic
2744 // layout
2745 for (unsigned int i1 = 0; i1 < nn; ++i1)
2746 for (unsigned int i0 = 0; i0 < nn; ++i0)
2747 for (unsigned int c = 0; c < n_comp; ++c)
2748 renumber_func(
2749 dof_indices[r.first + r.second * c + i1 +
2750 i0 * nn],
2751 owned_dofs,
2752 dof_numbers_mf_order,
2753 counter_dof_numbers);
2754 else
2755 for (unsigned int i = 0; i < r.second; ++i)
2756 for (unsigned int c = 0; c < n_comp; ++c)
2757 renumber_func(
2758 dof_indices[r.first + r.second * c + i],
2759 owned_dofs,
2760 dof_numbers_mf_order,
2761 counter_dof_numbers);
2762 }
2763 else
2764 for (const types::global_dof_index i : dof_indices)
2765 renumber_func(i,
2766 owned_dofs,
2767 dof_numbers_mf_order,
2768 counter_dof_numbers);
2769 }
2770
2771 // part (b): increment the touch count of a dof appearing in the
2772 // current cell batch if it was last touched by another than the
2773 // present cell batch range (we track them via storing the last
2774 // cell batch range that touched a particular dof)
2775 data.get_dof_info(component).get_dof_indices_on_cell_batch(
2776 dofs_extracted, cell);
2777 for (const unsigned int dof_index : dofs_extracted)
2778 if (dof_index < n_owned_dofs &&
2779 last_touch_by_cell_batch_range[dof_index] !=
2780 cell_range.first)
2781 {
2782 ++touch_count[dof_index];
2783 last_touch_by_cell_batch_range[dof_index] =
2784 cell_range.first;
2785 }
2786 }
2787 };
2788
2789 // Finally run the matrix-free loop.
2790 Assert(matrix_free.get_task_info().scheme ==
2792 ExcNotImplemented("Renumbering only available for non-threaded "
2793 "version of MatrixFree::cell_loop"));
2794
2795 matrix_free.template cell_loop<unsigned int, unsigned int>(
2796 operation_on_cell_range, counter_dof_numbers, counter_dof_numbers);
2797
2798 AssertDimension(counter_dof_numbers, n_owned_dofs);
2799
2800 return std::make_pair(dof_numbers_mf_order, touch_count);
2801 }
2802
2803 } // namespace
2804
2805
2806
2807 template <int dim,
2808 int spacedim,
2809 typename Number,
2810 typename VectorizedArrayType>
2811 std::vector<types::global_dof_index>
2813 const DoFHandler<dim, spacedim> &dof_handler,
2815 {
2816 Assert(matrix_free.indices_initialized(),
2817 ExcMessage("You need to set up indices in MatrixFree "
2818 "to be able to compute a renumbering!"));
2819
2820 // try to locate the `DoFHandler` within the given MatrixFree object.
2821 unsigned int component = 0;
2822 for (; component < matrix_free.n_components(); ++component)
2823 if (&matrix_free.get_dof_handler(component) == &dof_handler)
2824 break;
2825
2826 Assert(component < matrix_free.n_components(),
2827 ExcMessage("Could not locate the given DoFHandler in MatrixFree"));
2828
2829 // Summary of the algorithm below:
2830 // (a) compute renumbering of each DoF in the order the corresponding
2831 // object appears in the mf loop -> local_numbering
2832 // (b) determine by how many cell groups (call-back places in the loop) a
2833 // dof is touched -> touch_count
2834 // (c) determine by how many MPI processes a dof is touched
2835 // -> dofs_by_rank_access
2836 // (d) combine both category types of (b) and (c) and list the indices
2837 // according to four categories
2838
2839 const std::vector<std::vector<unsigned int>> dofs_by_rank_access =
2840 group_dofs_by_rank_access(
2841 *matrix_free.get_dof_info(component).vector_partitioner);
2842
2843 const auto &[local_numbering, touch_count] =
2844 compute_mf_numbering(matrix_free, component);
2845
2846 const IndexSet &owned_dofs = matrix_free.get_dof_info(component)
2847 .vector_partitioner->locally_owned_range();
2849 AssertDimension(locally_owned_size, local_numbering.size());
2850 AssertDimension(locally_owned_size, touch_count.size());
2851
2852 // Create a second permutation to group the unknowns into the following
2853 // four categories, to be eventually composed with 'local_renumbering'
2854 enum class Category : unsigned char
2855 {
2856 single, // DoFs only referring to a single batch and MPI process
2857 multiple, // DoFs referring to multiple batches (long-range)
2858 mpi_dof, // DoFs referring to DoFs in exchange with other MPI processes
2859 constrained // DoFs without any reference (constrained DoFs)
2860 };
2861 std::vector<Category> categories(locally_owned_size);
2862 for (unsigned int i = 0; i < locally_owned_size; ++i)
2863 switch (touch_count[i])
2864 {
2865 case 0:
2866 categories[local_numbering[i]] = Category::constrained;
2867 break;
2868 case 1:
2869 categories[local_numbering[i]] = Category::single;
2870 break;
2871 default:
2872 categories[local_numbering[i]] = Category::multiple;
2873 break;
2874 }
2875 for (unsigned int chunk = 1; chunk < dofs_by_rank_access.size(); ++chunk)
2876 for (const auto i : dofs_by_rank_access[chunk])
2877 categories[local_numbering[i]] = Category::mpi_dof;
2878
2879 // Assign numbers to the categories in the order listed in the enum, which
2880 // is done by starting 'counter' for each category at an offset
2881 std::array<unsigned int, 4> n_entries_per_category{};
2882 for (const Category i : categories)
2883 {
2884 AssertIndexRange(static_cast<int>(i), 4);
2885 ++n_entries_per_category[static_cast<int>(i)];
2886 }
2887 std::array<unsigned int, 4> counters{};
2888 for (unsigned int i = 1; i < 4; ++i)
2889 counters[i] = counters[i - 1] + n_entries_per_category[i - 1];
2890 std::vector<unsigned int> numbering_categories;
2891 numbering_categories.reserve(locally_owned_size);
2892 for (const Category category : categories)
2893 {
2894 numbering_categories.push_back(counters[static_cast<int>(category)]);
2895 ++counters[static_cast<int>(category)];
2896 }
2897
2898 // The final numbers are given by the composition of the two permutations
2899 std::vector<::types::global_dof_index> new_global_numbers(
2901 for (unsigned int i = 0; i < locally_owned_size; ++i)
2902 new_global_numbers[i] =
2903 owned_dofs.nth_index_in_set(numbering_categories[local_numbering[i]]);
2904
2905 return new_global_numbers;
2906 }
2907
2908} // namespace DoFRenumbering
2909
2910
2911
2912/*-------------- Explicit Instantiations -------------------------------*/
2913#include "dofs/dof_renumbering.inst"
2914
2915
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
cell_iterator end() const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const Triangulation< dim, spacedim > & get_triangulation() const
const IndexSet & locally_owned_dofs() const
active_cell_iterator begin_active(const unsigned int level=0) const
types::global_dof_index n_dofs() const
typename LevelSelector::cell_iterator level_cell_iterator
cell_iterator begin(const unsigned int level=0) const
MPI_Comm get_mpi_communicator() const
types::global_dof_index n_locally_owned_dofs() const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
size_type row_length(const size_type row) const
size_type column_number(const size_type row, const size_type index) const
const std::vector< Point< spacedim > > & get_quadrature_points() const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
unsigned int n_dofs_per_cell() const
unsigned int n_components() const
const unsigned int dofs_per_cell
Definition fe_data.h:436
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool has_support_points() const
const std::vector< Point< dim > > & get_unit_support_points() const
bool is_primitive() const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
size_type size() const
Definition index_set.h:1787
size_type index_within_set(const size_type global_index) const
Definition index_set.h:2013
size_type n_elements() const
Definition index_set.h:1945
bool is_element(const size_type index) const
Definition index_set.h:1905
ElementIterator begin() const
Definition index_set.h:1721
void subtract_set(const IndexSet &other)
Definition index_set.cc:478
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1994
void compress() const
Definition index_set.h:1795
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1842
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int get_mg_level() const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
bool indices_initialized() const
unsigned int n_components() const
Definition point.h:113
size_type n_rows() const
unsigned int n_active_cells() const
unsigned int n_cells() const
pointer data()
iterator end()
iterator begin()
unsigned int size() const
Definition collection.h:316
void push_back(const FiniteElement< dim, spacedim > &new_fe)
unsigned int n_blocks() const
unsigned int n_components() const
const FEValuesType & get_present_fe_values() const
Definition fe_values.h:695
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, lda > > &cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
Definition fe_values.cc:296
void push_back(const Quadrature< dim_in > &new_quadrature)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > second
Definition grid_out.cc:4633
Point< 2 > first
Definition grid_out.cc:4632
unsigned int level
Definition grid_out.cc:4635
const unsigned int v1
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDoFHandlerNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
Definition mpi.cc:746
std::size_t size
Definition mpi.cc:745
types::global_dof_index locally_owned_size
Definition mpi.cc:833
graph_traits< Graph >::vertices_size_type size_type
graph_traits< Graph >::vertex_descriptor Vertex
adjacency_list< vecS, vecS, undirectedS, property< vertex_color_t, default_color_type, property< vertex_degree_t, int > > > Graph
std::pair< size_type, size_type > Pair
void create_graph(const DoFHandler< dim, spacedim > &dof_handler, const bool use_constraints, boosttypes::Graph &graph, boosttypes::property_map< boosttypes::Graph, boosttypes::vertex_degree_t >::type &graph_degree)
void compute_Cuthill_McKee(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &, const bool reversed_numbering=false, const bool use_constraints=false)
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
void king_ordering(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
void compute_king_ordering(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &, const bool reversed_numbering=false, const bool use_constraints=false)
void compute_minimum_degree(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &, const bool reversed_numbering=false, const bool use_constraints=false)
void minimum_degree(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
void compute_support_point_wise(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &dof_handler)
void compute_subdomain_wise(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &dof_handler)
void matrix_free_data_locality(DoFHandler< dim, spacedim > &dof_handler, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free)
void subdomain_wise(DoFHandler< dim, spacedim > &dof_handler)
types::global_dof_index compute_block_wise(std::vector< types::global_dof_index > &new_dof_indices, const IteratorType &start, const EndIteratorType &end, const bool is_level_operation)
void hierarchical(DoFHandler< dim, spacedim > &dof_handler)
void compute_random(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &dof_handler)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void compute_Cuthill_McKee(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >(), const unsigned int level=numbers::invalid_unsigned_int)
void block_wise(DoFHandler< dim, spacedim > &dof_handler)
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
void compute_sort_selected_dofs_back(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &selected_dofs)
void sort_selected_dofs_back(DoFHandler< dim, spacedim > &dof_handler, const std::vector< bool > &selected_dofs)
void support_point_wise(DoFHandler< dim, spacedim > &dof_handler)
void compute_cell_wise(std::vector< types::global_dof_index > &renumbering, std::vector< types::global_dof_index > &inverse_renumbering, const DoFHandler< dim, spacedim > &dof_handler, const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &cell_order)
std::vector< types::global_dof_index > compute_matrix_free_data_locality(const DoFHandler< dim, spacedim > &dof_handler, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free)
void clockwise_dg(DoFHandler< dim, spacedim > &dof_handler, const Point< spacedim > &center, const bool counter=false)
void random(DoFHandler< dim, spacedim > &dof_handler)
void compute_downstream(std::vector< types::global_dof_index > &new_dof_indices, std::vector< types::global_dof_index > &reverse, const DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering)
void cell_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< typename DoFHandler< dim, spacedim >::active_cell_iterator > &cell_order)
types::global_dof_index compute_component_wise(std::vector< types::global_dof_index > &new_dof_indices, const CellIterator &start, const std_cxx20::type_identity_t< CellIterator > &end, const std::vector< unsigned int > &target_component, const bool is_level_operation)
void compute_clockwise_dg(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandler< dim, spacedim > &dof_handler, const Point< spacedim > &center, const bool counter)
void get_subdomain_association(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::subdomain_id > &subdomain)
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
IndexSet extract_locally_active_dofs(const DoFHandler< dim, spacedim > &dof_handler)
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
IndexSet extract_locally_active_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity, const unsigned int level, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true)
Definition mg_tools.cc:577
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void reorder_Cuthill_McKee(const DynamicSparsityPattern &sparsity, std::vector< DynamicSparsityPattern::size_type > &new_indices, const std::vector< DynamicSparsityPattern::size_type > &starting_indices=std::vector< DynamicSparsityPattern::size_type >())
T sum(const T &t, const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
std::vector< Integer > reverse_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1685
constexpr types::global_dof_index invalid_dof_index
Definition types.h:269
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
typename type_identity< T >::type type_identity_t
Definition type_traits.h:95
unsigned short int fe_index
Definition types.h:72
bool compare(const DHCellIterator &c1, const DHCellIterator &c2, std::integral_constant< int, xdim >) const
bool compare(const DHCellIterator &, const DHCellIterator &, std::integral_constant< int, 1 >) const
ClockCells(const Point< dim > &center, bool counter)
bool operator()(const DHCellIterator &c1, const DHCellIterator &c2) const
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:593