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