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