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