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