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