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