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