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