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