Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
dof_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_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 Assert(result == 0 || result == dof_handler.n_dofs(level),
728
729 if (Utilities::MPI::max(renumbering.size(),
730 dof_handler.get_communicator()) > 0)
731 dof_handler.renumber_dofs(level, renumbering);
732 }
733
734
735
736 template <int dim, int spacedim, typename CellIterator>
738 compute_component_wise(std::vector<types::global_dof_index> &new_indices,
739 const CellIterator &start,
741 const std::vector<unsigned int> &component_order_arg,
742 const bool is_level_operation)
743 {
744 const hp::FECollection<dim, spacedim> &fe_collection =
745 start->get_dof_handler().get_fe_collection();
746
747 // do nothing if the FE has only
748 // one component
749 if (fe_collection.n_components() == 1)
750 {
751 new_indices.resize(0);
752 return 0;
753 }
754
755 // Get a reference to the set of dofs. Note that we assume that all cells
756 // are assumed to be on the same level, otherwise the operation doesn't make
757 // much sense (we will assert this below).
758 const IndexSet &locally_owned_dofs =
759 is_level_operation ?
760 start->get_dof_handler().locally_owned_mg_dofs(start->level()) :
761 start->get_dof_handler().locally_owned_dofs();
762
763 // Copy last argument into a
764 // writable vector.
765 std::vector<unsigned int> component_order(component_order_arg);
766 // If the last argument was an
767 // empty vector, set up things to
768 // store components in the order
769 // found in the system.
770 if (component_order.empty())
771 for (unsigned int i = 0; i < fe_collection.n_components(); ++i)
772 component_order.push_back(i);
773
774 Assert(component_order.size() == fe_collection.n_components(),
775 ExcDimensionMismatch(component_order.size(),
776 fe_collection.n_components()));
777
778 for (const unsigned int component : component_order)
779 {
780 (void)component;
781 AssertIndexRange(component, fe_collection.n_components());
782 }
783
784 // vector to hold the dof indices on
785 // the cell we visit at a time
786 std::vector<types::global_dof_index> local_dof_indices;
787
788 // prebuilt list to which component
789 // a given dof on a cell
790 // should go. note that we get into
791 // trouble here if the shape
792 // function is not primitive, since
793 // then there is no single vector
794 // component to which it
795 // belongs. in this case, assign it
796 // to the first vector component to
797 // which it belongs
798 std::vector<std::vector<unsigned int>> component_list(fe_collection.size());
799 for (unsigned int f = 0; f < fe_collection.size(); ++f)
800 {
801 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
802 const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
803 component_list[f].resize(dofs_per_cell);
804 for (unsigned int i = 0; i < dofs_per_cell; ++i)
805 if (fe.is_primitive(i))
806 component_list[f][i] =
807 component_order[fe.system_to_component_index(i).first];
808 else
809 {
810 const unsigned int comp =
812
813 // then associate this degree
814 // of freedom with this
815 // component
816 component_list[f][i] = component_order[comp];
817 }
818 }
819
820 // set up a map where for each
821 // component the respective degrees
822 // of freedom are collected.
823 //
824 // note that this map is sorted by
825 // component but that within each
826 // component it is NOT sorted by
827 // dof index. note also that some
828 // dof indices are entered
829 // multiply, so we will have to
830 // take care of that
831 std::vector<std::vector<types::global_dof_index>> component_to_dof_map(
832 fe_collection.n_components());
833 for (CellIterator cell = start; cell != end; ++cell)
834 {
835 if (is_level_operation)
836 {
837 // we are dealing with mg dofs, skip foreign level cells:
838 if (!cell->is_locally_owned_on_level())
839 continue;
840 }
841 else
842 {
843 // we are dealing with active dofs, skip the loop if not locally
844 // owned:
845 if (!cell->is_locally_owned())
846 continue;
847 }
848
849 if (is_level_operation)
850 Assert(
851 cell->level() == start->level(),
853 "Multigrid renumbering in compute_component_wise() needs to be applied to a single level!"));
854
855 // on each cell: get dof indices
856 // and insert them into the global
857 // list using their component
858 const types::fe_index fe_index = cell->active_fe_index();
859 const unsigned int dofs_per_cell =
860 fe_collection[fe_index].n_dofs_per_cell();
861 local_dof_indices.resize(dofs_per_cell);
862 cell->get_active_or_mg_dof_indices(local_dof_indices);
863
864 for (unsigned int i = 0; i < dofs_per_cell; ++i)
865 if (locally_owned_dofs.is_element(local_dof_indices[i]))
866 component_to_dof_map[component_list[fe_index][i]].push_back(
867 local_dof_indices[i]);
868 }
869
870 // now we've got all indices sorted
871 // into buckets labeled by their
872 // target component number. we've
873 // only got to traverse this list
874 // and assign the new indices
875 //
876 // however, we first want to sort
877 // the indices entered into the
878 // buckets to preserve the order
879 // within each component and during
880 // this also remove duplicate
881 // entries
882 //
883 // note that we no longer have to
884 // care about non-primitive shape
885 // functions since the buckets
886 // corresponding to the second and
887 // following vector components of a
888 // non-primitive FE will simply be
889 // empty, everything being shoved
890 // into the first one. The same
891 // holds if several components were
892 // joined into a single target.
893 for (unsigned int component = 0; component < fe_collection.n_components();
894 ++component)
895 {
896 std::sort(component_to_dof_map[component].begin(),
897 component_to_dof_map[component].end());
898 component_to_dof_map[component].erase(
899 std::unique(component_to_dof_map[component].begin(),
900 component_to_dof_map[component].end()),
901 component_to_dof_map[component].end());
902 }
903
904 // calculate the number of locally owned
905 // DoFs per bucket
906 const unsigned int n_buckets = fe_collection.n_components();
907 std::vector<types::global_dof_index> shifts(n_buckets);
908
910 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
911 &start->get_dof_handler().get_triangulation())))
912 {
913#ifdef DEAL_II_WITH_MPI
914 std::vector<types::global_dof_index> local_dof_count(n_buckets);
915
916 for (unsigned int c = 0; c < n_buckets; ++c)
917 local_dof_count[c] = component_to_dof_map[c].size();
918
919 std::vector<types::global_dof_index> prefix_dof_count(n_buckets);
920 const int ierr = MPI_Exscan(local_dof_count.data(),
921 prefix_dof_count.data(),
922 n_buckets,
924 MPI_SUM,
925 tria->get_communicator());
926 AssertThrowMPI(ierr);
927
928 std::vector<types::global_dof_index> global_dof_count(n_buckets);
929 Utilities::MPI::sum(local_dof_count,
930 tria->get_communicator(),
931 global_dof_count);
932
933 // calculate shifts
934 types::global_dof_index cumulated = 0;
935 for (unsigned int c = 0; c < n_buckets; ++c)
936 {
937 shifts[c] = prefix_dof_count[c] + cumulated;
938 cumulated += global_dof_count[c];
939 }
940#else
941 (void)tria;
943#endif
944 }
945 else
946 {
947 shifts[0] = 0;
948 for (unsigned int c = 1; c < fe_collection.n_components(); ++c)
949 shifts[c] = shifts[c - 1] + component_to_dof_map[c - 1].size();
950 }
951
952
953
954 // now concatenate all the
955 // components in the order the user
956 // desired to see
957 types::global_dof_index next_free_index = 0;
958 for (unsigned int component = 0; component < fe_collection.n_components();
959 ++component)
960 {
961 next_free_index = shifts[component];
962
963 for (const types::global_dof_index dof_index :
964 component_to_dof_map[component])
965 {
966 Assert(locally_owned_dofs.index_within_set(dof_index) <
967 new_indices.size(),
969 new_indices[locally_owned_dofs.index_within_set(dof_index)] =
970 next_free_index;
971
972 ++next_free_index;
973 }
974 }
975
976 return next_free_index;
977 }
978
979
980
981 template <int dim, int spacedim>
982 void
984 {
985 std::vector<types::global_dof_index> renumbering(
987
989 dim,
990 spacedim,
993 renumbering, dof_handler.begin_active(), dof_handler.end(), false);
994 if (result == 0)
995 return;
996
997 // verify that the last numbered
998 // degree of freedom is either
999 // equal to the number of degrees
1000 // of freedom in total (the
1001 // sequential case) or in the
1002 // distributed case at least
1003 // makes sense
1004 Assert((result == dof_handler.n_locally_owned_dofs()) ||
1005 ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
1006 (result <= dof_handler.n_dofs())),
1008
1009 dof_handler.renumber_dofs(renumbering);
1010 }
1011
1012
1013
1014 template <int dim, int spacedim>
1015 void
1016 block_wise(DoFHandler<dim, spacedim> &dof_handler, const unsigned int level)
1017 {
1020
1021 std::vector<types::global_dof_index> renumbering(
1022 dof_handler.locally_owned_mg_dofs(level).n_elements(),
1024
1026 dof_handler.begin(level);
1028 dof_handler.end(level);
1029
1031 dim,
1032 spacedim,
1035 start,
1036 end,
1037 true);
1038 (void)result;
1039
1040 Assert(result == 0 || result == dof_handler.n_dofs(level),
1042
1043 if (Utilities::MPI::max(renumbering.size(),
1044 dof_handler.get_communicator()) > 0)
1045 dof_handler.renumber_dofs(level, renumbering);
1046 }
1047
1048
1049
1050 template <int dim, int spacedim, class IteratorType, class EndIteratorType>
1052 compute_block_wise(std::vector<types::global_dof_index> &new_indices,
1053 const IteratorType &start,
1054 const EndIteratorType &end,
1055 const bool is_level_operation)
1056 {
1057 const hp::FECollection<dim, spacedim> &fe_collection =
1058 start->get_dof_handler().get_fe_collection();
1059
1060 // do nothing if the FE has only
1061 // one component
1062 if (fe_collection.n_blocks() == 1)
1063 {
1064 new_indices.resize(0);
1065 return 0;
1066 }
1067
1068 // Get a reference to the set of dofs. Note that we assume that all cells
1069 // are assumed to be on the same level, otherwise the operation doesn't make
1070 // much sense (we will assert this below).
1071 const IndexSet &locally_owned_dofs =
1072 is_level_operation ?
1073 start->get_dof_handler().locally_owned_mg_dofs(start->level()) :
1074 start->get_dof_handler().locally_owned_dofs();
1075
1076 // vector to hold the dof indices on
1077 // the cell we visit at a time
1078 std::vector<types::global_dof_index> local_dof_indices;
1079
1080 // prebuilt list to which block
1081 // a given dof on a cell
1082 // should go.
1083 std::vector<std::vector<types::global_dof_index>> block_list(
1084 fe_collection.size());
1085 for (unsigned int f = 0; f < fe_collection.size(); ++f)
1086 {
1087 const FiniteElement<dim, spacedim> &fe = fe_collection[f];
1088 block_list[f].resize(fe.n_dofs_per_cell());
1089 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
1090 block_list[f][i] = fe.system_to_block_index(i).first;
1091 }
1092
1093 // set up a map where for each
1094 // block the respective degrees
1095 // of freedom are collected.
1096 //
1097 // note that this map is sorted by
1098 // block but that within each
1099 // block it is NOT sorted by
1100 // dof index. note also that some
1101 // dof indices are entered
1102 // multiply, so we will have to
1103 // take care of that
1104 std::vector<std::vector<types::global_dof_index>> block_to_dof_map(
1105 fe_collection.n_blocks());
1106 for (IteratorType cell = start; cell != end; ++cell)
1107 {
1108 if (is_level_operation)
1109 {
1110 // we are dealing with mg dofs, skip foreign level cells:
1111 if (!cell->is_locally_owned_on_level())
1112 continue;
1113 }
1114 else
1115 {
1116 // we are dealing with active dofs, skip the loop if not locally
1117 // owned:
1118 if (!cell->is_locally_owned())
1119 continue;
1120 }
1121
1122 if (is_level_operation)
1123 Assert(
1124 cell->level() == start->level(),
1125 ExcMessage(
1126 "Multigrid renumbering in compute_block_wise() needs to be applied to a single level!"));
1127
1128 // on each cell: get dof indices
1129 // and insert them into the global
1130 // list using their component
1131 const types::fe_index fe_index = cell->active_fe_index();
1132 const unsigned int dofs_per_cell =
1133 fe_collection[fe_index].n_dofs_per_cell();
1134 local_dof_indices.resize(dofs_per_cell);
1135 cell->get_active_or_mg_dof_indices(local_dof_indices);
1136
1137 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1138 if (locally_owned_dofs.is_element(local_dof_indices[i]))
1139 block_to_dof_map[block_list[fe_index][i]].push_back(
1140 local_dof_indices[i]);
1141 }
1142
1143 // now we've got all indices sorted
1144 // into buckets labeled by their
1145 // target block number. we've
1146 // only got to traverse this list
1147 // and assign the new indices
1148 //
1149 // however, we first want to sort
1150 // the indices entered into the
1151 // buckets to preserve the order
1152 // within each component and during
1153 // this also remove duplicate
1154 // entries
1155 for (unsigned int block = 0; block < fe_collection.n_blocks(); ++block)
1156 {
1157 std::sort(block_to_dof_map[block].begin(),
1158 block_to_dof_map[block].end());
1159 block_to_dof_map[block].erase(
1160 std::unique(block_to_dof_map[block].begin(),
1161 block_to_dof_map[block].end()),
1162 block_to_dof_map[block].end());
1163 }
1164
1165 // calculate the number of locally owned
1166 // DoFs per bucket
1167 const unsigned int n_buckets = fe_collection.n_blocks();
1168 std::vector<types::global_dof_index> shifts(n_buckets);
1169
1171 (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1172 &start->get_dof_handler().get_triangulation())))
1173 {
1174#ifdef DEAL_II_WITH_MPI
1175 std::vector<types::global_dof_index> local_dof_count(n_buckets);
1176
1177 for (unsigned int c = 0; c < n_buckets; ++c)
1178 local_dof_count[c] = block_to_dof_map[c].size();
1179
1180 std::vector<types::global_dof_index> prefix_dof_count(n_buckets);
1181 const int ierr = MPI_Exscan(local_dof_count.data(),
1182 prefix_dof_count.data(),
1183 n_buckets,
1185 MPI_SUM,
1186 tria->get_communicator());
1187 AssertThrowMPI(ierr);
1188
1189 std::vector<types::global_dof_index> global_dof_count(n_buckets);
1190 Utilities::MPI::sum(local_dof_count,
1191 tria->get_communicator(),
1192 global_dof_count);
1193
1194 // calculate shifts
1195 types::global_dof_index cumulated = 0;
1196 for (unsigned int c = 0; c < n_buckets; ++c)
1197 {
1198 shifts[c] = prefix_dof_count[c] + cumulated;
1199 cumulated += global_dof_count[c];
1200 }
1201#else
1202 (void)tria;
1204#endif
1205 }
1206 else
1207 {
1208 shifts[0] = 0;
1209 for (unsigned int c = 1; c < fe_collection.n_blocks(); ++c)
1210 shifts[c] = shifts[c - 1] + block_to_dof_map[c - 1].size();
1211 }
1212
1213
1214
1215 // now concatenate all the
1216 // components in the order the user
1217 // desired to see
1218 types::global_dof_index next_free_index = 0;
1219 for (unsigned int block = 0; block < fe_collection.n_blocks(); ++block)
1220 {
1221 const typename std::vector<types::global_dof_index>::const_iterator
1222 begin_of_component = block_to_dof_map[block].begin(),
1223 end_of_component = block_to_dof_map[block].end();
1224
1225 next_free_index = shifts[block];
1226
1227 for (typename std::vector<types::global_dof_index>::const_iterator
1228 dof_index = begin_of_component;
1229 dof_index != end_of_component;
1230 ++dof_index)
1231 {
1232 Assert(locally_owned_dofs.index_within_set(*dof_index) <
1233 new_indices.size(),
1235 new_indices[locally_owned_dofs.index_within_set(*dof_index)] =
1236 next_free_index++;
1237 }
1238 }
1239
1240 return next_free_index;
1241 }
1242
1243
1244
1245 namespace
1246 {
1247 // Helper function for DoFRenumbering::hierarchical(). This function
1248 // recurses into the given cell or, if that should be an active (terminal)
1249 // cell, renumbers DoF indices on it. The function starts renumbering with
1250 // 'next_free_dof_index' and returns the first still unused DoF index at the
1251 // end of its operation.
1252 template <int dim, typename CellIteratorType>
1254 compute_hierarchical_recursive(
1255 const types::global_dof_index next_free_dof_offset,
1256 const types::global_dof_index my_starting_index,
1257 const CellIteratorType &cell,
1258 const IndexSet &locally_owned_dof_indices,
1259 std::vector<types::global_dof_index> &new_indices)
1260 {
1261 types::global_dof_index current_next_free_dof_offset =
1262 next_free_dof_offset;
1263
1264 if (cell->has_children())
1265 {
1266 // recursion
1267 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1268 ++c)
1269 current_next_free_dof_offset =
1270 compute_hierarchical_recursive<dim>(current_next_free_dof_offset,
1271 my_starting_index,
1272 cell->child(c),
1273 locally_owned_dof_indices,
1274 new_indices);
1275 }
1276 else
1277 {
1278 // this is a terminal cell. we need to renumber its DoF indices. there
1279 // are now three cases to decide:
1280 // - this is a sequential triangulation: we can just go ahead and
1281 // number
1282 // the DoFs in the order in which we encounter cells. in this case,
1283 // all cells are actually locally owned
1284 // - if this is a parallel::distributed::Triangulation, then we only
1285 // need to work on the locally owned cells since they contain
1286 // all locally owned DoFs.
1287 // - if this is a parallel::shared::Triangulation, then the same
1288 // applies
1289 //
1290 // in all cases, each processor starts new indices so that we get
1291 // a consecutive numbering on each processor, and disjoint ownership
1292 // of the global range of DoF indices
1293 if (cell->is_locally_owned())
1294 {
1295 // first get the existing DoF indices
1296 const unsigned int dofs_per_cell =
1297 cell->get_fe().n_dofs_per_cell();
1298 std::vector<types::global_dof_index> local_dof_indices(
1299 dofs_per_cell);
1300 cell->get_dof_indices(local_dof_indices);
1301
1302 // then loop over the existing DoF indices on this cell
1303 // and see whether it has already been re-numbered (it
1304 // may have been on a face or vertex to a neighboring
1305 // cell that we have encountered before). if not,
1306 // give it a new number and store that number in the
1307 // output array (new_indices)
1308 //
1309 // if this is a parallel triangulation and a DoF index is
1310 // not locally owned, then don't touch it. since
1311 // we don't actually *change* DoF indices (just record new
1312 // numbers in an array), we don't need to worry about
1313 // the decision whether a DoF is locally owned or not changing
1314 // as we progress in renumbering DoFs -- all adjacent cells
1315 // will always agree that a DoF is locally owned or not.
1316 // that said, the first cell to encounter a locally owned DoF
1317 // gets to number it, so the order in which we traverse cells
1318 // matters
1319 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1320 if (locally_owned_dof_indices.is_element(local_dof_indices[i]))
1321 {
1322 // this is a locally owned DoF, assign new number if not
1323 // assigned a number yet
1324 const unsigned int idx =
1325 locally_owned_dof_indices.index_within_set(
1326 local_dof_indices[i]);
1327 if (new_indices[idx] == numbers::invalid_dof_index)
1328 {
1329 new_indices[idx] =
1330 my_starting_index + current_next_free_dof_offset;
1331 ++current_next_free_dof_offset;
1332 }
1333 }
1334 }
1335 }
1336
1337 return current_next_free_dof_offset;
1338 }
1339 } // namespace
1340
1341
1342
1343 template <int dim, int spacedim>
1344 void
1346 {
1347 std::vector<types::global_dof_index> renumbering(
1349
1350 types::global_dof_index next_free_dof_offset = 0;
1351 const IndexSet locally_owned = dof_handler.locally_owned_dofs();
1352
1353 // in the function we call recursively, we want to number DoFs so
1354 // that global cell zero starts with DoF zero, regardless of how
1355 // DoFs were previously numbered. to this end, we need to figure
1356 // out which DoF index the current processor should start with.
1357 //
1358 // if this is a sequential triangulation, then obviously the starting
1359 // index is zero. otherwise, make sure we get contiguous, successive
1360 // ranges on each processor. note that since the number of locally owned
1361 // DoFs is an invariant under renumbering, we can easily compute this
1362 // starting index by just accumulating over the number of locally owned
1363 // DoFs for all previous processes
1364 types::global_dof_index my_starting_index = 0;
1365
1367 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1368 &dof_handler.get_triangulation()))
1369 {
1370#ifdef DEAL_II_WITH_MPI
1371 types::global_dof_index locally_owned_size =
1372 dof_handler.locally_owned_dofs().n_elements();
1373 const int ierr = MPI_Exscan(&locally_owned_size,
1374 &my_starting_index,
1375 1,
1377 MPI_SUM,
1378 tria->get_communicator());
1379 AssertThrowMPI(ierr);
1380#endif
1381 }
1382
1385 *>(&dof_handler.get_triangulation()))
1386 {
1387#ifdef DEAL_II_WITH_P4EST
1388 // this is a distributed Triangulation. we need to traverse the coarse
1389 // cells in the order p4est does to match the z-order actually used
1390 // by p4est. this requires using the renumbering of coarse cells
1391 // we do before we hand things off to p4est
1392 for (unsigned int c = 0; c < tria->n_cells(0); ++c)
1393 {
1394 const unsigned int coarse_cell_index =
1395 tria->get_p4est_tree_to_coarse_cell_permutation()[c];
1396
1398 this_cell(tria, 0, coarse_cell_index, &dof_handler);
1399
1400 next_free_dof_offset =
1401 compute_hierarchical_recursive<dim>(next_free_dof_offset,
1402 my_starting_index,
1403 this_cell,
1404 locally_owned,
1405 renumbering);
1406 }
1407#else
1409#endif
1410 }
1411 else
1412 {
1413 // this is not a distributed Triangulation, so we can traverse coarse
1414 // cells in the normal order
1415 for (typename DoFHandler<dim, spacedim>::cell_iterator cell =
1416 dof_handler.begin(0);
1417 cell != dof_handler.end(0);
1418 ++cell)
1419 next_free_dof_offset =
1420 compute_hierarchical_recursive<dim>(next_free_dof_offset,
1421 my_starting_index,
1422 cell,
1423 locally_owned,
1424 renumbering);
1425 }
1426
1427 // verify that the last numbered degree of freedom is either
1428 // equal to the number of degrees of freedom in total (the
1429 // sequential case) or in the distributed case at least
1430 // makes sense
1431 Assert((next_free_dof_offset == dof_handler.n_locally_owned_dofs()) ||
1432 ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
1433 (next_free_dof_offset <= dof_handler.n_dofs())),
1435
1436 // make sure that all local DoFs got new numbers assigned
1437 Assert(std::find(renumbering.begin(),
1438 renumbering.end(),
1439 numbers::invalid_dof_index) == renumbering.end(),
1441
1442 dof_handler.renumber_dofs(renumbering);
1443 }
1444
1445
1446
1447 template <int dim, int spacedim>
1448 void
1450 const std::vector<bool> &selected_dofs)
1451 {
1452 std::vector<types::global_dof_index> renumbering(
1453 dof_handler.n_dofs(), numbers::invalid_dof_index);
1454 compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs);
1455
1456 dof_handler.renumber_dofs(renumbering);
1457 }
1458
1459
1460
1461 template <int dim, int spacedim>
1462 void
1464 const std::vector<bool> &selected_dofs,
1465 const unsigned int level)
1466 {
1469
1470 std::vector<types::global_dof_index> renumbering(
1473 dof_handler,
1474 selected_dofs,
1475 level);
1476
1477 dof_handler.renumber_dofs(level, renumbering);
1478 }
1479
1480
1481
1482 template <int dim, int spacedim>
1483 void
1485 std::vector<types::global_dof_index> &new_indices,
1486 const DoFHandler<dim, spacedim> &dof_handler,
1487 const std::vector<bool> &selected_dofs)
1488 {
1489 const types::global_dof_index n_dofs = dof_handler.n_dofs();
1490 Assert(selected_dofs.size() == n_dofs,
1491 ExcDimensionMismatch(selected_dofs.size(), n_dofs));
1492
1493 // re-sort the dofs according to
1494 // their selection state
1495 Assert(new_indices.size() == n_dofs,
1496 ExcDimensionMismatch(new_indices.size(), n_dofs));
1497
1498 const types::global_dof_index n_selected_dofs =
1499 std::count(selected_dofs.begin(), selected_dofs.end(), false);
1500
1501 types::global_dof_index next_unselected = 0;
1502 types::global_dof_index next_selected = n_selected_dofs;
1503 for (types::global_dof_index i = 0; i < n_dofs; ++i)
1504 if (selected_dofs[i] == false)
1505 {
1506 new_indices[i] = next_unselected;
1507 ++next_unselected;
1508 }
1509 else
1510 {
1511 new_indices[i] = next_selected;
1512 ++next_selected;
1513 }
1514 Assert(next_unselected == n_selected_dofs, ExcInternalError());
1515 Assert(next_selected == n_dofs, ExcInternalError());
1516 }
1517
1518
1519
1520 template <int dim, int spacedim>
1521 void
1523 std::vector<types::global_dof_index> &new_indices,
1524 const DoFHandler<dim, spacedim> &dof_handler,
1525 const std::vector<bool> &selected_dofs,
1526 const unsigned int level)
1527 {
1530
1531 const unsigned int n_dofs = dof_handler.n_dofs(level);
1532 Assert(selected_dofs.size() == n_dofs,
1533 ExcDimensionMismatch(selected_dofs.size(), n_dofs));
1534
1535 // re-sort the dofs according to
1536 // their selection state
1537 Assert(new_indices.size() == n_dofs,
1538 ExcDimensionMismatch(new_indices.size(), n_dofs));
1539
1540 const unsigned int n_selected_dofs =
1541 std::count(selected_dofs.begin(), selected_dofs.end(), false);
1542
1543 unsigned int next_unselected = 0;
1544 unsigned int next_selected = n_selected_dofs;
1545 for (unsigned int i = 0; i < n_dofs; ++i)
1546 if (selected_dofs[i] == false)
1547 {
1548 new_indices[i] = next_unselected;
1549 ++next_unselected;
1550 }
1551 else
1552 {
1553 new_indices[i] = next_selected;
1554 ++next_selected;
1555 }
1556 Assert(next_unselected == n_selected_dofs, ExcInternalError());
1557 Assert(next_selected == n_dofs, ExcInternalError());
1558 }
1559
1560
1561
1562 template <int dim, int spacedim>
1563 void
1566 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
1567 &cells)
1568 {
1569 std::vector<types::global_dof_index> renumbering(
1570 dof.n_locally_owned_dofs());
1571 std::vector<types::global_dof_index> reverse(dof.n_locally_owned_dofs());
1572 compute_cell_wise(renumbering, reverse, dof, cells);
1573
1574 dof.renumber_dofs(renumbering);
1575 }
1576
1577
1578 template <int dim, int spacedim>
1579 void
1581 std::vector<types::global_dof_index> &new_indices,
1582 std::vector<types::global_dof_index> &reverse,
1583 const DoFHandler<dim, spacedim> &dof,
1584 const typename std::vector<
1586 {
1588 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1589 &dof.get_triangulation()))
1590 {
1591 AssertDimension(cells.size(), p->n_locally_owned_active_cells());
1592 }
1593 else
1594 {
1595 AssertDimension(cells.size(), dof.get_triangulation().n_active_cells());
1596 }
1597
1598 const auto n_owned_dofs = dof.n_locally_owned_dofs();
1599
1600 // Actually, we compute the inverse of the reordering vector, called reverse
1601 // here. Later, its inverse is computed into new_indices, which is the
1602 // return argument.
1603
1604 AssertDimension(new_indices.size(), n_owned_dofs);
1605 AssertDimension(reverse.size(), n_owned_dofs);
1606
1607 // For continuous elements, we must make sure, that each dof is reordered
1608 // only once.
1609 std::vector<bool> already_sorted(n_owned_dofs, false);
1610 std::vector<types::global_dof_index> cell_dofs;
1611
1612 const auto &owned_dofs = dof.locally_owned_dofs();
1613
1614 unsigned int index = 0;
1615
1616 for (const auto &cell : cells)
1617 {
1618 // Determine the number of dofs on this cell and reinit the
1619 // vector storing these numbers.
1620 const unsigned int n_cell_dofs = cell->get_fe().n_dofs_per_cell();
1621 cell_dofs.resize(n_cell_dofs);
1622
1623 cell->get_active_or_mg_dof_indices(cell_dofs);
1624
1625 // Sort here to make sure that degrees of freedom inside a single cell
1626 // are in the same order after renumbering.
1627 std::sort(cell_dofs.begin(), cell_dofs.end());
1628
1629 for (const auto dof : cell_dofs)
1630 {
1631 const auto local_dof = owned_dofs.index_within_set(dof);
1632 if (local_dof != numbers::invalid_dof_index &&
1633 !already_sorted[local_dof])
1634 {
1635 already_sorted[local_dof] = true;
1636 reverse[index++] = local_dof;
1637 }
1638 }
1639 }
1640 Assert(index == n_owned_dofs,
1641 ExcMessage(
1642 "Traversing over the given set of cells did not cover all "
1643 "degrees of freedom in the DoFHandler. Does the set of cells "
1644 "not include all active cells?"));
1645
1646 for (types::global_dof_index i = 0; i < reverse.size(); ++i)
1647 new_indices[reverse[i]] = owned_dofs.nth_index_in_set(i);
1648 }
1649
1650
1651
1652 template <int dim, int spacedim>
1653 void
1655 const unsigned int level,
1656 const typename std::vector<
1658 {
1661
1662 std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1663 std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1664
1665 compute_cell_wise(renumbering, reverse, dof, level, cells);
1666 dof.renumber_dofs(level, renumbering);
1667 }
1668
1669
1670
1671 template <int dim, int spacedim>
1672 void
1674 std::vector<types::global_dof_index> &new_order,
1675 std::vector<types::global_dof_index> &reverse,
1676 const DoFHandler<dim, spacedim> &dof,
1677 const unsigned int level,
1678 const typename std::vector<
1680 {
1681 Assert(cells.size() == dof.get_triangulation().n_cells(level),
1682 ExcDimensionMismatch(cells.size(),
1684 Assert(new_order.size() == dof.n_dofs(level),
1685 ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
1686 Assert(reverse.size() == dof.n_dofs(level),
1687 ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
1688
1689 unsigned int n_global_dofs = dof.n_dofs(level);
1690 unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
1691
1692 std::vector<bool> already_sorted(n_global_dofs, false);
1693 std::vector<types::global_dof_index> cell_dofs(n_cell_dofs);
1694
1695 unsigned int global_index = 0;
1696
1697 for (const auto &cell : cells)
1698 {
1699 Assert(cell->level() == static_cast<int>(level), ExcInternalError());
1700
1701 cell->get_active_or_mg_dof_indices(cell_dofs);
1702 std::sort(cell_dofs.begin(), cell_dofs.end());
1703
1704 for (unsigned int i = 0; i < n_cell_dofs; ++i)
1705 {
1706 if (!already_sorted[cell_dofs[i]])
1707 {
1708 already_sorted[cell_dofs[i]] = true;
1709 reverse[global_index++] = cell_dofs[i];
1710 }
1711 }
1712 }
1713 Assert(global_index == n_global_dofs,
1714 ExcMessage(
1715 "Traversing over the given set of cells did not cover all "
1716 "degrees of freedom in the DoFHandler. Does the set of cells "
1717 "not include all cells of the specified level?"));
1718
1719 for (types::global_dof_index i = 0; i < new_order.size(); ++i)
1720 new_order[reverse[i]] = i;
1721 }
1722
1723
1724
1725 template <int dim, int spacedim>
1726 void
1728 const Tensor<1, spacedim> &direction,
1729 const bool dof_wise_renumbering)
1730 {
1731 std::vector<types::global_dof_index> renumbering(dof.n_dofs());
1732 std::vector<types::global_dof_index> reverse(dof.n_dofs());
1734 renumbering, reverse, dof, direction, dof_wise_renumbering);
1735
1736 dof.renumber_dofs(renumbering);
1737 }
1738
1739
1740
1741 template <int dim, int spacedim>
1742 void
1743 compute_downstream(std::vector<types::global_dof_index> &new_indices,
1744 std::vector<types::global_dof_index> &reverse,
1745 const DoFHandler<dim, spacedim> &dof,
1746 const Tensor<1, spacedim> &direction,
1747 const bool dof_wise_renumbering)
1748 {
1750 &dof.get_triangulation()) == nullptr),
1752
1753 if (dof_wise_renumbering == false)
1754 {
1755 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
1756 ordered_cells;
1757 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
1758 const CompareDownstream<
1760 spacedim>
1761 comparator(direction);
1762
1763 for (const auto &cell : dof.active_cell_iterators())
1764 ordered_cells.push_back(cell);
1765
1766 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
1767
1768 compute_cell_wise(new_indices, reverse, dof, ordered_cells);
1769 }
1770 else
1771 {
1772 // similar code as for
1773 // DoFTools::map_dofs_to_support_points, but
1774 // need to do this for general DoFHandler<dim, spacedim> classes and
1775 // want to be able to sort the result
1776 // (otherwise, could use something like
1777 // DoFTools::map_support_points_to_dofs)
1778 const unsigned int n_dofs = dof.n_dofs();
1779 std::vector<std::pair<Point<spacedim>, unsigned int>>
1780 support_point_list(n_dofs);
1781
1782 const hp::FECollection<dim> &fe_collection = dof.get_fe_collection();
1783 Assert(fe_collection[0].has_support_points(),
1785 hp::QCollection<dim> quadrature_collection;
1786 for (unsigned int comp = 0; comp < fe_collection.size(); ++comp)
1787 {
1788 Assert(fe_collection[comp].has_support_points(),
1790 quadrature_collection.push_back(
1791 Quadrature<dim>(fe_collection[comp].get_unit_support_points()));
1792 }
1793 hp::FEValues<dim, spacedim> hp_fe_values(fe_collection,
1794 quadrature_collection,
1796
1797 std::vector<bool> already_touched(n_dofs, false);
1798
1799 std::vector<types::global_dof_index> local_dof_indices;
1800
1801 for (const auto &cell : dof.active_cell_iterators())
1802 {
1803 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1804 local_dof_indices.resize(dofs_per_cell);
1805 hp_fe_values.reinit(cell);
1806 const FEValues<dim> &fe_values =
1807 hp_fe_values.get_present_fe_values();
1808 cell->get_active_or_mg_dof_indices(local_dof_indices);
1809 const std::vector<Point<spacedim>> &points =
1810 fe_values.get_quadrature_points();
1811 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1812 if (!already_touched[local_dof_indices[i]])
1813 {
1814 support_point_list[local_dof_indices[i]].first = points[i];
1815 support_point_list[local_dof_indices[i]].second =
1816 local_dof_indices[i];
1817 already_touched[local_dof_indices[i]] = true;
1818 }
1819 }
1820
1821 ComparePointwiseDownstream<spacedim> comparator(direction);
1822 std::sort(support_point_list.begin(),
1823 support_point_list.end(),
1824 comparator);
1825 for (types::global_dof_index i = 0; i < n_dofs; ++i)
1826 new_indices[support_point_list[i].second] = i;
1827 }
1828 }
1829
1830
1831
1832 template <int dim, int spacedim>
1833 void
1835 const unsigned int level,
1836 const Tensor<1, spacedim> &direction,
1837 const bool dof_wise_renumbering)
1838 {
1839 std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1840 std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1842 renumbering, reverse, dof, level, direction, dof_wise_renumbering);
1843
1844 dof.renumber_dofs(level, renumbering);
1845 }
1846
1847
1848
1849 template <int dim, int spacedim>
1850 void
1851 compute_downstream(std::vector<types::global_dof_index> &new_indices,
1852 std::vector<types::global_dof_index> &reverse,
1853 const DoFHandler<dim, spacedim> &dof,
1854 const unsigned int level,
1855 const Tensor<1, spacedim> &direction,
1856 const bool dof_wise_renumbering)
1857 {
1858 if (dof_wise_renumbering == false)
1859 {
1860 std::vector<typename DoFHandler<dim, spacedim>::level_cell_iterator>
1861 ordered_cells;
1862 ordered_cells.reserve(dof.get_triangulation().n_cells(level));
1863 const CompareDownstream<
1865 spacedim>
1866 comparator(direction);
1867
1869 dof.begin(level);
1871 dof.end(level);
1872
1873 while (p != end)
1874 {
1875 ordered_cells.push_back(p);
1876 ++p;
1877 }
1878 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
1879
1880 compute_cell_wise(new_indices, reverse, dof, level, ordered_cells);
1881 }
1882 else
1883 {
1886 const unsigned int n_dofs = dof.n_dofs(level);
1887 std::vector<std::pair<Point<spacedim>, unsigned int>>
1888 support_point_list(n_dofs);
1889
1891 FEValues<dim, spacedim> fe_values(dof.get_fe(),
1892 q_dummy,
1894
1895 std::vector<bool> already_touched(dof.n_dofs(), false);
1896
1897 const unsigned int dofs_per_cell = dof.get_fe().n_dofs_per_cell();
1898 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1900 dof.begin(level);
1902 dof.end(level);
1903 for (; begin != end; ++begin)
1904 {
1906 &begin_tria = begin;
1907 begin->get_active_or_mg_dof_indices(local_dof_indices);
1908 fe_values.reinit(begin_tria);
1909 const std::vector<Point<spacedim>> &points =
1910 fe_values.get_quadrature_points();
1911 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1912 if (!already_touched[local_dof_indices[i]])
1913 {
1914 support_point_list[local_dof_indices[i]].first = points[i];
1915 support_point_list[local_dof_indices[i]].second =
1916 local_dof_indices[i];
1917 already_touched[local_dof_indices[i]] = true;
1918 }
1919 }
1920
1921 ComparePointwiseDownstream<spacedim> comparator(direction);
1922 std::sort(support_point_list.begin(),
1923 support_point_list.end(),
1924 comparator);
1925 for (types::global_dof_index i = 0; i < n_dofs; ++i)
1926 new_indices[support_point_list[i].second] = i;
1927 }
1928 }
1929
1930
1931
1935 namespace internal
1936 {
1937 template <int dim>
1939 {
1948
1953 : center(center)
1954 , counter(counter)
1955 {}
1956
1960 template <class DHCellIterator>
1961 bool
1962 operator()(const DHCellIterator &c1, const DHCellIterator &c2) const
1963 {
1964 // dispatch to
1965 // dimension-dependent functions
1966 return compare(c1, c2, std::integral_constant<int, dim>());
1967 }
1968
1969 private:
1973 template <class DHCellIterator, int xdim>
1974 bool
1975 compare(const DHCellIterator &c1,
1976 const DHCellIterator &c2,
1977 std::integral_constant<int, xdim>) const
1978 {
1979 const Tensor<1, dim> v1 = c1->center() - center;
1980 const Tensor<1, dim> v2 = c2->center() - center;
1981 const double s1 = std::atan2(v1[0], v1[1]);
1982 const double s2 = std::atan2(v2[0], v2[1]);
1983 return (counter ? (s1 > s2) : (s2 > s1));
1984 }
1985
1986
1991 template <class DHCellIterator>
1992 bool
1993 compare(const DHCellIterator &,
1994 const DHCellIterator &,
1995 std::integral_constant<int, 1>) const
1996 {
1997 Assert(dim >= 2,
1998 ExcMessage("This operation only makes sense for dim>=2."));
1999 return false;
2000 }
2001 };
2002 } // namespace internal
2003
2004
2005
2006 template <int dim, int spacedim>
2007 void
2009 const Point<spacedim> &center,
2010 const bool counter)
2011 {
2012 std::vector<types::global_dof_index> renumbering(dof.n_dofs());
2013 compute_clockwise_dg(renumbering, dof, center, counter);
2014
2015 dof.renumber_dofs(renumbering);
2016 }
2017
2018
2019
2020 template <int dim, int spacedim>
2021 void
2022 compute_clockwise_dg(std::vector<types::global_dof_index> &new_indices,
2023 const DoFHandler<dim, spacedim> &dof,
2024 const Point<spacedim> &center,
2025 const bool counter)
2026 {
2027 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
2028 ordered_cells;
2029 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
2030 internal::ClockCells<spacedim> comparator(center, counter);
2031
2032 for (const auto &cell : dof.active_cell_iterators())
2033 ordered_cells.push_back(cell);
2034
2035 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
2036
2037 std::vector<types::global_dof_index> reverse(new_indices.size());
2038 compute_cell_wise(new_indices, reverse, dof, ordered_cells);
2039 }
2040
2041
2042
2043 template <int dim, int spacedim>
2044 void
2046 const unsigned int level,
2047 const Point<spacedim> &center,
2048 const bool counter)
2049 {
2050 std::vector<typename DoFHandler<dim, spacedim>::level_cell_iterator>
2051 ordered_cells;
2052 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
2053 internal::ClockCells<spacedim> comparator(center, counter);
2054
2056 dof.begin(level);
2058 dof.end(level);
2059
2060 while (p != end)
2061 {
2062 ordered_cells.push_back(p);
2063 ++p;
2064 }
2065 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
2066
2067 cell_wise(dof, level, ordered_cells);
2068 }
2069
2070
2071
2072 template <int dim, int spacedim>
2073 void
2075 {
2076 std::vector<types::global_dof_index> renumbering(
2077 dof_handler.n_dofs(), numbers::invalid_dof_index);
2078 compute_random(renumbering, dof_handler);
2079
2080 dof_handler.renumber_dofs(renumbering);
2081 }
2082
2083
2084
2085 template <int dim, int spacedim>
2086 void
2087 random(DoFHandler<dim, spacedim> &dof_handler, const unsigned int level)
2088 {
2091
2092 std::vector<types::global_dof_index> renumbering(
2093 dof_handler.locally_owned_mg_dofs(level).n_elements(),
2095
2096 compute_random(renumbering, dof_handler, level);
2097
2098 dof_handler.renumber_dofs(level, renumbering);
2099 }
2100
2101
2102
2103 template <int dim, int spacedim>
2104 void
2105 compute_random(std::vector<types::global_dof_index> &new_indices,
2106 const DoFHandler<dim, spacedim> &dof_handler)
2107 {
2108 const types::global_dof_index n_dofs = dof_handler.n_dofs();
2109 Assert(new_indices.size() == n_dofs,
2110 ExcDimensionMismatch(new_indices.size(), n_dofs));
2111
2112 std::iota(new_indices.begin(),
2113 new_indices.end(),
2115
2116 // shuffle the elements; the following is essentially std::shuffle (which
2117 // is new in C++11) but with a boost URNG
2118 // we could use std::mt19937 here but doing so results in compiler-dependent
2119 // output
2120 ::boost::mt19937 random_number_generator;
2121 for (unsigned int i = 1; i < n_dofs; ++i)
2122 {
2123 // get a random number between 0 and i (inclusive)
2124 const unsigned int j =
2125 ::boost::random::uniform_int_distribution<>(0, i)(
2126 random_number_generator);
2127
2128 // if possible, swap the elements
2129 if (i != j)
2130 std::swap(new_indices[i], new_indices[j]);
2131 }
2132 }
2133
2134
2135
2136 template <int dim, int spacedim>
2137 void
2138 compute_random(std::vector<types::global_dof_index> &new_indices,
2139 const DoFHandler<dim, spacedim> &dof_handler,
2140 const unsigned int level)
2141 {
2142 const types::global_dof_index n_dofs = dof_handler.n_dofs(level);
2143 Assert(new_indices.size() == n_dofs,
2144 ExcDimensionMismatch(new_indices.size(), n_dofs));
2145
2146 std::iota(new_indices.begin(),
2147 new_indices.end(),
2149
2150 // shuffle the elements; the following is essentially std::shuffle (which
2151 // is new in C++11) but with a boost URNG
2152 // we could use std::mt19937 here but doing so results in
2153 // compiler-dependent output
2154 ::boost::mt19937 random_number_generator;
2155 for (unsigned int i = 1; i < n_dofs; ++i)
2156 {
2157 // get a random number between 0 and i (inclusive)
2158 const unsigned int j =
2159 ::boost::random::uniform_int_distribution<>(0, i)(
2160 random_number_generator);
2161
2162 // if possible, swap the elements
2163 if (i != j)
2164 std::swap(new_indices[i], new_indices[j]);
2165 }
2166 }
2167
2168
2169
2170 template <int dim, int spacedim>
2171 void
2173 {
2174 Assert(
2175 (!dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
2176 &dof_handler.get_triangulation())),
2177 ExcMessage(
2178 "Parallel triangulations are already enumerated according to their MPI process id."));
2179
2180 std::vector<types::global_dof_index> renumbering(
2181 dof_handler.n_dofs(), numbers::invalid_dof_index);
2182 compute_subdomain_wise(renumbering, dof_handler);
2183
2184 dof_handler.renumber_dofs(renumbering);
2185 }
2186
2187
2188
2189 template <int dim, int spacedim>
2190 void
2191 compute_subdomain_wise(std::vector<types::global_dof_index> &new_dof_indices,
2192 const DoFHandler<dim, spacedim> &dof_handler)
2193 {
2194 const types::global_dof_index n_dofs = dof_handler.n_dofs();
2195 Assert(new_dof_indices.size() == n_dofs,
2196 ExcDimensionMismatch(new_dof_indices.size(), n_dofs));
2197
2198 // first get the association of each dof
2199 // with a subdomain and determine the total
2200 // number of subdomain ids used
2201 std::vector<types::subdomain_id> subdomain_association(n_dofs);
2202 DoFTools::get_subdomain_association(dof_handler, subdomain_association);
2203 const unsigned int n_subdomains =
2204 *std::max_element(subdomain_association.begin(),
2205 subdomain_association.end()) +
2206 1;
2207
2208 // then renumber the subdomains by first
2209 // looking at those belonging to subdomain
2210 // 0, then those of subdomain 1, etc. note
2211 // that the algorithm is stable, i.e. if
2212 // two dofs i,j have i<j and belong to the
2213 // same subdomain, then they will be in
2214 // this order also after reordering
2215 std::fill(new_dof_indices.begin(),
2216 new_dof_indices.end(),
2218 types::global_dof_index next_free_index = 0;
2219 for (types::subdomain_id subdomain = 0; subdomain < n_subdomains;
2220 ++subdomain)
2221 for (types::global_dof_index i = 0; i < n_dofs; ++i)
2222 if (subdomain_association[i] == subdomain)
2223 {
2224 Assert(new_dof_indices[i] == numbers::invalid_dof_index,
2226 new_dof_indices[i] = next_free_index;
2227 ++next_free_index;
2228 }
2229
2230 // we should have numbered all dofs
2231 Assert(next_free_index == n_dofs, ExcInternalError());
2232 Assert(std::find(new_dof_indices.begin(),
2233 new_dof_indices.end(),
2234 numbers::invalid_dof_index) == new_dof_indices.end(),
2236 }
2237
2238
2239
2240 template <int dim, int spacedim>
2241 void
2243 {
2244 std::vector<types::global_dof_index> renumbering(
2246 compute_support_point_wise(renumbering, dof_handler);
2247
2248 // If there is only one component then there is nothing to do, so check
2249 // first:
2250 if (Utilities::MPI::max(renumbering.size(),
2251 dof_handler.get_communicator()) > 0)
2252 dof_handler.renumber_dofs(renumbering);
2253 }
2254
2255
2256
2257 template <int dim, int spacedim>
2258 void
2260 std::vector<types::global_dof_index> &new_dof_indices,
2261 const DoFHandler<dim, spacedim> &dof_handler)
2262 {
2263 const types::global_dof_index n_dofs = dof_handler.n_locally_owned_dofs();
2264 Assert(new_dof_indices.size() == n_dofs,
2265 ExcDimensionMismatch(new_dof_indices.size(), n_dofs));
2266
2267 // This renumbering occurs in three steps:
2268 // 1. Compute the component-wise renumbering so that all DoFs of component
2269 // i are less than the DoFs of component i + 1.
2270 // 2. Compute a second renumbering component_to_nodal in which the
2271 // renumbering is now, for two components, [u0, v0, u1, v1, ...]: i.e.,
2272 // DoFs are first sorted by component and then by support point.
2273 // 3. Compose the two renumberings to obtain the final result.
2274
2275 // Step 1:
2276 std::vector<types::global_dof_index> component_renumbering(
2278 compute_component_wise<dim, spacedim>(component_renumbering,
2279 dof_handler.begin_active(),
2280 dof_handler.end(),
2281 std::vector<unsigned int>(),
2282 false);
2283
2284#ifdef DEBUG
2285 {
2286 const std::vector<types::global_dof_index> dofs_per_component =
2287 DoFTools::count_dofs_per_fe_component(dof_handler, true);
2288 for (const auto &dpc : dofs_per_component)
2289 Assert(dofs_per_component[0] == dpc, ExcNotImplemented());
2290 }
2291#endif
2292 const unsigned int n_components =
2293 dof_handler.get_fe_collection().n_components();
2294 Assert(dof_handler.n_dofs() % n_components == 0, ExcInternalError());
2295 const types::global_dof_index dofs_per_component =
2296 dof_handler.n_dofs() / n_components;
2297 const types::global_dof_index local_dofs_per_component =
2298 dof_handler.n_locally_owned_dofs() / n_components;
2299
2300 // At this point we have no more communication to do - simplify things by
2301 // returning early if possible
2302 if (component_renumbering.empty())
2303 {
2304 new_dof_indices.resize(0);
2305 return;
2306 }
2307 std::fill(new_dof_indices.begin(),
2308 new_dof_indices.end(),
2310 // This index set equals what dof_handler.locally_owned_dofs() would be if
2311 // we executed the componentwise renumbering.
2312 IndexSet component_renumbered_dofs(dof_handler.n_dofs());
2313 // DoFs in each component are now consecutive, which IndexSet::add_indices()
2314 // can exploit by avoiding calls to sort. Make use of that by adding DoFs
2315 // one component at a time:
2316 std::vector<types::global_dof_index> component_dofs(
2317 local_dofs_per_component);
2318 for (unsigned int component = 0; component < n_components; ++component)
2319 {
2320 for (std::size_t i = 0; i < local_dofs_per_component; ++i)
2321 component_dofs[i] =
2322 component_renumbering[n_components * i + component];
2323 component_renumbered_dofs.add_indices(component_dofs.begin(),
2324 component_dofs.end());
2325 }
2326 component_renumbered_dofs.compress();
2327#ifdef DEBUG
2328 {
2329 IndexSet component_renumbered_dofs2(dof_handler.n_dofs());
2330 component_renumbered_dofs2.add_indices(component_renumbering.begin(),
2331 component_renumbering.end());
2332 Assert(component_renumbered_dofs2 == component_renumbered_dofs,
2334 }
2335#endif
2336 for (const FiniteElement<dim, spacedim> &fe :
2337 dof_handler.get_fe_collection())
2338 {
2339 AssertThrow(fe.dofs_per_cell == 0 || fe.has_support_points(),
2341 for (unsigned int i = 0; i < fe.n_base_elements(); ++i)
2343 fe.base_element(0).get_unit_support_points() ==
2344 fe.base_element(i).get_unit_support_points(),
2345 ExcMessage(
2346 "All base elements should have the same support points."));
2347 }
2348
2349 std::vector<types::global_dof_index> component_to_nodal(
2351
2352 // Step 2:
2353 std::vector<types::global_dof_index> cell_dofs;
2354 std::vector<types::global_dof_index> component_renumbered_cell_dofs;
2355 const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
2356 // Reuse the original index space for the new renumbering: it is the right
2357 // size and is contiguous on the current processor
2358 auto next_dof_it = locally_owned_dofs.begin();
2359 for (const auto &cell : dof_handler.active_cell_iterators())
2360 if (cell->is_locally_owned())
2361 {
2362 const FiniteElement<dim, spacedim> &fe = cell->get_fe();
2363 cell_dofs.resize(fe.dofs_per_cell);
2364 component_renumbered_cell_dofs.resize(fe.dofs_per_cell);
2365 cell->get_dof_indices(cell_dofs);
2366 // Apply the component renumbering while skipping any ghost dofs. This
2367 // algorithm assumes that all locally owned DoFs before the component
2368 // renumbering are still locally owned afterwards (just with a new
2369 // index).
2370 for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
2371 {
2372 if (locally_owned_dofs.is_element(cell_dofs[i]))
2373 {
2374 const auto local_index =
2375 locally_owned_dofs.index_within_set(cell_dofs[i]);
2376 component_renumbered_cell_dofs[i] =
2377 component_renumbering[local_index];
2378 }
2379 else
2380 {
2381 component_renumbered_cell_dofs[i] =
2383 }
2384 }
2385
2386 for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
2387 {
2388 if (fe.system_to_component_index(i).first == 0 &&
2389 component_renumbered_dofs.is_element(
2390 component_renumbered_cell_dofs[i]))
2391 {
2392 for (unsigned int component = 0;
2393 component < fe.n_components();
2394 ++component)
2395 {
2396 // Since we are indexing in an odd way here it is much
2397 // simpler to compute the permutation separately and
2398 // combine it at the end instead of doing both at once
2399 const auto local_index =
2400 component_renumbered_dofs.index_within_set(
2401 component_renumbered_cell_dofs[i] +
2402 dofs_per_component * component);
2403
2404 if (component_to_nodal[local_index] ==
2406 {
2407 component_to_nodal[local_index] = *next_dof_it;
2408 ++next_dof_it;
2409 }
2410 }
2411 }
2412 }
2413 }
2414
2415 // Step 3:
2416 for (std::size_t i = 0; i < dof_handler.n_locally_owned_dofs(); ++i)
2417 {
2418 const auto local_index =
2419 component_renumbered_dofs.index_within_set(component_renumbering[i]);
2420 new_dof_indices[i] = component_to_nodal[local_index];
2421 }
2422 }
2423
2424
2425
2426 template <int dim,
2427 int spacedim,
2428 typename Number,
2429 typename VectorizedArrayType>
2430 void
2432 DoFHandler<dim, spacedim> &dof_handler,
2434 {
2435 const std::vector<types::global_dof_index> new_global_numbers =
2436 compute_matrix_free_data_locality(dof_handler, matrix_free);
2437 if (matrix_free.get_mg_level() == ::numbers::invalid_unsigned_int)
2438 dof_handler.renumber_dofs(new_global_numbers);
2439 else
2440 dof_handler.renumber_dofs(matrix_free.get_mg_level(), new_global_numbers);
2441 }
2442
2443
2444
2445 template <int dim, int spacedim, typename Number, typename AdditionalDataType>
2446 void
2448 const AffineConstraints<Number> &constraints,
2449 const AdditionalDataType &matrix_free_data)
2450 {
2451 const std::vector<types::global_dof_index> new_global_numbers =
2453 constraints,
2454 matrix_free_data);
2455 if (matrix_free_data.mg_level == ::numbers::invalid_unsigned_int)
2456 dof_handler.renumber_dofs(new_global_numbers);
2457 else
2458 dof_handler.renumber_dofs(matrix_free_data.mg_level, new_global_numbers);
2459 }
2460
2461
2462
2463 template <int dim, int spacedim, typename Number, typename AdditionalDataType>
2464 std::vector<types::global_dof_index>
2466 const DoFHandler<dim, spacedim> &dof_handler,
2467 const AffineConstraints<Number> &constraints,
2468 const AdditionalDataType &matrix_free_data)
2469 {
2470 AdditionalDataType my_mf_data = matrix_free_data;
2471 my_mf_data.initialize_mapping = false;
2472 my_mf_data.tasks_parallel_scheme = AdditionalDataType::none;
2473
2474 typename AdditionalDataType::MatrixFreeType separate_matrix_free;
2475 separate_matrix_free.reinit(::MappingQ1<dim>(),
2476 dof_handler,
2477 constraints,
2478 ::QGauss<1>(2),
2479 my_mf_data);
2480 return compute_matrix_free_data_locality(dof_handler, separate_matrix_free);
2481 }
2482
2483
2484
2485 // Implementation details for matrix-free renumbering
2486 namespace
2487 {
2488 // Compute a vector of lists with number of unknowns of the same category
2489 // in terms of influence from other MPI processes, starting with unknowns
2490 // touched only by the local process and finally a new set of indices
2491 // going to different MPI neighbors. Later passes of the algorithm will
2492 // re-order unknowns within each of these sets.
2493 std::vector<std::vector<unsigned int>>
2494 group_dofs_by_rank_access(
2495 const ::Utilities::MPI::Partitioner &partitioner)
2496 {
2497 // count the number of times a locally owned DoF is accessed by the
2498 // remote ghost data
2499 std::vector<unsigned int> touch_count(partitioner.locally_owned_size());
2500 for (const auto &p : partitioner.import_indices())
2501 for (unsigned int i = p.first; i < p.second; ++i)
2502 touch_count[i]++;
2503
2504 // category 0: DoFs never touched by ghosts
2505 std::vector<std::vector<unsigned int>> result(1);
2506 for (unsigned int i = 0; i < touch_count.size(); ++i)
2507 if (touch_count[i] == 0)
2508 result.back().push_back(i);
2509
2510 // DoFs with 1 appearance can be simply categorized by their (single)
2511 // MPI rank, whereas we need to go an extra round for the remaining DoFs
2512 // by collecting the owning processes by hand
2513 std::map<unsigned int, std::vector<unsigned int>>
2514 multiple_ranks_access_dof;
2515 const std::vector<std::pair<unsigned int, unsigned int>> &import_targets =
2516 partitioner.import_targets();
2517 auto it = partitioner.import_indices().begin();
2518 for (const std::pair<unsigned int, unsigned int> &proc : import_targets)
2519 {
2520 result.emplace_back();
2521 unsigned int count_dofs = 0;
2522 while (count_dofs < proc.second)
2523 {
2524 for (unsigned int i = it->first; i < it->second;
2525 ++i, ++count_dofs)
2526 {
2527 if (touch_count[i] == 1)
2528 result.back().push_back(i);
2529 else
2530 multiple_ranks_access_dof[i].push_back(proc.first);
2531 }
2532 ++it;
2533 }
2534 }
2535 Assert(it == partitioner.import_indices().end(), ExcInternalError());
2536
2537 // Now go from the computed map on DoFs to a map on the processes for
2538 // DoFs with multiple owners, and append the DoFs found this way to our
2539 // global list
2540 std::map<std::vector<unsigned int>,
2541 std::vector<unsigned int>,
2542 std::function<bool(const std::vector<unsigned int> &,
2543 const std::vector<unsigned int> &)>>
2544 dofs_by_rank{[](const std::vector<unsigned int> &a,
2545 const std::vector<unsigned int> &b) {
2546 if (a.size() < b.size())
2547 return true;
2548 if (a.size() == b.size())
2549 {
2550 for (unsigned int i = 0; i < a.size(); ++i)
2551 if (a[i] < b[i])
2552 return true;
2553 else if (a[i] > b[i])
2554 return false;
2555 }
2556 return false;
2557 }};
2558 for (const auto &entry : multiple_ranks_access_dof)
2559 dofs_by_rank[entry.second].push_back(entry.first);
2560
2561 for (const auto &procs : dofs_by_rank)
2562 result.push_back(procs.second);
2563
2564 return result;
2565 }
2566
2567
2568
2569 // Compute two vectors, the first indicating the best numbers for a
2570 // MatrixFree::cell_loop and the second the count of how often a DoF is
2571 // touched by different cell groups, in order to later figure out DoFs
2572 // with far reach and those with only local influence.
2573 template <int dim, typename Number, typename VectorizedArrayType>
2574 std::pair<std::vector<unsigned int>, std::vector<unsigned char>>
2575 compute_mf_numbering(
2577 const unsigned int component)
2578 {
2579 const IndexSet &owned_dofs = matrix_free.get_dof_info(component)
2580 .vector_partitioner->locally_owned_range();
2581 const unsigned int n_comp =
2582 matrix_free.get_dof_handler(component).get_fe().n_components();
2583 Assert(
2584 matrix_free.get_dof_handler(component).get_fe().n_base_elements() == 1,
2586 Assert(dynamic_cast<const FE_Q_Base<dim> *>(
2587 &matrix_free.get_dof_handler(component).get_fe().base_element(
2588 0)),
2589 ExcNotImplemented("Matrix-free renumbering only works for "
2590 "FE_Q elements"));
2591
2592 const unsigned int fe_degree =
2593 matrix_free.get_dof_handler(component).get_fe().degree;
2594 const unsigned int nn = fe_degree - 1;
2595
2596 // Data structure used further down for the identification of various
2597 // entities in the hierarchical numbering of unknowns. The first number
2598 // indicates the offset from which a given object starts its range of
2599 // numbers in the hierarchical DoF numbering of FE_Q, and the second the
2600 // number of unknowns per component on that particular component. The
2601 // numbers are group by the 3^dim possible objects, listed in
2602 // lexicographic order.
2603 std::array<std::pair<unsigned int, unsigned int>,
2604 ::Utilities::pow(3, dim)>
2605 dofs_on_objects;
2606 if (dim == 1)
2607 {
2608 dofs_on_objects[0] = std::make_pair(0U, 1U);
2609 dofs_on_objects[1] = std::make_pair(2 * n_comp, nn);
2610 dofs_on_objects[2] = std::make_pair(n_comp, 1U);
2611 }
2612 else if (dim == 2)
2613 {
2614 dofs_on_objects[0] = std::make_pair(0U, 1U);
2615 dofs_on_objects[1] = std::make_pair(n_comp * (4 + 2 * nn), nn);
2616 dofs_on_objects[2] = std::make_pair(n_comp, 1U);
2617 dofs_on_objects[3] = std::make_pair(n_comp * 4, nn);
2618 dofs_on_objects[4] = std::make_pair(n_comp * (4 + 4 * nn), nn * nn);
2619 dofs_on_objects[5] = std::make_pair(n_comp * (4 + 1 * nn), nn);
2620 dofs_on_objects[6] = std::make_pair(2 * n_comp, 1U);
2621 dofs_on_objects[7] = std::make_pair(n_comp * (4 + 3 * nn), nn);
2622 dofs_on_objects[8] = std::make_pair(3 * n_comp, 1U);
2623 }
2624 else if (dim == 3)
2625 {
2626 dofs_on_objects[0] = std::make_pair(0U, 1U);
2627 dofs_on_objects[1] = std::make_pair(n_comp * (8 + 2 * nn), nn);
2628 dofs_on_objects[2] = std::make_pair(n_comp, 1U);
2629 dofs_on_objects[3] = std::make_pair(n_comp * 8, nn);
2630 dofs_on_objects[4] =
2631 std::make_pair(n_comp * (8 + 12 * nn + 4 * nn * nn), nn * nn);
2632 dofs_on_objects[5] = std::make_pair(n_comp * (8 + 1 * nn), nn);
2633 dofs_on_objects[6] = std::make_pair(n_comp * 2, 1U);
2634 dofs_on_objects[7] = std::make_pair(n_comp * (8 + 3 * nn), nn);
2635 dofs_on_objects[8] = std::make_pair(n_comp * 3, 1U);
2636 dofs_on_objects[9] = std::make_pair(n_comp * (8 + 8 * nn), nn);
2637 dofs_on_objects[10] =
2638 std::make_pair(n_comp * (8 + 12 * nn + 2 * nn * nn), nn * nn);
2639 dofs_on_objects[11] = std::make_pair(n_comp * (8 + 9 * nn), nn);
2640 dofs_on_objects[12] = std::make_pair(n_comp * (8 + 12 * nn), nn * nn);
2641 dofs_on_objects[13] =
2642 std::make_pair(n_comp * (8 + 12 * nn + 6 * nn * nn), nn * nn * nn);
2643 dofs_on_objects[14] =
2644 std::make_pair(n_comp * (8 + 12 * nn + 1 * nn * nn), nn * nn);
2645 dofs_on_objects[15] = std::make_pair(n_comp * (8 + 10 * nn), nn);
2646 dofs_on_objects[16] =
2647 std::make_pair(n_comp * (8 + 12 * nn + 3 * nn * nn), nn * nn);
2648 dofs_on_objects[17] = std::make_pair(n_comp * (8 + 11 * nn), nn);
2649 dofs_on_objects[18] = std::make_pair(n_comp * 4, 1U);
2650 dofs_on_objects[19] = std::make_pair(n_comp * (8 + 6 * nn), nn);
2651 dofs_on_objects[20] = std::make_pair(n_comp * 5, 1U);
2652 dofs_on_objects[21] = std::make_pair(n_comp * (8 + 4 * nn), nn);
2653 dofs_on_objects[22] =
2654 std::make_pair(n_comp * (8 + 12 * nn + 5 * nn * nn), nn * nn);
2655 dofs_on_objects[23] = std::make_pair(n_comp * (8 + 5 * nn), nn);
2656 dofs_on_objects[24] = std::make_pair(n_comp * 6, 1U);
2657 dofs_on_objects[25] = std::make_pair(n_comp * (8 + 7 * nn), nn);
2658 dofs_on_objects[26] = std::make_pair(n_comp * 7, 1U);
2659 }
2660
2661 const auto renumber_func = [](const types::global_dof_index dof_index,
2662 const IndexSet &owned_dofs,
2663 std::vector<unsigned int> &result,
2664 unsigned int &counter_dof_numbers) {
2665 const types::global_dof_index local_dof_index =
2666 owned_dofs.index_within_set(dof_index);
2667 if (local_dof_index != numbers::invalid_dof_index)
2668 {
2669 AssertIndexRange(local_dof_index, result.size());
2670 if (result[local_dof_index] == numbers::invalid_unsigned_int)
2671 result[local_dof_index] = counter_dof_numbers++;
2672 }
2673 };
2674
2675 unsigned int counter_dof_numbers = 0;
2676 std::vector<unsigned int> dofs_extracted;
2677 std::vector<types::global_dof_index> dof_indices(
2678 matrix_free.get_dof_handler(component).get_fe().dofs_per_cell);
2679
2680 // We now define a lambda function that does two things: (a) determine
2681 // DoF numbers in a way that fit with the order a MatrixFree loop
2682 // travels through the cells (variable 'dof_numbers_mf_order'), and (b)
2683 // determine which unknowns are only touched from within a single range
2684 // of cells and which ones span multiple ranges (variable
2685 // 'touch_count'). Note that this process is done by calling into
2686 // MatrixFree::cell_loop, which gives the right level of granularity
2687 // (when executed in serial) for the scheduled vector operations. Note
2688 // that we pick the unconstrained indices in the hierarchical order for
2689 // part (a) as this makes it easy to identify the DoFs on the individual
2690 // entities, whereas we pick the numbers with constraints eliminated for
2691 // part (b). For the latter, we keep track of each DoF's interaction
2692 // with different ranges of cell batches, i.e., call-backs into the
2693 // operation_on_cell_range() function.
2694 const unsigned int n_owned_dofs = owned_dofs.n_elements();
2695 std::vector<unsigned int> dof_numbers_mf_order(
2696 n_owned_dofs, ::numbers::invalid_unsigned_int);
2697 std::vector<unsigned int> last_touch_by_cell_batch_range(
2698 n_owned_dofs, ::numbers::invalid_unsigned_int);
2699 std::vector<unsigned char> touch_count(n_owned_dofs);
2700
2701 const auto operation_on_cell_range =
2703 unsigned int &,
2704 const unsigned int &,
2705 const std::pair<unsigned int, unsigned int> &cell_range) {
2706 for (unsigned int cell = cell_range.first; cell < cell_range.second;
2707 ++cell)
2708 {
2709 // part (a): assign beneficial numbers
2710 for (unsigned int v = 0;
2711 v < data.n_active_entries_per_cell_batch(cell);
2712 ++v)
2713 {
2714 // get the indices for the dofs in cell_batch
2716 data.get_cell_iterator(cell, v, component)
2717 ->get_dof_indices(dof_indices);
2718 else
2719 data.get_cell_iterator(cell, v, component)
2720 ->get_mg_dof_indices(dof_indices);
2721
2722 for (unsigned int a = 0; a < dofs_on_objects.size(); ++a)
2723 {
2724 const auto &r = dofs_on_objects[a];
2725 if (a == 10 || a == 16)
2726 // switch order x-z for y faces in 3d to lexicographic
2727 // layout
2728 for (unsigned int i1 = 0; i1 < nn; ++i1)
2729 for (unsigned int i0 = 0; i0 < nn; ++i0)
2730 for (unsigned int c = 0; c < n_comp; ++c)
2731 renumber_func(dof_indices[r.first + r.second * c +
2732 i1 + i0 * nn],
2733 owned_dofs,
2734 dof_numbers_mf_order,
2735 counter_dof_numbers);
2736 else
2737 for (unsigned int i = 0; i < r.second; ++i)
2738 for (unsigned int c = 0; c < n_comp; ++c)
2739 renumber_func(
2740 dof_indices[r.first + r.second * c + i],
2741 owned_dofs,
2742 dof_numbers_mf_order,
2743 counter_dof_numbers);
2744 }
2745 }
2746
2747 // part (b): increment the touch count of a dof appearing in the
2748 // current cell batch if it was last touched by another than the
2749 // present cell batch range (we track them via storing the last
2750 // cell batch range that touched a particular dof)
2752 dofs_extracted, cell);
2753 for (const unsigned int dof_index : dofs_extracted)
2754 if (dof_index < n_owned_dofs &&
2755 last_touch_by_cell_batch_range[dof_index] !=
2756 cell_range.first)
2757 {
2758 ++touch_count[dof_index];
2759 last_touch_by_cell_batch_range[dof_index] =
2760 cell_range.first;
2761 }
2762 }
2763 };
2764
2765 // Finally run the matrix-free loop.
2766 Assert(matrix_free.get_task_info().scheme ==
2768 ExcNotImplemented("Renumbering only available for non-threaded "
2769 "version of MatrixFree::cell_loop"));
2770
2771 matrix_free.template cell_loop<unsigned int, unsigned int>(
2772 operation_on_cell_range, counter_dof_numbers, counter_dof_numbers);
2773
2774 AssertDimension(counter_dof_numbers, n_owned_dofs);
2775
2776 return std::make_pair(dof_numbers_mf_order, touch_count);
2777 }
2778
2779 } // namespace
2780
2781
2782
2783 template <int dim,
2784 int spacedim,
2785 typename Number,
2786 typename VectorizedArrayType>
2787 std::vector<types::global_dof_index>
2789 const DoFHandler<dim, spacedim> &dof_handler,
2791 {
2792 Assert(matrix_free.indices_initialized(),
2793 ExcMessage("You need to set up indices in MatrixFree "
2794 "to be able to compute a renumbering!"));
2795
2796 // try to locate the `DoFHandler` within the given MatrixFree object.
2797 unsigned int component = 0;
2798 for (; component < matrix_free.n_components(); ++component)
2799 if (&matrix_free.get_dof_handler(component) == &dof_handler)
2800 break;
2801
2802 Assert(component < matrix_free.n_components(),
2803 ExcMessage("Could not locate the given DoFHandler in MatrixFree"));
2804
2805 // Summary of the algorithm below:
2806 // (a) compute renumbering of each DoF in the order the corresponding
2807 // object appears in the mf loop -> local_numbering
2808 // (b) determine by how many cell groups (call-back places in the loop) a
2809 // dof is touched -> touch_count
2810 // (c) determine by how many MPI processes a dof is touched
2811 // -> dofs_by_rank_access
2812 // (d) combine both category types of (b) and (c) and list the indices
2813 // according to four categories
2814
2815 const std::vector<std::vector<unsigned int>> dofs_by_rank_access =
2816 group_dofs_by_rank_access(
2817 *matrix_free.get_dof_info(component).vector_partitioner);
2818
2819 const auto &[local_numbering, touch_count] =
2820 compute_mf_numbering(matrix_free, component);
2821
2822 const IndexSet &owned_dofs = matrix_free.get_dof_info(component)
2823 .vector_partitioner->locally_owned_range();
2824 const types::global_dof_index locally_owned_size = owned_dofs.n_elements();
2825 AssertDimension(locally_owned_size, local_numbering.size());
2826 AssertDimension(locally_owned_size, touch_count.size());
2827
2828 // Create a second permutation to group the unknowns into the following
2829 // four categories, to be eventually composed with 'local_renumbering'
2830 enum class Category : unsigned char
2831 {
2832 single, // DoFs only referring to a single batch and MPI process
2833 multiple, // DoFs referring to multiple batches (long-range)
2834 mpi_dof, // DoFs referring to DoFs in exchange with other MPI processes
2835 constrained // DoFs without any reference (constrained DoFs)
2836 };
2837 std::vector<Category> categories(locally_owned_size);
2838 for (unsigned int i = 0; i < locally_owned_size; ++i)
2839 switch (touch_count[i])
2840 {
2841 case 0:
2842 categories[local_numbering[i]] = Category::constrained;
2843 break;
2844 case 1:
2845 categories[local_numbering[i]] = Category::single;
2846 break;
2847 default:
2848 categories[local_numbering[i]] = Category::multiple;
2849 break;
2850 }
2851 for (unsigned int chunk = 1; chunk < dofs_by_rank_access.size(); ++chunk)
2852 for (const auto i : dofs_by_rank_access[chunk])
2853 categories[local_numbering[i]] = Category::mpi_dof;
2854
2855 // Assign numbers to the categories in the order listed in the enum, which
2856 // is done by starting 'counter' for each category at an offset
2857 std::array<unsigned int, 4> n_entries_per_category{};
2858 for (const Category i : categories)
2859 {
2860 AssertIndexRange(static_cast<int>(i), 4);
2861 ++n_entries_per_category[static_cast<int>(i)];
2862 }
2863 std::array<unsigned int, 4> counters{};
2864 for (unsigned int i = 1; i < 4; ++i)
2865 counters[i] = counters[i - 1] + n_entries_per_category[i - 1];
2866 std::vector<unsigned int> numbering_categories;
2867 numbering_categories.reserve(locally_owned_size);
2868 for (const Category category : categories)
2869 {
2870 numbering_categories.push_back(counters[static_cast<int>(category)]);
2871 ++counters[static_cast<int>(category)];
2872 }
2873
2874 // The final numbers are given by the composition of the two permutations
2875 std::vector<::types::global_dof_index> new_global_numbers(
2876 locally_owned_size);
2877 for (unsigned int i = 0; i < locally_owned_size; ++i)
2878 new_global_numbers[i] =
2879 owned_dofs.nth_index_in_set(numbering_categories[local_numbering[i]]);
2880
2881 return new_global_numbers;
2882 }
2883
2884} // namespace DoFRenumbering
2885
2886
2887
2888/*-------------- Explicit Instantiations -------------------------------*/
2889#include "dof_renumbering.inst"
2890
2891
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_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:1765
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1990
size_type n_elements() const
Definition index_set.h:1923
bool is_element(const size_type index) const
Definition index_set.h:1883
ElementIterator begin() const
Definition index_set.h:1699
void subtract_set(const IndexSet &other)
Definition index_set.cc:470
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1971
void compress() const
Definition index_set.h:1773
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1820
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) 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
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
bool indices_initialized() const
unsigned int n_components() const
Definition point.h:111
size_type n_rows() const
virtual MPI_Comm get_communicator() 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:264
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:694
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:293
void push_back(const Quadrature< dim_in > &new_quadrature)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > center
Point< 3 > vertices[4]
Point< 2 > second
Definition grid_out.cc:4614
Point< 2 > first
Definition grid_out.cc:4613
unsigned int level
Definition grid_out.cc:4616
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.
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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={})
if(marked_vertices.size() !=0) for(auto it
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:575
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:963
std::vector< Integer > reverse_permutation(const std::vector< Integer > &permutation)
Definition utilities.h:1681
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::global_dof_index invalid_dof_index
Definition types.h:252
typename type_identity< T >::type type_identity_t
Definition type_traits.h:95
unsigned short int fe_index
Definition types.h:59
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
void get_dof_indices_on_cell_batch(std::vector< unsigned int > &local_indices, const unsigned int cell_batch, const bool with_constraints=true) const
Definition dof_info.cc:82
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition dof_info.h:592
#define DEAL_II_DOF_INDEX_MPI_TYPE
Definition types.h:102