Reference documentation for deal.II version Git a39983f358 2020-01-20 20:47:10 -0500
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
dof_renumbering.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 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 
16 #include <deal.II/base/quadrature_lib.h>
17 #include <deal.II/base/template_constraints.h>
18 #include <deal.II/base/types.h>
19 #include <deal.II/base/utilities.h>
20 
21 #include <deal.II/distributed/tria.h>
22 
23 #include <deal.II/dofs/dof_accessor.h>
24 #include <deal.II/dofs/dof_handler.h>
25 #include <deal.II/dofs/dof_renumbering.h>
26 #include <deal.II/dofs/dof_tools.h>
27 
28 #include <deal.II/fe/fe.h>
29 
30 #include <deal.II/grid/tria.h>
31 #include <deal.II/grid/tria_iterator.h>
32 
33 #include <deal.II/hp/dof_handler.h>
34 #include <deal.II/hp/fe_collection.h>
35 #include <deal.II/hp/fe_values.h>
36 
37 #include <deal.II/lac/affine_constraints.h>
38 #include <deal.II/lac/dynamic_sparsity_pattern.h>
39 #include <deal.II/lac/sparsity_pattern.h>
40 #include <deal.II/lac/sparsity_tools.h>
41 
42 #include <deal.II/multigrid/mg_tools.h>
43 
44 #include <boost/config.hpp>
45 #include <boost/graph/adjacency_list.hpp>
46 #include <boost/graph/bandwidth.hpp>
47 #include <boost/graph/cuthill_mckee_ordering.hpp>
48 #include <boost/graph/king_ordering.hpp>
49 #include <boost/graph/minimum_degree_ordering.hpp>
50 #include <boost/graph/properties.hpp>
51 #include <boost/random.hpp>
52 #include <boost/random/uniform_int_distribution.hpp>
53 
54 #include <algorithm>
55 #include <cmath>
56 #include <functional>
57 #include <map>
58 #include <vector>
59 
60 
61 DEAL_II_NAMESPACE_OPEN
62 
63 namespace DoFRenumbering
64 {
65  namespace boost
66  {
67  namespace boosttypes
68  {
69  using namespace ::boost;
70  using namespace std;
71 
72  using Graph = adjacency_list<vecS,
73  vecS,
74  undirectedS,
75  property<vertex_color_t,
76  default_color_type,
77  property<vertex_degree_t, int>>>;
78  using Vertex = graph_traits<Graph>::vertex_descriptor;
79  using size_type = graph_traits<Graph>::vertices_size_type;
80 
81  using Pair = std::pair<size_type, size_type>;
82  } // namespace boosttypes
83 
84 
85  namespace internal
86  {
87  template <typename DoFHandlerType>
88  void
89  create_graph(const DoFHandlerType &dof_handler,
90  const bool use_constraints,
91  boosttypes::Graph & graph,
92  boosttypes::property_map<boosttypes::Graph,
93  boosttypes::vertex_degree_t>::type
94  &graph_degree)
95  {
96  {
97  // create intermediate sparsity pattern (faster than directly
98  // submitting indices)
99  AffineConstraints<double> constraints;
100  if (use_constraints)
101  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
102  constraints.close();
103  DynamicSparsityPattern dsp(dof_handler.n_dofs(),
104  dof_handler.n_dofs());
105  DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
106 
107  // submit the entries to the boost graph
108  for (unsigned int row = 0; row < dsp.n_rows(); ++row)
109  for (unsigned int col = 0; col < dsp.row_length(row); ++col)
110  add_edge(row, dsp.column_number(row, col), graph);
111  }
112 
113  boosttypes::graph_traits<boosttypes::Graph>::vertex_iterator ui, ui_end;
114 
115  graph_degree = get(::boost::vertex_degree, graph);
116  for (::boost::tie(ui, ui_end) = vertices(graph); ui != ui_end; ++ui)
117  graph_degree[*ui] = degree(*ui, graph);
118  }
119  } // namespace internal
120 
121 
122  template <typename DoFHandlerType>
123  void
124  Cuthill_McKee(DoFHandlerType &dof_handler,
125  const bool reversed_numbering,
126  const bool use_constraints)
127  {
128  std::vector<types::global_dof_index> renumbering(
129  dof_handler.n_dofs(), numbers::invalid_dof_index);
130  compute_Cuthill_McKee(renumbering,
131  dof_handler,
132  reversed_numbering,
133  use_constraints);
134 
135  // actually perform renumbering;
136  // this is dimension specific and
137  // thus needs an own function
138  dof_handler.renumber_dofs(renumbering);
139  }
140 
141 
142  template <typename DoFHandlerType>
143  void
144  compute_Cuthill_McKee(std::vector<types::global_dof_index> &new_dof_indices,
145  const DoFHandlerType & dof_handler,
146  const bool reversed_numbering,
147  const bool use_constraints)
148  {
149  boosttypes::Graph graph(dof_handler.n_dofs());
150  boosttypes::property_map<boosttypes::Graph,
151  boosttypes::vertex_degree_t>::type graph_degree;
152 
153  internal::create_graph(dof_handler, use_constraints, graph, graph_degree);
154 
155  boosttypes::property_map<boosttypes::Graph,
156  boosttypes::vertex_index_t>::type index_map =
157  get(::boost::vertex_index, graph);
158 
159 
160  std::vector<boosttypes::Vertex> inv_perm(num_vertices(graph));
161 
162  if (reversed_numbering == false)
163  ::boost::cuthill_mckee_ordering(graph,
164  inv_perm.rbegin(),
165  get(::boost::vertex_color, graph),
166  make_degree_map(graph));
167  else
168  ::boost::cuthill_mckee_ordering(graph,
169  inv_perm.begin(),
170  get(::boost::vertex_color, graph),
171  make_degree_map(graph));
172 
173  for (boosttypes::size_type c = 0; c != inv_perm.size(); ++c)
174  new_dof_indices[index_map[inv_perm[c]]] = c;
175 
176  Assert(std::find(new_dof_indices.begin(),
177  new_dof_indices.end(),
178  numbers::invalid_dof_index) == new_dof_indices.end(),
179  ExcInternalError());
180  }
181 
182 
183 
184  template <typename DoFHandlerType>
185  void
186  king_ordering(DoFHandlerType &dof_handler,
187  const bool reversed_numbering,
188  const bool use_constraints)
189  {
190  std::vector<types::global_dof_index> renumbering(
191  dof_handler.n_dofs(), numbers::invalid_dof_index);
192  compute_king_ordering(renumbering,
193  dof_handler,
194  reversed_numbering,
195  use_constraints);
196 
197  // actually perform renumbering;
198  // this is dimension specific and
199  // thus needs an own function
200  dof_handler.renumber_dofs(renumbering);
201  }
202 
203 
204  template <typename DoFHandlerType>
205  void
206  compute_king_ordering(std::vector<types::global_dof_index> &new_dof_indices,
207  const DoFHandlerType & dof_handler,
208  const bool reversed_numbering,
209  const bool use_constraints)
210  {
211  boosttypes::Graph graph(dof_handler.n_dofs());
212  boosttypes::property_map<boosttypes::Graph,
213  boosttypes::vertex_degree_t>::type graph_degree;
214 
215  internal::create_graph(dof_handler, use_constraints, graph, graph_degree);
216 
217  boosttypes::property_map<boosttypes::Graph,
218  boosttypes::vertex_index_t>::type index_map =
219  get(::boost::vertex_index, graph);
220 
221 
222  std::vector<boosttypes::Vertex> inv_perm(num_vertices(graph));
223 
224  if (reversed_numbering == false)
225  ::boost::king_ordering(graph, inv_perm.rbegin());
226  else
227  ::boost::king_ordering(graph, inv_perm.begin());
228 
229  for (boosttypes::size_type c = 0; c != inv_perm.size(); ++c)
230  new_dof_indices[index_map[inv_perm[c]]] = c;
231 
232  Assert(std::find(new_dof_indices.begin(),
233  new_dof_indices.end(),
234  numbers::invalid_dof_index) == new_dof_indices.end(),
235  ExcInternalError());
236  }
237 
238 
239 
240  template <typename DoFHandlerType>
241  void
242  minimum_degree(DoFHandlerType &dof_handler,
243  const bool reversed_numbering,
244  const bool use_constraints)
245  {
246  std::vector<types::global_dof_index> renumbering(
247  dof_handler.n_dofs(), numbers::invalid_dof_index);
248  compute_minimum_degree(renumbering,
249  dof_handler,
250  reversed_numbering,
251  use_constraints);
252 
253  // actually perform renumbering;
254  // this is dimension specific and
255  // thus needs an own function
256  dof_handler.renumber_dofs(renumbering);
257  }
258 
259 
260  template <typename DoFHandlerType>
261  void
263  std::vector<types::global_dof_index> &new_dof_indices,
264  const DoFHandlerType & dof_handler,
265  const bool reversed_numbering,
266  const bool use_constraints)
267  {
268  (void)use_constraints;
269  Assert(use_constraints == false, ExcNotImplemented());
270 
271  // the following code is pretty
272  // much a verbatim copy of the
273  // sample code for the
274  // minimum_degree_ordering manual
275  // page from the BOOST Graph
276  // Library
277  using namespace ::boost;
278 
279  int delta = 0;
280 
281  // must be BGL directed graph now
282  using Graph = adjacency_list<vecS, vecS, directedS>;
283 
284  int n = dof_handler.n_dofs();
285 
286  Graph G(n);
287 
288  std::vector<::types::global_dof_index> dofs_on_this_cell;
289 
290  typename DoFHandlerType::active_cell_iterator cell = dof_handler
291  .begin_active(),
292  endc = dof_handler.end();
293 
294  for (; cell != endc; ++cell)
295  {
296  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
297 
298  dofs_on_this_cell.resize(dofs_per_cell);
299 
300  cell->get_active_or_mg_dof_indices(dofs_on_this_cell);
301  for (unsigned int i = 0; i < dofs_per_cell; ++i)
302  for (unsigned int j = 0; j < dofs_per_cell; ++j)
303  if (dofs_on_this_cell[i] > dofs_on_this_cell[j])
304  {
305  add_edge(dofs_on_this_cell[i], dofs_on_this_cell[j], G);
306  add_edge(dofs_on_this_cell[j], dofs_on_this_cell[i], G);
307  }
308  }
309 
310 
311  using Vector = std::vector<int>;
312 
313 
314  Vector inverse_perm(n, 0);
315 
316  Vector perm(n, 0);
317 
318 
319  Vector supernode_sizes(n, 1);
320  // init has to be 1
321 
322  ::boost::property_map<Graph, vertex_index_t>::type id =
323  get(vertex_index, G);
324 
325 
326  Vector degree(n, 0);
327 
328 
329  minimum_degree_ordering(
330  G,
331  make_iterator_property_map(degree.begin(), id, degree[0]),
332  inverse_perm.data(),
333  perm.data(),
334  make_iterator_property_map(supernode_sizes.begin(),
335  id,
336  supernode_sizes[0]),
337  delta,
338  id);
339 
340 
341  for (int i = 0; i < n; ++i)
342  {
343  Assert(std::find(perm.begin(), perm.end(), i) != perm.end(),
344  ExcInternalError());
345  Assert(std::find(inverse_perm.begin(), inverse_perm.end(), i) !=
346  inverse_perm.end(),
347  ExcInternalError());
348  Assert(inverse_perm[perm[i]] == i, ExcInternalError());
349  }
350 
351  if (reversed_numbering == true)
352  std::copy(perm.begin(), perm.end(), new_dof_indices.begin());
353  else
354  std::copy(inverse_perm.begin(),
355  inverse_perm.end(),
356  new_dof_indices.begin());
357  }
358 
359  } // namespace boost
360 
361 
362 
363  template <typename DoFHandlerType>
364  void
365  Cuthill_McKee(DoFHandlerType & dof_handler,
366  const bool reversed_numbering,
367  const bool use_constraints,
368  const std::vector<types::global_dof_index> &starting_indices)
369  {
370  std::vector<types::global_dof_index> renumbering(
371  dof_handler.locally_owned_dofs().n_elements(),
373  compute_Cuthill_McKee(renumbering,
374  dof_handler,
375  reversed_numbering,
376  use_constraints,
377  starting_indices);
378 
379  // actually perform renumbering;
380  // this is dimension specific and
381  // thus needs an own function
382  dof_handler.renumber_dofs(renumbering);
383  }
384 
385 
386 
387  template <typename DoFHandlerType>
388  void
390  std::vector<types::global_dof_index> & new_indices,
391  const DoFHandlerType & dof_handler,
392  const bool reversed_numbering,
393  const bool use_constraints,
394  const std::vector<types::global_dof_index> &starting_indices)
395  {
396  // see if there is anything to do at all or whether we can skip the work on
397  // this processor
398  if (dof_handler.locally_owned_dofs().n_elements() == 0)
399  {
400  Assert(new_indices.size() == 0, ExcInternalError());
401  return;
402  }
403 
404  // make the connection graph
405  //
406  // note that if constraints are not requested, then the 'constraints'
407  // object will be empty and using it has no effect
408  IndexSet locally_relevant_dofs;
409  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
410 
411  AffineConstraints<double> constraints;
412  if (use_constraints)
413  {
414  constraints.reinit(locally_relevant_dofs);
415  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
416  }
417  constraints.close();
418 
419  const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
420 
421  // see if we can get away with the sequential algorithm
422  if (locally_owned_dofs.n_elements() == locally_owned_dofs.size())
423  {
424  AssertDimension(new_indices.size(), dof_handler.n_dofs());
425 
426  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
427  DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
428 
430  new_indices,
431  starting_indices);
432  if (reversed_numbering)
433  new_indices = Utilities::reverse_permutation(new_indices);
434  }
435  else
436  {
437  // we are in the parallel case where we need to work in the
438  // local index space, i.e., the locally owned part of the
439  // sparsity pattern.
440  //
441  // first figure out whether the user only gave us starting
442  // indices that are locally owned, or that are only locally
443  // relevant. in the process, also check that all indices
444  // really belong to at least the locally relevant ones
445  IndexSet locally_active_dofs;
446  DoFTools::extract_locally_active_dofs(dof_handler, locally_active_dofs);
447 
448  bool needs_locally_active = false;
449  for (const auto starting_index : starting_indices)
450  {
451  if ((needs_locally_active ==
452  /* previously already set to */ true) ||
453  (locally_owned_dofs.is_element(starting_index) == false))
454  {
455  Assert(
456  locally_active_dofs.is_element(starting_index),
457  ExcMessage(
458  "You specified global degree of freedom " +
459  std::to_string(starting_index) +
460  " as a starting index, but this index is not among the "
461  "locally active ones on this processor, as required "
462  "for this function."));
463  needs_locally_active = true;
464  }
465  }
466 
467  const IndexSet index_set_to_use =
468  (needs_locally_active ? locally_active_dofs : locally_owned_dofs);
469 
470  // then create first the global sparsity pattern, and then the local
471  // sparsity pattern from the global one by transferring its indices to
472  // processor-local (locally owned or locally active) index space
473  DynamicSparsityPattern dsp(dof_handler.n_dofs(),
474  dof_handler.n_dofs(),
475  index_set_to_use);
476  DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
477 
478  DynamicSparsityPattern local_sparsity(index_set_to_use.n_elements(),
479  index_set_to_use.n_elements());
480  std::vector<types::global_dof_index> row_entries;
481  for (unsigned int i = 0; i < index_set_to_use.n_elements(); ++i)
482  {
483  const types::global_dof_index row =
484  index_set_to_use.nth_index_in_set(i);
485  const unsigned int row_length = dsp.row_length(row);
486  row_entries.clear();
487  for (unsigned int j = 0; j < row_length; ++j)
488  {
489  const unsigned int col = dsp.column_number(row, j);
490  if (col != row && index_set_to_use.is_element(col))
491  row_entries.push_back(index_set_to_use.index_within_set(col));
492  }
493  local_sparsity.add_entries(i,
494  row_entries.begin(),
495  row_entries.end(),
496  true);
497  }
498 
499  // translate starting indices from global to local indices
500  std::vector<types::global_dof_index> local_starting_indices(
501  starting_indices.size());
502  for (unsigned int i = 0; i < starting_indices.size(); ++i)
503  local_starting_indices[i] =
504  index_set_to_use.index_within_set(starting_indices[i]);
505 
506  // then do the renumbering on the locally owned portion
507  AssertDimension(new_indices.size(), locally_owned_dofs.n_elements());
508  std::vector<types::global_dof_index> my_new_indices(
509  index_set_to_use.n_elements());
511  my_new_indices,
512  local_starting_indices);
513  if (reversed_numbering)
514  my_new_indices = Utilities::reverse_permutation(my_new_indices);
515 
516  // now that we have a re-enumeration of all DoFs, we need to throw
517  // out the ones that are not locally owned in case we have worked
518  // with the locally active ones. that's because the renumbering
519  // functions only want new indices for the locally owned DoFs (other
520  // processors are responsible for renumbering the ones that are
521  // on cell interfaces)
522  if (needs_locally_active == true)
523  {
524  // first step: figure out which DoF indices to eliminate
525  IndexSet active_but_not_owned_dofs = locally_active_dofs;
526  active_but_not_owned_dofs.subtract_set(locally_owned_dofs);
527 
528  std::set<types::global_dof_index> erase_these_indices;
529  for (const auto p : active_but_not_owned_dofs)
530  {
531  const auto index = index_set_to_use.index_within_set(p);
532  Assert(index < index_set_to_use.n_elements(),
533  ExcInternalError());
534  erase_these_indices.insert(my_new_indices[index]);
535  my_new_indices[index] = numbers::invalid_dof_index;
536  }
537  Assert(erase_these_indices.size() ==
538  active_but_not_owned_dofs.n_elements(),
539  ExcInternalError());
540  Assert(static_cast<unsigned int>(
541  std::count(my_new_indices.begin(),
542  my_new_indices.end(),
544  active_but_not_owned_dofs.n_elements(),
545  ExcInternalError());
546 
547  // then compute a renumbering of the remaining ones
548  std::vector<types::global_dof_index> translate_indices(
549  my_new_indices.size());
550  {
551  std::set<types::global_dof_index>::const_iterator
552  next_erased_index = erase_these_indices.begin();
553  types::global_dof_index next_new_index = 0;
554  for (unsigned int i = 0; i < translate_indices.size(); ++i)
555  if ((next_erased_index != erase_these_indices.end()) &&
556  (*next_erased_index == i))
557  {
558  translate_indices[i] = numbers::invalid_dof_index;
559  ++next_erased_index;
560  }
561  else
562  {
563  translate_indices[i] = next_new_index;
564  ++next_new_index;
565  }
566  Assert(next_new_index == locally_owned_dofs.n_elements(),
567  ExcInternalError());
568  }
569 
570  // and then do the renumbering of the result of the
571  // Cuthill-McKee algorithm above, right into the output array
572  new_indices.clear();
573  new_indices.reserve(locally_owned_dofs.n_elements());
574  for (const auto &p : my_new_indices)
576  {
577  Assert(translate_indices[p] != numbers::invalid_dof_index,
578  ExcInternalError());
579  new_indices.push_back(translate_indices[p]);
580  }
581  Assert(new_indices.size() == locally_owned_dofs.n_elements(),
582  ExcInternalError());
583  }
584  else
585  new_indices = std::move(my_new_indices);
586 
587  // convert indices back to global index space. in both of the branches
588  // above, we ended up with new_indices only containing the local
589  // indices of the locally-owned DoFs. so that's where we get the
590  // indices
591  for (types::global_dof_index &new_index : new_indices)
592  new_index = locally_owned_dofs.nth_index_in_set(new_index);
593  }
594  }
595 
596 
597 
598  template <typename DoFHandlerType>
599  void
600  Cuthill_McKee(DoFHandlerType & dof_handler,
601  const unsigned int level,
602  const bool reversed_numbering,
603  const std::vector<types::global_dof_index> &starting_indices)
604  {
605  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
607 
608  // make the connection graph
609  DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
610  dof_handler.n_dofs(level));
611  MGTools::make_sparsity_pattern(dof_handler, dsp, level);
612 
613  std::vector<types::global_dof_index> new_indices(dsp.n_rows());
614  SparsityTools::reorder_Cuthill_McKee(dsp, new_indices, starting_indices);
615 
616  if (reversed_numbering)
617  new_indices = Utilities::reverse_permutation(new_indices);
618 
619  // actually perform renumbering;
620  // this is dimension specific and
621  // thus needs an own function
622  dof_handler.renumber_dofs(level, new_indices);
623  }
624 
625 
626 
627  template <typename DoFHandlerType>
628  void
629  component_wise(DoFHandlerType & dof_handler,
630  const std::vector<unsigned int> &component_order_arg)
631  {
632  std::vector<types::global_dof_index> renumbering(
633  dof_handler.n_locally_owned_dofs(), numbers::invalid_dof_index);
634 
635  typename DoFHandlerType::active_cell_iterator start =
636  dof_handler.begin_active();
637  const typename DoFHandlerType::level_cell_iterator end = dof_handler.end();
638 
639  const types::global_dof_index result =
640  compute_component_wise<DoFHandlerType::dimension,
641  DoFHandlerType::space_dimension>(
642  renumbering, start, end, component_order_arg, false);
643  if (result == 0)
644  return;
645 
646  // verify that the last numbered
647  // degree of freedom is either
648  // equal to the number of degrees
649  // of freedom in total (the
650  // sequential case) or in the
651  // distributed case at least
652  // makes sense
653  Assert((result == dof_handler.n_locally_owned_dofs()) ||
654  ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
655  (result <= dof_handler.n_dofs())),
656  ExcInternalError());
657 
658  dof_handler.renumber_dofs(renumbering);
659  }
660 
661 
662 
663  template <typename DoFHandlerType>
664  void
665  component_wise(DoFHandlerType & dof_handler,
666  const unsigned int level,
667  const std::vector<unsigned int> &component_order_arg)
668  {
669  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
671 
672  std::vector<types::global_dof_index> renumbering(
673  dof_handler.locally_owned_mg_dofs(level).n_elements(),
675 
676  typename DoFHandlerType::level_cell_iterator start =
677  dof_handler.begin(level);
678  typename DoFHandlerType::level_cell_iterator end = dof_handler.end(level);
679 
680  const types::global_dof_index result =
681  compute_component_wise<DoFHandlerType::dimension,
682  DoFHandlerType::space_dimension>(
683  renumbering, start, end, component_order_arg, true);
684 
685  if (result == 0)
686  return;
687 
688  Assert(result == dof_handler.n_dofs(level), ExcInternalError());
689 
690  if (renumbering.size() != 0)
691  dof_handler.renumber_dofs(level, renumbering);
692  }
693 
694 
695 
696  template <int dim, int spacedim, typename CellIterator>
698  compute_component_wise(std::vector<types::global_dof_index> &new_indices,
699  const CellIterator & start,
700  const typename identity<CellIterator>::type &end,
701  const std::vector<unsigned int> &component_order_arg,
702  const bool is_level_operation)
703  {
704  const hp::FECollection<dim, spacedim> fe_collection(
705  start->get_dof_handler().get_fe_collection());
706 
707  // do nothing if the FE has only
708  // one component
709  if (fe_collection.n_components() == 1)
710  {
711  new_indices.resize(0);
712  return 0;
713  }
714 
715  // Get a reference to the set of dofs. Note that we assume that all cells
716  // are assumed to be on the same level, otherwise the operation doesn't make
717  // much sense (we will assert this below).
718  const IndexSet &locally_owned_dofs =
719  is_level_operation ?
720  start->get_dof_handler().locally_owned_mg_dofs(start->level()) :
721  start->get_dof_handler().locally_owned_dofs();
722 
723  // Copy last argument into a
724  // writable vector.
725  std::vector<unsigned int> component_order(component_order_arg);
726  // If the last argument was an
727  // empty vector, set up things to
728  // store components in the order
729  // found in the system.
730  if (component_order.size() == 0)
731  for (unsigned int i = 0; i < fe_collection.n_components(); ++i)
732  component_order.push_back(i);
733 
734  Assert(component_order.size() == fe_collection.n_components(),
735  ExcDimensionMismatch(component_order.size(),
736  fe_collection.n_components()));
737 
738  for (const unsigned int component : component_order)
739  {
740  (void)component;
741  AssertIndexRange(component, fe_collection.n_components());
742  }
743 
744  // vector to hold the dof indices on
745  // the cell we visit at a time
746  std::vector<types::global_dof_index> local_dof_indices;
747 
748  // prebuilt list to which component
749  // a given dof on a cell
750  // should go. note that we get into
751  // trouble here if the shape
752  // function is not primitive, since
753  // then there is no single vector
754  // component to which it
755  // belongs. in this case, assign it
756  // to the first vector component to
757  // which it belongs
758  std::vector<std::vector<unsigned int>> component_list(fe_collection.size());
759  for (unsigned int f = 0; f < fe_collection.size(); ++f)
760  {
761  const FiniteElement<dim, spacedim> &fe = fe_collection[f];
762  const unsigned int dofs_per_cell = fe.dofs_per_cell;
763  component_list[f].resize(dofs_per_cell);
764  for (unsigned int i = 0; i < dofs_per_cell; ++i)
765  if (fe.is_primitive(i))
766  component_list[f][i] =
767  component_order[fe.system_to_component_index(i).first];
768  else
769  {
770  const unsigned int comp =
772 
773  // then associate this degree
774  // of freedom with this
775  // component
776  component_list[f][i] = component_order[comp];
777  }
778  }
779 
780  // set up a map where for each
781  // component the respective degrees
782  // of freedom are collected.
783  //
784  // note that this map is sorted by
785  // component but that within each
786  // component it is NOT sorted by
787  // dof index. note also that some
788  // dof indices are entered
789  // multiply, so we will have to
790  // take care of that
791  std::vector<std::vector<types::global_dof_index>> component_to_dof_map(
792  fe_collection.n_components());
793  for (CellIterator cell = start; cell != end; ++cell)
794  {
795  if (is_level_operation)
796  {
797  // we are dealing with mg dofs, skip foreign level cells:
798  if (!cell->is_locally_owned_on_level())
799  continue;
800  }
801  else
802  {
803  // we are dealing with active dofs, skip the loop if not locally
804  // owned:
805  if (!cell->is_locally_owned())
806  continue;
807  }
808 
809  if (is_level_operation)
810  Assert(
811  cell->level() == start->level(),
812  ExcMessage(
813  "Multigrid renumbering in compute_component_wise() needs to be applied to a single level!"));
814 
815  // on each cell: get dof indices
816  // and insert them into the global
817  // list using their component
818  const unsigned int fe_index = cell->active_fe_index();
819  const unsigned int dofs_per_cell =
820  fe_collection[fe_index].dofs_per_cell;
821  local_dof_indices.resize(dofs_per_cell);
822  cell->get_active_or_mg_dof_indices(local_dof_indices);
823 
824  for (unsigned int i = 0; i < dofs_per_cell; ++i)
825  if (locally_owned_dofs.is_element(local_dof_indices[i]))
826  component_to_dof_map[component_list[fe_index][i]].push_back(
827  local_dof_indices[i]);
828  }
829 
830  // now we've got all indices sorted
831  // into buckets labeled by their
832  // target component number. we've
833  // only got to traverse this list
834  // and assign the new indices
835  //
836  // however, we first want to sort
837  // the indices entered into the
838  // buckets to preserve the order
839  // within each component and during
840  // this also remove duplicate
841  // entries
842  //
843  // note that we no longer have to
844  // care about non-primitive shape
845  // functions since the buckets
846  // corresponding to the second and
847  // following vector components of a
848  // non-primitive FE will simply be
849  // empty, everything being shoved
850  // into the first one. The same
851  // holds if several components were
852  // joined into a single target.
853  for (unsigned int component = 0; component < fe_collection.n_components();
854  ++component)
855  {
856  std::sort(component_to_dof_map[component].begin(),
857  component_to_dof_map[component].end());
858  component_to_dof_map[component].erase(
859  std::unique(component_to_dof_map[component].begin(),
860  component_to_dof_map[component].end()),
861  component_to_dof_map[component].end());
862  }
863 
864  // calculate the number of locally owned
865  // DoFs per bucket
866  const unsigned int n_buckets = fe_collection.n_components();
867  std::vector<types::global_dof_index> shifts(n_buckets);
868 
870  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
871  &start->get_dof_handler().get_triangulation())))
872  {
873 #ifdef DEAL_II_WITH_MPI
874  std::vector<types::global_dof_index> local_dof_count(n_buckets);
875 
876  for (unsigned int c = 0; c < n_buckets; ++c)
877  local_dof_count[c] = component_to_dof_map[c].size();
878 
879  std::vector<types::global_dof_index> prefix_dof_count(n_buckets);
880  const int ierr = MPI_Exscan(local_dof_count.data(),
881  prefix_dof_count.data(),
882  n_buckets,
883  DEAL_II_DOF_INDEX_MPI_TYPE,
884  MPI_SUM,
885  tria->get_communicator());
886  AssertThrowMPI(ierr);
887 
888  std::vector<types::global_dof_index> global_dof_count(n_buckets);
889  Utilities::MPI::sum(local_dof_count,
890  tria->get_communicator(),
891  global_dof_count);
892 
893  // calculate shifts
894  types::global_dof_index cumulated = 0;
895  for (unsigned int c = 0; c < n_buckets; ++c)
896  {
897  shifts[c] = prefix_dof_count[c] + cumulated;
898  cumulated += global_dof_count[c];
899  }
900 #else
901  (void)tria;
902  Assert(false, ExcInternalError());
903 #endif
904  }
905  else
906  {
907  shifts[0] = 0;
908  for (unsigned int c = 1; c < fe_collection.n_components(); ++c)
909  shifts[c] = shifts[c - 1] + component_to_dof_map[c - 1].size();
910  }
911 
912 
913 
914  // now concatenate all the
915  // components in the order the user
916  // desired to see
917  types::global_dof_index next_free_index = 0;
918  for (unsigned int component = 0; component < fe_collection.n_components();
919  ++component)
920  {
921  const typename std::vector<types::global_dof_index>::const_iterator
922  begin_of_component = component_to_dof_map[component].begin(),
923  end_of_component = component_to_dof_map[component].end();
924 
925  next_free_index = shifts[component];
926 
927  for (typename std::vector<types::global_dof_index>::const_iterator
928  dof_index = begin_of_component;
929  dof_index != end_of_component;
930  ++dof_index)
931  {
932  Assert(locally_owned_dofs.index_within_set(*dof_index) <
933  new_indices.size(),
934  ExcInternalError());
935  new_indices[locally_owned_dofs.index_within_set(*dof_index)] =
936  next_free_index++;
937  }
938  }
939 
940  return next_free_index;
941  }
942 
943 
944 
945  template <int dim, int spacedim>
946  void
948  {
949  std::vector<types::global_dof_index> renumbering(
951 
953  dof_handler.begin_active();
954  const typename DoFHandler<dim, spacedim>::level_cell_iterator end =
955  dof_handler.end();
956 
958  dim,
959  spacedim,
961  typename DoFHandler<dim, spacedim>::level_cell_iterator>(renumbering,
962  start,
963  end,
964  false);
965  if (result == 0)
966  return;
967 
968  // verify that the last numbered
969  // degree of freedom is either
970  // equal to the number of degrees
971  // of freedom in total (the
972  // sequential case) or in the
973  // distributed case at least
974  // makes sense
975  Assert((result == dof_handler.n_locally_owned_dofs()) ||
976  ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
977  (result <= dof_handler.n_dofs())),
978  ExcInternalError());
979 
980  dof_handler.renumber_dofs(renumbering);
981  }
982 
983 
984 
985  template <int dim, int spacedim>
986  void
988  {
989  std::vector<types::global_dof_index> renumbering(
990  dof_handler.n_dofs(), numbers::invalid_dof_index);
991 
993  dof_handler.begin_active();
994  const typename hp::DoFHandler<dim, spacedim>::level_cell_iterator end =
995  dof_handler.end();
996 
998  dim,
999  spacedim,
1001  typename hp::DoFHandler<dim, spacedim>::level_cell_iterator>(renumbering,
1002  start,
1003  end,
1004  false);
1005 
1006  if (result == 0)
1007  return;
1008 
1009  Assert(result == dof_handler.n_dofs(), ExcInternalError());
1010 
1011  dof_handler.renumber_dofs(renumbering);
1012  }
1013 
1014 
1015 
1016  template <int dim, int spacedim>
1017  void
1018  block_wise(DoFHandler<dim, spacedim> &dof_handler, const unsigned int level)
1019  {
1020  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
1022 
1023  std::vector<types::global_dof_index> renumbering(
1024  dof_handler.locally_owned_mg_dofs(level).n_elements(),
1026 
1027  typename DoFHandler<dim, spacedim>::level_cell_iterator start =
1028  dof_handler.begin(level);
1029  typename DoFHandler<dim, spacedim>::level_cell_iterator end =
1030  dof_handler.end(level);
1031 
1033  dim,
1034  spacedim,
1035  typename DoFHandler<dim, spacedim>::level_cell_iterator,
1036  typename DoFHandler<dim, spacedim>::level_cell_iterator>(renumbering,
1037  start,
1038  end,
1039  true);
1040 
1041  if (result == 0)
1042  return;
1043 
1044  Assert(result == dof_handler.n_dofs(level), ExcInternalError());
1045 
1046  if (renumbering.size() != 0)
1047  dof_handler.renumber_dofs(level, renumbering);
1048  }
1049 
1050 
1051 
1052  template <int dim, int spacedim, class ITERATOR, class ENDITERATOR>
1054  compute_block_wise(std::vector<types::global_dof_index> &new_indices,
1055  const ITERATOR & start,
1056  const ENDITERATOR & end,
1057  const bool is_level_operation)
1058  {
1059  const hp::FECollection<dim, spacedim> fe_collection(
1060  start->get_dof_handler().get_fe_collection());
1061 
1062  // do nothing if the FE has only
1063  // one component
1064  if (fe_collection.n_blocks() == 1)
1065  {
1066  new_indices.resize(0);
1067  return 0;
1068  }
1069 
1070  // Get a reference to the set of dofs. Note that we assume that all cells
1071  // are assumed to be on the same level, otherwise the operation doesn't make
1072  // much sense (we will assert this below).
1073  const IndexSet &locally_owned_dofs =
1074  is_level_operation ?
1075  start->get_dof_handler().locally_owned_mg_dofs(start->level()) :
1076  start->get_dof_handler().locally_owned_dofs();
1077 
1078  // vector to hold the dof indices on
1079  // the cell we visit at a time
1080  std::vector<types::global_dof_index> local_dof_indices;
1081 
1082  // prebuilt list to which block
1083  // a given dof on a cell
1084  // should go.
1085  std::vector<std::vector<types::global_dof_index>> block_list(
1086  fe_collection.size());
1087  for (unsigned int f = 0; f < fe_collection.size(); ++f)
1088  {
1089  const FiniteElement<dim, spacedim> &fe = fe_collection[f];
1090  block_list[f].resize(fe.dofs_per_cell);
1091  for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
1092  block_list[f][i] = fe.system_to_block_index(i).first;
1093  }
1094 
1095  // set up a map where for each
1096  // block the respective degrees
1097  // of freedom are collected.
1098  //
1099  // note that this map is sorted by
1100  // block but that within each
1101  // block it is NOT sorted by
1102  // dof index. note also that some
1103  // dof indices are entered
1104  // multiply, so we will have to
1105  // take care of that
1106  std::vector<std::vector<types::global_dof_index>> block_to_dof_map(
1107  fe_collection.n_blocks());
1108  for (ITERATOR cell = start; cell != end; ++cell)
1109  {
1110  if (is_level_operation)
1111  {
1112  // we are dealing with mg dofs, skip foreign level cells:
1113  if (!cell->is_locally_owned_on_level())
1114  continue;
1115  }
1116  else
1117  {
1118  // we are dealing with active dofs, skip the loop if not locally
1119  // owned:
1120  if (!cell->is_locally_owned())
1121  continue;
1122  }
1123 
1124  if (is_level_operation)
1125  Assert(
1126  cell->level() == start->level(),
1127  ExcMessage(
1128  "Multigrid renumbering in compute_block_wise() needs to be applied to a single level!"));
1129 
1130  // on each cell: get dof indices
1131  // and insert them into the global
1132  // list using their component
1133  const unsigned int fe_index = cell->active_fe_index();
1134  const unsigned int dofs_per_cell =
1135  fe_collection[fe_index].dofs_per_cell;
1136  local_dof_indices.resize(dofs_per_cell);
1137  cell->get_active_or_mg_dof_indices(local_dof_indices);
1138 
1139  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1140  if (locally_owned_dofs.is_element(local_dof_indices[i]))
1141  block_to_dof_map[block_list[fe_index][i]].push_back(
1142  local_dof_indices[i]);
1143  }
1144 
1145  // now we've got all indices sorted
1146  // into buckets labeled by their
1147  // target block number. we've
1148  // only got to traverse this list
1149  // and assign the new indices
1150  //
1151  // however, we first want to sort
1152  // the indices entered into the
1153  // buckets to preserve the order
1154  // within each component and during
1155  // this also remove duplicate
1156  // entries
1157  for (unsigned int block = 0; block < fe_collection.n_blocks(); ++block)
1158  {
1159  std::sort(block_to_dof_map[block].begin(),
1160  block_to_dof_map[block].end());
1161  block_to_dof_map[block].erase(
1162  std::unique(block_to_dof_map[block].begin(),
1163  block_to_dof_map[block].end()),
1164  block_to_dof_map[block].end());
1165  }
1166 
1167  // calculate the number of locally owned
1168  // DoFs per bucket
1169  const unsigned int n_buckets = fe_collection.n_blocks();
1170  std::vector<types::global_dof_index> shifts(n_buckets);
1171 
1173  (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1174  &start->get_dof_handler().get_triangulation())))
1175  {
1176 #ifdef DEAL_II_WITH_MPI
1177  std::vector<types::global_dof_index> local_dof_count(n_buckets);
1178 
1179  for (unsigned int c = 0; c < n_buckets; ++c)
1180  local_dof_count[c] = block_to_dof_map[c].size();
1181 
1182  std::vector<types::global_dof_index> prefix_dof_count(n_buckets);
1183  const int ierr = MPI_Exscan(local_dof_count.data(),
1184  prefix_dof_count.data(),
1185  n_buckets,
1186  DEAL_II_DOF_INDEX_MPI_TYPE,
1187  MPI_SUM,
1188  tria->get_communicator());
1189  AssertThrowMPI(ierr);
1190 
1191  std::vector<types::global_dof_index> global_dof_count(n_buckets);
1192  Utilities::MPI::sum(local_dof_count,
1193  tria->get_communicator(),
1194  global_dof_count);
1195 
1196  // calculate shifts
1197  types::global_dof_index cumulated = 0;
1198  for (unsigned int c = 0; c < n_buckets; ++c)
1199  {
1200  shifts[c] = prefix_dof_count[c] + cumulated;
1201  cumulated += global_dof_count[c];
1202  }
1203 #else
1204  (void)tria;
1205  Assert(false, ExcInternalError());
1206 #endif
1207  }
1208  else
1209  {
1210  shifts[0] = 0;
1211  for (unsigned int c = 1; c < fe_collection.n_blocks(); ++c)
1212  shifts[c] = shifts[c - 1] + block_to_dof_map[c - 1].size();
1213  }
1214 
1215 
1216 
1217  // now concatenate all the
1218  // components in the order the user
1219  // desired to see
1220  types::global_dof_index next_free_index = 0;
1221  for (unsigned int block = 0; block < fe_collection.n_blocks(); ++block)
1222  {
1223  const typename std::vector<types::global_dof_index>::const_iterator
1224  begin_of_component = block_to_dof_map[block].begin(),
1225  end_of_component = block_to_dof_map[block].end();
1226 
1227  next_free_index = shifts[block];
1228 
1229  for (typename std::vector<types::global_dof_index>::const_iterator
1230  dof_index = begin_of_component;
1231  dof_index != end_of_component;
1232  ++dof_index)
1233  {
1234  Assert(locally_owned_dofs.index_within_set(*dof_index) <
1235  new_indices.size(),
1236  ExcInternalError());
1237  new_indices[locally_owned_dofs.index_within_set(*dof_index)] =
1238  next_free_index++;
1239  }
1240  }
1241 
1242  return next_free_index;
1243  }
1244 
1245 
1246 
1247  namespace
1248  {
1249  // Helper function for DoFRenumbering::hierarchical(). This function
1250  // recurses into the given cell or, if that should be an active (terminal)
1251  // cell, renumbers DoF indices on it. The function starts renumbering with
1252  // 'next_free_dof_index' and returns the first still unused DoF index at the
1253  // end of its operation.
1254  template <int dim, class CellIteratorType>
1256  compute_hierarchical_recursive(
1257  const types::global_dof_index next_free_dof_offset,
1258  const types::global_dof_index my_starting_index,
1259  const CellIteratorType & cell,
1260  const IndexSet & locally_owned_dof_indices,
1261  std::vector<types::global_dof_index> &new_indices)
1262  {
1263  types::global_dof_index current_next_free_dof_offset =
1264  next_free_dof_offset;
1265 
1266  if (cell->has_children())
1267  {
1268  // recursion
1269  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1270  ++c)
1271  current_next_free_dof_offset =
1272  compute_hierarchical_recursive<dim>(current_next_free_dof_offset,
1273  my_starting_index,
1274  cell->child(c),
1275  locally_owned_dof_indices,
1276  new_indices);
1277  }
1278  else
1279  {
1280  // this is a terminal cell. we need to renumber its DoF indices. there
1281  // are now three cases to decide:
1282  // - this is a sequential triangulation: we can just go ahead and
1283  // number
1284  // the DoFs in the order in which we encounter cells. in this case,
1285  // all cells are actually locally owned
1286  // - if this is a parallel::distributed::Triangulation, then we only
1287  // need to work on the locally owned cells since they contain
1288  // all locally owned DoFs.
1289  // - if this is a parallel::shared::Triangulation, then the same
1290  // applies
1291  //
1292  // in all cases, each processor starts new indices so that we get
1293  // a consecutive numbering on each processor, and disjoint ownership
1294  // of the global range of DoF indices
1295  if (cell->is_locally_owned())
1296  {
1297  // first get the existing DoF indices
1298  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
1299  std::vector<types::global_dof_index> local_dof_indices(
1300  dofs_per_cell);
1301  cell->get_dof_indices(local_dof_indices);
1302 
1303  // then loop over the existing DoF indices on this cell
1304  // and see whether it has already been re-numbered (it
1305  // may have been on a face or vertex to a neighboring
1306  // cell that we have encountered before). if not,
1307  // give it a new number and store that number in the
1308  // output array (new_indices)
1309  //
1310  // if this is a parallel triangulation and a DoF index is
1311  // not locally owned, then don't touch it. since
1312  // we don't actually *change* DoF indices (just record new
1313  // numbers in an array), we don't need to worry about
1314  // the decision whether a DoF is locally owned or not changing
1315  // as we progress in renumbering DoFs -- all adjacent cells
1316  // will always agree that a DoF is locally owned or not.
1317  // that said, the first cell to encounter a locally owned DoF
1318  // gets to number it, so the order in which we traverse cells
1319  // matters
1320  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1321  if (locally_owned_dof_indices.is_element(local_dof_indices[i]))
1322  {
1323  // this is a locally owned DoF, assign new number if not
1324  // assigned a number yet
1325  const unsigned int idx =
1326  locally_owned_dof_indices.index_within_set(
1327  local_dof_indices[i]);
1328  if (new_indices[idx] == numbers::invalid_dof_index)
1329  {
1330  new_indices[idx] =
1331  my_starting_index + current_next_free_dof_offset;
1332  ++current_next_free_dof_offset;
1333  }
1334  }
1335  }
1336  }
1337 
1338  return current_next_free_dof_offset;
1339  }
1340  } // namespace
1341 
1342 
1343 
1344  template <typename DoFHandlerType>
1345  void
1346  hierarchical(DoFHandlerType &dof_handler)
1347  {
1348  const int dim = DoFHandlerType::dimension;
1349  const int spacedim = DoFHandlerType::space_dimension;
1350 
1351  std::vector<types::global_dof_index> renumbering(
1352  dof_handler.n_locally_owned_dofs(), numbers::invalid_dof_index);
1353 
1354  types::global_dof_index next_free_dof_offset = 0;
1355  const IndexSet locally_owned = dof_handler.locally_owned_dofs();
1356 
1357  // in the function we call recursively, we want to number DoFs so
1358  // that global cell zero starts with DoF zero, regardless of how
1359  // DoFs were previously numbered. to this end, we need to figure
1360  // out which DoF index the current processor should start with.
1361  //
1362  // if this is a sequential triangulation, then obviously the starting
1363  // index is zero. otherwise, make sure we get contiguous, successive
1364  // ranges on each processor. note that since the number of locally owned
1365  // DoFs is an invariant under renumbering, we can easily compute this
1366  // starting index by just accumulating over the number of locally owned
1367  // DoFs for all previous processes
1368  types::global_dof_index my_starting_index = 0;
1369 
1371  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1372  &dof_handler.get_triangulation()))
1373  {
1374 #ifdef DEAL_II_WITH_MPI
1375  types::global_dof_index local_size =
1376  dof_handler.locally_owned_dofs().n_elements();
1377  MPI_Exscan(&local_size,
1378  &my_starting_index,
1379  1,
1380  DEAL_II_DOF_INDEX_MPI_TYPE,
1381  MPI_SUM,
1382  tria->get_communicator());
1383 #endif
1384  }
1385 
1388  *>(&dof_handler.get_triangulation()))
1389  {
1390 #ifdef DEAL_II_WITH_P4EST
1391  // this is a distributed Triangulation. we need to traverse the coarse
1392  // cells in the order p4est does to match the z-order actually used
1393  // by p4est. this requires using the renumbering of coarse cells
1394  // we do before we hand things off to p4est
1395  for (unsigned int c = 0; c < tria->n_cells(0); ++c)
1396  {
1397  const unsigned int coarse_cell_index =
1398  tria->get_p4est_tree_to_coarse_cell_permutation()[c];
1399 
1400  const typename DoFHandlerType::level_cell_iterator this_cell(
1401  tria, 0, coarse_cell_index, &dof_handler);
1402 
1403  next_free_dof_offset =
1404  compute_hierarchical_recursive<dim>(next_free_dof_offset,
1405  my_starting_index,
1406  this_cell,
1407  locally_owned,
1408  renumbering);
1409  }
1410 #else
1411  Assert(false, ExcNotImplemented());
1412 #endif
1413  }
1414  else
1415  {
1416  // this is not a distributed Triangulation, so we can traverse coarse
1417  // cells in the normal order
1418  for (typename DoFHandlerType::cell_iterator cell = dof_handler.begin(0);
1419  cell != dof_handler.end(0);
1420  ++cell)
1421  next_free_dof_offset =
1422  compute_hierarchical_recursive<dim>(next_free_dof_offset,
1423  my_starting_index,
1424  cell,
1425  locally_owned,
1426  renumbering);
1427  }
1428 
1429  // verify that the last numbered degree of freedom is either
1430  // equal to the number of degrees of freedom in total (the
1431  // sequential case) or in the distributed case at least
1432  // makes sense
1433  Assert((next_free_dof_offset == dof_handler.n_locally_owned_dofs()) ||
1434  ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
1435  (next_free_dof_offset <= dof_handler.n_dofs())),
1436  ExcInternalError());
1437 
1438  // make sure that all local DoFs got new numbers assigned
1439  Assert(std::find(renumbering.begin(),
1440  renumbering.end(),
1441  numbers::invalid_dof_index) == renumbering.end(),
1442  ExcInternalError());
1443 
1444  dof_handler.renumber_dofs(renumbering);
1445  }
1446 
1447 
1448 
1449  template <typename DoFHandlerType>
1450  void
1451  sort_selected_dofs_back(DoFHandlerType & dof_handler,
1452  const std::vector<bool> &selected_dofs)
1453  {
1454  std::vector<types::global_dof_index> renumbering(
1455  dof_handler.n_dofs(), numbers::invalid_dof_index);
1456  compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs);
1457 
1458  dof_handler.renumber_dofs(renumbering);
1459  }
1460 
1461 
1462 
1463  template <typename DoFHandlerType>
1464  void
1465  sort_selected_dofs_back(DoFHandlerType & dof_handler,
1466  const std::vector<bool> &selected_dofs,
1467  const unsigned int level)
1468  {
1469  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
1471 
1472  std::vector<types::global_dof_index> renumbering(
1473  dof_handler.n_dofs(level), numbers::invalid_dof_index);
1474  compute_sort_selected_dofs_back(renumbering,
1475  dof_handler,
1476  selected_dofs,
1477  level);
1478 
1479  dof_handler.renumber_dofs(level, renumbering);
1480  }
1481 
1482 
1483 
1484  template <typename DoFHandlerType>
1485  void
1487  std::vector<types::global_dof_index> &new_indices,
1488  const DoFHandlerType & dof_handler,
1489  const std::vector<bool> & selected_dofs)
1490  {
1491  const types::global_dof_index n_dofs = dof_handler.n_dofs();
1492  Assert(selected_dofs.size() == n_dofs,
1493  ExcDimensionMismatch(selected_dofs.size(), n_dofs));
1494 
1495  // re-sort the dofs according to
1496  // their selection state
1497  Assert(new_indices.size() == n_dofs,
1498  ExcDimensionMismatch(new_indices.size(), n_dofs));
1499 
1500  const types::global_dof_index n_selected_dofs =
1501  std::count(selected_dofs.begin(), selected_dofs.end(), false);
1502 
1503  types::global_dof_index next_unselected = 0;
1504  types::global_dof_index next_selected = n_selected_dofs;
1505  for (types::global_dof_index i = 0; i < n_dofs; ++i)
1506  if (selected_dofs[i] == false)
1507  {
1508  new_indices[i] = next_unselected;
1509  ++next_unselected;
1510  }
1511  else
1512  {
1513  new_indices[i] = next_selected;
1514  ++next_selected;
1515  }
1516  Assert(next_unselected == n_selected_dofs, ExcInternalError());
1517  Assert(next_selected == n_dofs, ExcInternalError());
1518  }
1519 
1520 
1521 
1522  template <typename DoFHandlerType>
1523  void
1525  std::vector<types::global_dof_index> &new_indices,
1526  const DoFHandlerType & dof_handler,
1527  const std::vector<bool> & selected_dofs,
1528  const unsigned int level)
1529  {
1530  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
1532 
1533  const unsigned int n_dofs = dof_handler.n_dofs(level);
1534  Assert(selected_dofs.size() == n_dofs,
1535  ExcDimensionMismatch(selected_dofs.size(), n_dofs));
1536 
1537  // re-sort the dofs according to
1538  // their selection state
1539  Assert(new_indices.size() == n_dofs,
1540  ExcDimensionMismatch(new_indices.size(), n_dofs));
1541 
1542  const unsigned int n_selected_dofs =
1543  std::count(selected_dofs.begin(), selected_dofs.end(), false);
1544 
1545  unsigned int next_unselected = 0;
1546  unsigned int next_selected = n_selected_dofs;
1547  for (unsigned int i = 0; i < n_dofs; ++i)
1548  if (selected_dofs[i] == false)
1549  {
1550  new_indices[i] = next_unselected;
1551  ++next_unselected;
1552  }
1553  else
1554  {
1555  new_indices[i] = next_selected;
1556  ++next_selected;
1557  }
1558  Assert(next_unselected == n_selected_dofs, ExcInternalError());
1559  Assert(next_selected == n_dofs, ExcInternalError());
1560  }
1561 
1562 
1563 
1564  template <typename DoFHandlerType>
1565  void
1567  DoFHandlerType & dof,
1568  const std::vector<typename DoFHandlerType::active_cell_iterator> &cells)
1569  {
1570  std::vector<types::global_dof_index> renumbering(
1571  dof.n_locally_owned_dofs());
1572  std::vector<types::global_dof_index> reverse(dof.n_locally_owned_dofs());
1573  compute_cell_wise(renumbering, reverse, dof, cells);
1574 
1575  dof.renumber_dofs(renumbering);
1576  }
1577 
1578 
1579  template <typename DoFHandlerType>
1580  void
1582  std::vector<types::global_dof_index> &new_indices,
1583  std::vector<types::global_dof_index> &reverse,
1584  const DoFHandlerType & dof,
1585  const typename std::vector<typename DoFHandlerType::active_cell_iterator>
1586  &cells)
1587  {
1588  if (const parallel::TriangulationBase<DoFHandlerType::dimension,
1589  DoFHandlerType::space_dimension> *p =
1590  dynamic_cast<const parallel::TriangulationBase<
1591  DoFHandlerType::dimension,
1592  DoFHandlerType::space_dimension> *>(&dof.get_triangulation()))
1593  {
1594  AssertDimension(cells.size(), p->n_locally_owned_active_cells());
1595  }
1596  else
1597  {
1598  AssertDimension(cells.size(), dof.get_triangulation().n_active_cells());
1599  }
1600 
1601  const auto n_owned_dofs = dof.n_locally_owned_dofs();
1602 
1603  // Actually, we compute the inverse of the reordering vector, called reverse
1604  // here. Later, its inverse is computed into new_indices, which is the
1605  // return argument.
1606 
1607  AssertDimension(new_indices.size(), n_owned_dofs);
1608  AssertDimension(reverse.size(), n_owned_dofs);
1609 
1610  // For continuous elements, we must make sure, that each dof is reordered
1611  // only once.
1612  std::vector<bool> already_sorted(n_owned_dofs, false);
1613  std::vector<types::global_dof_index> cell_dofs;
1614 
1615  const auto &owned_dofs = dof.locally_owned_dofs();
1616 
1617  unsigned int index = 0;
1618 
1619  for (const auto &cell : cells)
1620  {
1621  // Determine the number of dofs on this cell and reinit the
1622  // vector storing these numbers.
1623  const unsigned int n_cell_dofs = cell->get_fe().n_dofs_per_cell();
1624  cell_dofs.resize(n_cell_dofs);
1625 
1626  cell->get_active_or_mg_dof_indices(cell_dofs);
1627 
1628  // Sort here to make sure that degrees of freedom inside a single cell
1629  // are in the same order after renumbering.
1630  std::sort(cell_dofs.begin(), cell_dofs.end());
1631 
1632  for (const auto dof : cell_dofs)
1633  {
1634  const auto local_dof = owned_dofs.index_within_set(dof);
1635  if (local_dof != numbers::invalid_dof_index &&
1636  !already_sorted[local_dof])
1637  {
1638  already_sorted[local_dof] = true;
1639  reverse[index++] = local_dof;
1640  }
1641  }
1642  }
1643  Assert(index == n_owned_dofs,
1644  ExcMessage(
1645  "Traversing over the given set of cells did not cover all "
1646  "degrees of freedom in the DoFHandler. Does the set of cells "
1647  "not include all active cells?"));
1648 
1649  for (types::global_dof_index i = 0; i < reverse.size(); ++i)
1650  new_indices[reverse[i]] = owned_dofs.nth_index_in_set(i);
1651  }
1652 
1653 
1654 
1655  template <typename DoFHandlerType>
1656  void
1657  cell_wise(
1658  DoFHandlerType & dof,
1659  const unsigned int level,
1660  const typename std::vector<typename DoFHandlerType::level_cell_iterator>
1661  &cells)
1662  {
1663  Assert(dof.n_dofs(level) != numbers::invalid_dof_index,
1665 
1666  std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1667  std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1668 
1669  compute_cell_wise(renumbering, reverse, dof, level, cells);
1670  dof.renumber_dofs(level, renumbering);
1671  }
1672 
1673 
1674 
1675  template <typename DoFHandlerType>
1676  void
1678  std::vector<types::global_dof_index> &new_order,
1679  std::vector<types::global_dof_index> &reverse,
1680  const DoFHandlerType & dof,
1681  const unsigned int level,
1682  const typename std::vector<typename DoFHandlerType::level_cell_iterator>
1683  &cells)
1684  {
1685  Assert(cells.size() == dof.get_triangulation().n_cells(level),
1686  ExcDimensionMismatch(cells.size(),
1687  dof.get_triangulation().n_cells(level)));
1688  Assert(new_order.size() == dof.n_dofs(level),
1689  ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
1690  Assert(reverse.size() == dof.n_dofs(level),
1691  ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
1692 
1693  unsigned int n_global_dofs = dof.n_dofs(level);
1694  unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
1695 
1696  std::vector<bool> already_sorted(n_global_dofs, false);
1697  std::vector<types::global_dof_index> cell_dofs(n_cell_dofs);
1698 
1699  unsigned int global_index = 0;
1700 
1701  for (const auto &cell : cells)
1702  {
1703  Assert(cell->level() == static_cast<int>(level), ExcInternalError());
1704 
1705  cell->get_active_or_mg_dof_indices(cell_dofs);
1706  std::sort(cell_dofs.begin(), cell_dofs.end());
1707 
1708  for (unsigned int i = 0; i < n_cell_dofs; ++i)
1709  {
1710  if (!already_sorted[cell_dofs[i]])
1711  {
1712  already_sorted[cell_dofs[i]] = true;
1713  reverse[global_index++] = cell_dofs[i];
1714  }
1715  }
1716  }
1717  Assert(global_index == n_global_dofs,
1718  ExcMessage(
1719  "Traversing over the given set of cells did not cover all "
1720  "degrees of freedom in the DoFHandler. Does the set of cells "
1721  "not include all cells of the specified level?"));
1722 
1723  for (types::global_dof_index i = 0; i < new_order.size(); ++i)
1724  new_order[reverse[i]] = i;
1725  }
1726 
1727 
1728 
1729  template <typename DoFHandlerType>
1730  void
1731  downstream(DoFHandlerType & dof,
1733  const bool dof_wise_renumbering)
1734  {
1735  std::vector<types::global_dof_index> renumbering(dof.n_dofs());
1736  std::vector<types::global_dof_index> reverse(dof.n_dofs());
1738  renumbering, reverse, dof, direction, dof_wise_renumbering);
1739 
1740  dof.renumber_dofs(renumbering);
1741  }
1742 
1743 
1744 
1745  template <typename DoFHandlerType>
1746  void
1748  std::vector<types::global_dof_index> & new_indices,
1749  std::vector<types::global_dof_index> & reverse,
1750  const DoFHandlerType & dof,
1752  const bool dof_wise_renumbering)
1753  {
1754  Assert(
1755  (dynamic_cast<
1756  const parallel::TriangulationBase<DoFHandlerType::dimension,
1757  DoFHandlerType::space_dimension> *>(
1758  &dof.get_triangulation()) == nullptr),
1759  ExcNotImplemented());
1760 
1761  if (dof_wise_renumbering == false)
1762  {
1763  std::vector<typename DoFHandlerType::active_cell_iterator>
1764  ordered_cells;
1765  ordered_cells.reserve(dof.get_triangulation().n_active_cells());
1766  const CompareDownstream<typename DoFHandlerType::active_cell_iterator,
1767  DoFHandlerType::space_dimension>
1768  comparator(direction);
1769 
1770  typename DoFHandlerType::active_cell_iterator p = dof.begin_active();
1771  typename DoFHandlerType::active_cell_iterator end = dof.end();
1772 
1773  while (p != end)
1774  {
1775  ordered_cells.push_back(p);
1776  ++p;
1777  }
1778  std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
1779 
1780  compute_cell_wise(new_indices, reverse, dof, ordered_cells);
1781  }
1782  else
1783  {
1784  // similar code as for
1785  // DoFTools::map_dofs_to_support_points, but
1786  // need to do this for general DoFHandlerType classes and
1787  // want to be able to sort the result
1788  // (otherwise, could use something like
1789  // DoFTools::map_support_points_to_dofs)
1790  const unsigned int n_dofs = dof.n_dofs();
1791  std::vector<
1792  std::pair<Point<DoFHandlerType::space_dimension>, unsigned int>>
1793  support_point_list(n_dofs);
1794 
1795  const hp::FECollection<DoFHandlerType::dimension> &fe_collection =
1796  dof.get_fe_collection();
1797  Assert(fe_collection[0].has_support_points(),
1798  typename FiniteElement<
1799  DoFHandlerType::dimension>::ExcFEHasNoSupportPoints());
1800  hp::QCollection<DoFHandlerType::dimension> quadrature_collection;
1801  for (unsigned int comp = 0; comp < fe_collection.size(); ++comp)
1802  {
1803  Assert(fe_collection[comp].has_support_points(),
1804  typename FiniteElement<
1805  DoFHandlerType::dimension>::ExcFEHasNoSupportPoints());
1806  quadrature_collection.push_back(
1808  fe_collection[comp].get_unit_support_points()));
1809  }
1811  hp_fe_values(fe_collection,
1812  quadrature_collection,
1814 
1815  std::vector<bool> already_touched(n_dofs, false);
1816 
1817  std::vector<types::global_dof_index> local_dof_indices;
1818  typename DoFHandlerType::active_cell_iterator begin =
1819  dof.begin_active();
1820  typename DoFHandlerType::active_cell_iterator end = dof.end();
1821  for (; begin != end; ++begin)
1822  {
1823  const unsigned int dofs_per_cell = begin->get_fe().dofs_per_cell;
1824  local_dof_indices.resize(dofs_per_cell);
1825  hp_fe_values.reinit(begin);
1826  const FEValues<DoFHandlerType::dimension> &fe_values =
1827  hp_fe_values.get_present_fe_values();
1828  begin->get_active_or_mg_dof_indices(local_dof_indices);
1829  const std::vector<Point<DoFHandlerType::space_dimension>> &points =
1830  fe_values.get_quadrature_points();
1831  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1832  if (!already_touched[local_dof_indices[i]])
1833  {
1834  support_point_list[local_dof_indices[i]].first = points[i];
1835  support_point_list[local_dof_indices[i]].second =
1836  local_dof_indices[i];
1837  already_touched[local_dof_indices[i]] = true;
1838  }
1839  }
1840 
1842  direction);
1843  std::sort(support_point_list.begin(),
1844  support_point_list.end(),
1845  comparator);
1846  for (types::global_dof_index i = 0; i < n_dofs; ++i)
1847  new_indices[support_point_list[i].second] = i;
1848  }
1849  }
1850 
1851 
1852 
1853  template <typename DoFHandlerType>
1854  void
1855  downstream(DoFHandlerType & dof,
1856  const unsigned int level,
1858  const bool dof_wise_renumbering)
1859  {
1860  std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1861  std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1863  renumbering, reverse, dof, level, direction, dof_wise_renumbering);
1864 
1865  dof.renumber_dofs(level, renumbering);
1866  }
1867 
1868 
1869 
1870  template <typename DoFHandlerType>
1871  void
1873  std::vector<types::global_dof_index> & new_indices,
1874  std::vector<types::global_dof_index> & reverse,
1875  const DoFHandlerType & dof,
1876  const unsigned int level,
1878  const bool dof_wise_renumbering)
1879  {
1880  if (dof_wise_renumbering == false)
1881  {
1882  std::vector<typename DoFHandlerType::level_cell_iterator> ordered_cells;
1883  ordered_cells.reserve(dof.get_triangulation().n_cells(level));
1884  const CompareDownstream<typename DoFHandlerType::level_cell_iterator,
1885  DoFHandlerType::space_dimension>
1886  comparator(direction);
1887 
1888  typename DoFHandlerType::level_cell_iterator p = dof.begin(level);
1889  typename DoFHandlerType::level_cell_iterator end = dof.end(level);
1890 
1891  while (p != end)
1892  {
1893  ordered_cells.push_back(p);
1894  ++p;
1895  }
1896  std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
1897 
1898  compute_cell_wise(new_indices, reverse, dof, level, ordered_cells);
1899  }
1900  else
1901  {
1902  Assert(dof.get_fe().has_support_points(),
1903  typename FiniteElement<
1904  DoFHandlerType::dimension>::ExcFEHasNoSupportPoints());
1905  const unsigned int n_dofs = dof.n_dofs(level);
1906  std::vector<
1907  std::pair<Point<DoFHandlerType::space_dimension>, unsigned int>>
1908  support_point_list(n_dofs);
1909 
1911  dof.get_fe().get_unit_support_points());
1913  fe_values(dof.get_fe(), q_dummy, update_quadrature_points);
1914 
1915  std::vector<bool> already_touched(dof.n_dofs(), false);
1916 
1917  const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
1918  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1919  typename DoFHandlerType::level_cell_iterator begin = dof.begin(level);
1920  typename DoFHandlerType::level_cell_iterator end = dof.end(level);
1921  for (; begin != end; ++begin)
1922  {
1923  const typename Triangulation<
1924  DoFHandlerType::dimension,
1925  DoFHandlerType::space_dimension>::cell_iterator &begin_tria =
1926  begin;
1927  begin->get_active_or_mg_dof_indices(local_dof_indices);
1928  fe_values.reinit(begin_tria);
1929  const std::vector<Point<DoFHandlerType::space_dimension>> &points =
1930  fe_values.get_quadrature_points();
1931  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1932  if (!already_touched[local_dof_indices[i]])
1933  {
1934  support_point_list[local_dof_indices[i]].first = points[i];
1935  support_point_list[local_dof_indices[i]].second =
1936  local_dof_indices[i];
1937  already_touched[local_dof_indices[i]] = true;
1938  }
1939  }
1940 
1942  direction);
1943  std::sort(support_point_list.begin(),
1944  support_point_list.end(),
1945  comparator);
1946  for (types::global_dof_index i = 0; i < n_dofs; ++i)
1947  new_indices[support_point_list[i].second] = i;
1948  }
1949  }
1950 
1951 
1952 
1956  namespace internal
1957  {
1958  template <int dim>
1959  struct ClockCells
1960  {
1964  const Point<dim> &center;
1968  bool counter;
1969 
1973  ClockCells(const Point<dim> &center, bool counter)
1974  : center(center)
1975  , counter(counter)
1976  {}
1977 
1981  template <class DHCellIterator>
1982  bool
1983  operator()(const DHCellIterator &c1, const DHCellIterator &c2) const
1984  {
1985  // dispatch to
1986  // dimension-dependent functions
1987  return compare(c1, c2, std::integral_constant<int, dim>());
1988  }
1989 
1990  private:
1994  template <class DHCellIterator, int xdim>
1995  bool
1996  compare(const DHCellIterator &c1,
1997  const DHCellIterator &c2,
1998  std::integral_constant<int, xdim>) const
1999  {
2000  const Tensor<1, dim> v1 = c1->center() - center;
2001  const Tensor<1, dim> v2 = c2->center() - center;
2002  const double s1 = std::atan2(v1[0], v1[1]);
2003  const double s2 = std::atan2(v2[0], v2[1]);
2004  return (counter ? (s1 > s2) : (s2 > s1));
2005  }
2006 
2007 
2012  template <class DHCellIterator>
2013  bool
2014  compare(const DHCellIterator &,
2015  const DHCellIterator &,
2016  std::integral_constant<int, 1>) const
2017  {
2018  Assert(dim >= 2,
2019  ExcMessage("This operation only makes sense for dim>=2."));
2020  return false;
2021  }
2022  };
2023  } // namespace internal
2024 
2025 
2026 
2027  template <typename DoFHandlerType>
2028  void
2029  clockwise_dg(DoFHandlerType & dof,
2031  const bool counter)
2032  {
2033  std::vector<types::global_dof_index> renumbering(dof.n_dofs());
2034  compute_clockwise_dg(renumbering, dof, center, counter);
2035 
2036  dof.renumber_dofs(renumbering);
2037  }
2038 
2039 
2040 
2041  template <typename DoFHandlerType>
2042  void
2043  compute_clockwise_dg(std::vector<types::global_dof_index> &new_indices,
2044  const DoFHandlerType & dof,
2046  const bool counter)
2047  {
2048  std::vector<typename DoFHandlerType::active_cell_iterator> ordered_cells;
2049  ordered_cells.reserve(dof.get_triangulation().n_active_cells());
2050  internal::ClockCells<DoFHandlerType::space_dimension> comparator(center,
2051  counter);
2052 
2053  typename DoFHandlerType::active_cell_iterator p = dof.begin_active();
2054  typename DoFHandlerType::active_cell_iterator end = dof.end();
2055 
2056  while (p != end)
2057  {
2058  ordered_cells.push_back(p);
2059  ++p;
2060  }
2061  std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
2062 
2063  std::vector<types::global_dof_index> reverse(new_indices.size());
2064  compute_cell_wise(new_indices, reverse, dof, ordered_cells);
2065  }
2066 
2067 
2068 
2069  template <typename DoFHandlerType>
2070  void
2071  clockwise_dg(DoFHandlerType & dof,
2072  const unsigned int level,
2074  const bool counter)
2075  {
2076  std::vector<typename DoFHandlerType::level_cell_iterator> ordered_cells;
2077  ordered_cells.reserve(dof.get_triangulation().n_active_cells());
2078  internal::ClockCells<DoFHandlerType::space_dimension> comparator(center,
2079  counter);
2080 
2081  typename DoFHandlerType::level_cell_iterator p = dof.begin(level);
2082  typename DoFHandlerType::level_cell_iterator end = dof.end(level);
2083 
2084  while (p != end)
2085  {
2086  ordered_cells.push_back(p);
2087  ++p;
2088  }
2089  std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
2090 
2091  cell_wise(dof, level, ordered_cells);
2092  }
2093 
2094 
2095 
2096  template <typename DoFHandlerType>
2097  void
2098  random(DoFHandlerType &dof_handler)
2099  {
2100  std::vector<types::global_dof_index> renumbering(
2101  dof_handler.n_dofs(), numbers::invalid_dof_index);
2102  compute_random(renumbering, dof_handler);
2103 
2104  dof_handler.renumber_dofs(renumbering);
2105  }
2106 
2107 
2108 
2109  template <typename DoFHandlerType>
2110  void
2111  random(DoFHandlerType &dof_handler, const unsigned int level)
2112  {
2113  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
2115 
2116  std::vector<types::global_dof_index> renumbering(
2117  dof_handler.locally_owned_mg_dofs(level).n_elements(),
2119 
2120  compute_random(renumbering, dof_handler, level);
2121 
2122  dof_handler.renumber_dofs(level, renumbering);
2123  }
2124 
2125 
2126 
2127  template <typename DoFHandlerType>
2128  void
2129  compute_random(std::vector<types::global_dof_index> &new_indices,
2130  const DoFHandlerType & dof_handler)
2131  {
2132  const types::global_dof_index n_dofs = dof_handler.n_dofs();
2133  Assert(new_indices.size() == n_dofs,
2134  ExcDimensionMismatch(new_indices.size(), n_dofs));
2135 
2136  std::iota(new_indices.begin(),
2137  new_indices.end(),
2139 
2140  // shuffle the elements; the following is essentially std::shuffle (which
2141  // is new in C++11) but with a boost URNG
2142  // we could use std::mt19937 here but doing so results in compiler-dependent
2143  // output
2144  ::boost::mt19937 random_number_generator;
2145  for (unsigned int i = 1; i < n_dofs; ++i)
2146  {
2147  // get a random number between 0 and i (inclusive)
2148  const unsigned int j =
2149  ::boost::random::uniform_int_distribution<>(0, i)(
2150  random_number_generator);
2151 
2152  // if possible, swap the elements
2153  if (i != j)
2154  std::swap(new_indices[i], new_indices[j]);
2155  }
2156  }
2157 
2158 
2159 
2160  template <typename DoFHandlerType>
2161  void
2162  compute_random(std::vector<types::global_dof_index> &new_indices,
2163  const DoFHandlerType & dof_handler,
2164  const unsigned int level)
2165  {
2166  const types::global_dof_index n_dofs = dof_handler.n_dofs(level);
2167  Assert(new_indices.size() == n_dofs,
2168  ExcDimensionMismatch(new_indices.size(), n_dofs));
2169 
2170  std::iota(new_indices.begin(),
2171  new_indices.end(),
2173 
2174  // shuffle the elements; the following is essentially std::shuffle (which
2175  // is new in C++11) but with a boost URNG
2176  // we could use std::mt19937 here but doing so results in
2177  // compiler-dependent output
2178  ::boost::mt19937 random_number_generator;
2179  for (unsigned int i = 1; i < n_dofs; ++i)
2180  {
2181  // get a random number between 0 and i (inclusive)
2182  const unsigned int j =
2183  ::boost::random::uniform_int_distribution<>(0, i)(
2184  random_number_generator);
2185 
2186  // if possible, swap the elements
2187  if (i != j)
2188  std::swap(new_indices[i], new_indices[j]);
2189  }
2190  }
2191 
2192 
2193 
2194  template <typename DoFHandlerType>
2195  void
2196  subdomain_wise(DoFHandlerType &dof_handler)
2197  {
2198  Assert(
2199  (!dynamic_cast<
2200  const parallel::TriangulationBase<DoFHandlerType::dimension,
2201  DoFHandlerType::space_dimension> *>(
2202  &dof_handler.get_triangulation())),
2203  ExcMessage(
2204  "Parallel triangulations are already enumerated according to their MPI process id."));
2205 
2206  std::vector<types::global_dof_index> renumbering(
2207  dof_handler.n_dofs(), numbers::invalid_dof_index);
2208  compute_subdomain_wise(renumbering, dof_handler);
2209 
2210  dof_handler.renumber_dofs(renumbering);
2211  }
2212 
2213 
2214 
2215  template <typename DoFHandlerType>
2216  void
2217  compute_subdomain_wise(std::vector<types::global_dof_index> &new_dof_indices,
2218  const DoFHandlerType & dof_handler)
2219  {
2220  const types::global_dof_index n_dofs = dof_handler.n_dofs();
2221  Assert(new_dof_indices.size() == n_dofs,
2222  ExcDimensionMismatch(new_dof_indices.size(), n_dofs));
2223 
2224  // first get the association of each dof
2225  // with a subdomain and determine the total
2226  // number of subdomain ids used
2227  std::vector<types::subdomain_id> subdomain_association(n_dofs);
2228  DoFTools::get_subdomain_association(dof_handler, subdomain_association);
2229  const unsigned int n_subdomains =
2230  *std::max_element(subdomain_association.begin(),
2231  subdomain_association.end()) +
2232  1;
2233 
2234  // then renumber the subdomains by first
2235  // looking at those belonging to subdomain
2236  // 0, then those of subdomain 1, etc. note
2237  // that the algorithm is stable, i.e. if
2238  // two dofs i,j have i<j and belong to the
2239  // same subdomain, then they will be in
2240  // this order also after reordering
2241  std::fill(new_dof_indices.begin(),
2242  new_dof_indices.end(),
2244  types::global_dof_index next_free_index = 0;
2245  for (types::subdomain_id subdomain = 0; subdomain < n_subdomains;
2246  ++subdomain)
2247  for (types::global_dof_index i = 0; i < n_dofs; ++i)
2248  if (subdomain_association[i] == subdomain)
2249  {
2250  Assert(new_dof_indices[i] == numbers::invalid_dof_index,
2251  ExcInternalError());
2252  new_dof_indices[i] = next_free_index;
2253  ++next_free_index;
2254  }
2255 
2256  // we should have numbered all dofs
2257  Assert(next_free_index == n_dofs, ExcInternalError());
2258  Assert(std::find(new_dof_indices.begin(),
2259  new_dof_indices.end(),
2260  numbers::invalid_dof_index) == new_dof_indices.end(),
2261  ExcInternalError());
2262  }
2263 
2264 } // namespace DoFRenumbering
2265 
2266 
2267 
2268 /*-------------- Explicit Instantiations -------------------------------*/
2269 #include "dof_renumbering.inst"
2270 
2271 
2272 DEAL_II_NAMESPACE_CLOSE
void compute_downstream(std::vector< types::global_dof_index > &new_dof_indices, std::vector< types::global_dof_index > &reverse, const DoFHandlerType &dof_handler, const Tensor< 1, DoFHandlerType::space_dimension > &direction, const bool dof_wise_renumbering)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1571
void downstream(DoFHandlerType &dof_handler, const Tensor< 1, DoFHandlerType::space_dimension > &direction, const bool dof_wise_renumbering=false)
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:931
void minimum_degree(DoFHandlerType &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
void compute_cell_wise(std::vector< types::global_dof_index > &renumbering, std::vector< types::global_dof_index > &inverse_renumbering, const DoFHandlerType &dof_handler, const std::vector< typename DoFHandlerType::active_cell_iterator > &cell_order)
void compute_king_ordering(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &, const bool reversed_numbering=false, const bool use_constraints=false)
void sort_selected_dofs_back(DoFHandlerType &dof_handler, const std::vector< bool > &selected_dofs)
void get_subdomain_association(const DoFHandlerType &dof_handler, std::vector< types::subdomain_id > &subdomain)
Definition: dof_tools.cc:1552
cell_iterator end() const
Definition: dof_handler.cc:960
#define AssertIndexRange(index, range)
Definition: exceptions.h:1641
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:567
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 >())
active_cell_iterator begin_active(const unsigned int level=0) const
const FiniteElement< dim, spacedim > & get_fe() const
void component_wise(DoFHandlerType &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
STL namespace.
unsigned int size() const
Transformed quadrature points.
types::global_dof_index compute_component_wise(std::vector< types::global_dof_index > &new_dof_indices, const CellIterator &start, const typename identity< CellIterator >::type &end, const std::vector< unsigned int > &target_component, const bool is_level_operation)
bool is_primitive() const
Definition: fe.h:3301
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
LinearAlgebra::distributed::Vector< Number > Vector
static ::ExceptionBase & ExcDoFHandlerNotInitialized()
size_type size() const
Definition: index_set.h:1625
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:312
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:944
const std::vector< Point< spacedim > > & get_quadrature_points() const
unsigned int n_locally_owned_dofs() const
void random(DoFHandlerType &dof_handler)
static ::ExceptionBase & ExcMessage(std::string arg1)
void compute_minimum_degree(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &, const bool reversed_numbering=false, const bool use_constraints=false)
unsigned int subdomain_id
Definition: types.h:43
void Cuthill_McKee(DoFHandlerType &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
T sum(const T &t, const MPI_Comm &mpi_communicator)
void subtract_set(const IndexSet &other)
Definition: index_set.cc:238
#define Assert(cond, exc)
Definition: exceptions.h:1411
void compute_subdomain_wise(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1912
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:3281
types::global_dof_index n_dofs() const
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1183
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType, 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:132
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
void compute_random(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler)
void compute_sort_selected_dofs_back(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler, const std::vector< bool > &selected_dofs)
const unsigned int dofs_per_cell
Definition: fe_base.h:282
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
void swap(Vector< Number > &u, Vector< Number > &v)
Definition: vector.h:1376
void reinit(const IndexSet &local_constraints=IndexSet())
unsigned int global_dof_index
Definition: types.h:89
void extract_locally_active_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1148
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:3110
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1699
void compute_clockwise_dg(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler, const Point< DoFHandlerType::space_dimension > &center, const bool counter)
Definition: tensor.h:422
void king_ordering(DoFHandlerType &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
Definition: fe.h:3253
void hierarchical(DoFHandlerType &dof_handler)
void block_wise(DoFHandler< dim, spacedim > &dof_handler)
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
Definition: fe.h:38
void compute_Cuthill_McKee(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &, const bool reversed_numbering=false, const bool use_constraints=false)
types::global_dof_index n_dofs() const
std::vector< unsigned int > reverse_permutation(const std::vector< unsigned int > &permutation)
Definition: utilities.cc:861
static ::ExceptionBase & ExcNotImplemented()
bool is_element(const size_type index) const
Definition: index_set.h:1756
void subdomain_wise(DoFHandlerType &dof_handler)
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:329
void clockwise_dg(DoFHandlerType &dof_handler, const Point< DoFHandlerType::space_dimension > &center, const bool counter=false)
const types::global_dof_index invalid_dof_index
Definition: types.h:202
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
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)
size_type n_elements() const
Definition: index_set.h:1823
void push_back(const Quadrature< dim > &new_quadrature)
Definition: q_collection.h:221
cell_iterator end() const
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1871
static ::ExceptionBase & ExcInternalError()
void cell_wise(DoFHandlerType &dof_handler, const std::vector< typename DoFHandlerType::active_cell_iterator > &cell_order)
const ::FEValues< dim, spacedim > & get_present_fe_values() const
Definition: fe_values.h:569