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