Reference documentation for deal.II version 8.4.1
dof_renumbering.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2015 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/thread_management.h>
17 #include <deal.II/base/utilities.h>
18 #include <deal.II/base/quadrature_lib.h>
19 #include <deal.II/base/types.h>
20 #include <deal.II/base/template_constraints.h>
21 
22 #include <deal.II/lac/sparsity_pattern.h>
23 #include <deal.II/lac/sparsity_tools.h>
24 #include <deal.II/lac/dynamic_sparsity_pattern.h>
25 #include <deal.II/lac/constraint_matrix.h>
26 
27 #include <deal.II/dofs/dof_accessor.h>
28 #include <deal.II/dofs/dof_handler.h>
29 #include <deal.II/dofs/dof_tools.h>
30 #include <deal.II/dofs/dof_renumbering.h>
31 
32 #include <deal.II/grid/tria_iterator.h>
33 #include <deal.II/grid/tria.h>
34 
35 #include <deal.II/fe/fe.h>
36 #include <deal.II/hp/dof_handler.h>
37 #include <deal.II/hp/fe_collection.h>
38 #include <deal.II/hp/fe_values.h>
39 
40 #include <deal.II/multigrid/mg_tools.h>
41 
42 #include <deal.II/distributed/tria.h>
43 
44 #include <boost/config.hpp>
45 #include <boost/graph/adjacency_list.hpp>
46 #include <boost/graph/cuthill_mckee_ordering.hpp>
47 #include <boost/graph/king_ordering.hpp>
48 #include <boost/graph/minimum_degree_ordering.hpp>
49 #include <boost/graph/properties.hpp>
50 #include <boost/graph/bandwidth.hpp>
51 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
52 #include <boost/random.hpp>
53 #include <boost/random/uniform_int_distribution.hpp>
55 
56 #include <vector>
57 #include <map>
58 #include <algorithm>
59 #include <cmath>
60 #include <functional>
61 
62 
63 DEAL_II_NAMESPACE_OPEN
64 
65 namespace DoFRenumbering
66 {
67  namespace boost
68  {
69  namespace boosttypes
70  {
71  using namespace ::boost;
72  using namespace std;
73 
74  typedef adjacency_list<vecS, vecS, undirectedS,
75  property<vertex_color_t, default_color_type,
76  property<vertex_degree_t,int> > > Graph;
77  typedef graph_traits<Graph>::vertex_descriptor Vertex;
78  typedef graph_traits<Graph>::vertices_size_type size_type;
79 
80  typedef std::pair<size_type, size_type> Pair;
81  }
82 
83 
84  namespace internal
85  {
86  template <typename DoFHandlerType>
87  void create_graph
88  (const DoFHandlerType &dof_handler,
89  const bool use_constraints,
90  boosttypes::Graph &graph,
91  boosttypes::property_map<boosttypes::Graph,boosttypes::vertex_degree_t>::type &graph_degree)
92  {
93  {
94  // create intermediate sparsity pattern
95  // (faster than directly submitting
96  // indices)
97  ConstraintMatrix constraints;
98  if (use_constraints)
99  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
100  constraints.close ();
101  DynamicSparsityPattern dsp (dof_handler.n_dofs(),
102  dof_handler.n_dofs());
103  DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints);
104 
105  // submit the entries to the boost graph
106  for (unsigned int row=0; row<dsp.n_rows(); ++row)
107  for (unsigned int col=0; col < dsp.row_length(row); ++col)
108  add_edge (row, dsp.column_number (row, col), graph);
109  }
110 
111  boosttypes::graph_traits<boosttypes::Graph>::vertex_iterator ui, ui_end;
112 
113  graph_degree = get(::boost::vertex_degree, graph);
114  for (::boost::tie(ui, ui_end) = vertices(graph); ui != ui_end; ++ui)
115  graph_degree[*ui] = degree(*ui, graph);
116  }
117  }
118 
119 
120  template <typename DoFHandlerType>
121  void
122  Cuthill_McKee (DoFHandlerType &dof_handler,
123  const bool reversed_numbering,
124  const bool use_constraints)
125  {
126  std::vector<types::global_dof_index> renumbering(dof_handler.n_dofs(),
127  DoFHandlerType::invalid_dof_index);
128  compute_Cuthill_McKee(renumbering, dof_handler, reversed_numbering,
129  use_constraints);
130 
131  // actually perform renumbering;
132  // this is dimension specific and
133  // thus needs an own function
134  dof_handler.renumber_dofs (renumbering);
135  }
136 
137 
138  template <typename DoFHandlerType>
139  void
140  compute_Cuthill_McKee (std::vector<types::global_dof_index> &new_dof_indices,
141  const DoFHandlerType &dof_handler,
142  const bool reversed_numbering,
143  const bool use_constraints)
144  {
145  boosttypes::Graph
146  graph(dof_handler.n_dofs());
147  boosttypes::property_map<boosttypes::Graph,boosttypes::vertex_degree_t>::type
148  graph_degree;
149 
150  internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
151 
152  boosttypes::property_map<boosttypes::Graph, boosttypes::vertex_index_t>::type
153  index_map = get(::boost::vertex_index, graph);
154 
155 
156  std::vector<boosttypes::Vertex> inv_perm(num_vertices(graph));
157 
158  if (reversed_numbering == false)
159  ::boost::cuthill_mckee_ordering(graph, inv_perm.rbegin(),
160  get(::boost::vertex_color, graph),
161  make_degree_map(graph));
162  else
163  ::boost::cuthill_mckee_ordering(graph, inv_perm.begin(),
164  get(::boost::vertex_color, graph),
165  make_degree_map(graph));
166 
167  for (boosttypes::size_type c = 0; c != inv_perm.size(); ++c)
168  new_dof_indices[index_map[inv_perm[c]]] = c;
169 
170  Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
171  DoFHandlerType::invalid_dof_index) == new_dof_indices.end(),
172  ExcInternalError());
173  }
174 
175 
176 
177  template <typename DoFHandlerType>
178  void
179  king_ordering (DoFHandlerType &dof_handler,
180  const bool reversed_numbering,
181  const bool use_constraints)
182  {
183  std::vector<types::global_dof_index> renumbering(dof_handler.n_dofs(),
184  DoFHandlerType::invalid_dof_index);
185  compute_king_ordering(renumbering, dof_handler, reversed_numbering,
186  use_constraints);
187 
188  // actually perform renumbering;
189  // this is dimension specific and
190  // thus needs an own function
191  dof_handler.renumber_dofs (renumbering);
192  }
193 
194 
195  template <typename DoFHandlerType>
196  void
197  compute_king_ordering (std::vector<types::global_dof_index> &new_dof_indices,
198  const DoFHandlerType &dof_handler,
199  const bool reversed_numbering,
200  const bool use_constraints)
201  {
202  boosttypes::Graph
203  graph(dof_handler.n_dofs());
204  boosttypes::property_map<boosttypes::Graph,boosttypes::vertex_degree_t>::type
205  graph_degree;
206 
207  internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
208 
209  boosttypes::property_map<boosttypes::Graph, boosttypes::vertex_index_t>::type
210  index_map = get(::boost::vertex_index, graph);
211 
212 
213  std::vector<boosttypes::Vertex> inv_perm(num_vertices(graph));
214 
215  if (reversed_numbering == false)
216  ::boost::king_ordering(graph, inv_perm.rbegin());
217  else
218  ::boost::king_ordering(graph, inv_perm.begin());
219 
220  for (boosttypes::size_type c = 0; c != inv_perm.size(); ++c)
221  new_dof_indices[index_map[inv_perm[c]]] = c;
222 
223  Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
224  DoFHandlerType::invalid_dof_index) == new_dof_indices.end(),
225  ExcInternalError());
226  }
227 
228 
229 
230  template <typename DoFHandlerType>
231  void
232  minimum_degree (DoFHandlerType &dof_handler,
233  const bool reversed_numbering,
234  const bool use_constraints)
235  {
236  std::vector<types::global_dof_index> renumbering(dof_handler.n_dofs(),
237  DoFHandlerType::invalid_dof_index);
238  compute_minimum_degree(renumbering, dof_handler, reversed_numbering,
239  use_constraints);
240 
241  // actually perform renumbering;
242  // this is dimension specific and
243  // thus needs an own function
244  dof_handler.renumber_dofs (renumbering);
245  }
246 
247 
248  template <typename DoFHandlerType>
249  void
250  compute_minimum_degree (std::vector<types::global_dof_index> &new_dof_indices,
251  const DoFHandlerType &dof_handler,
252  const bool reversed_numbering,
253  const bool use_constraints)
254  {
255  (void)use_constraints;
256  Assert (use_constraints == false, ExcNotImplemented());
257 
258  // the following code is pretty
259  // much a verbatim copy of the
260  // sample code for the
261  // minimum_degree_ordering manual
262  // page from the BOOST Graph
263  // Library
264  using namespace ::boost;
265 
266  int delta = 0;
267 
268  typedef double Type;
269 
270  // must be BGL directed graph now
271  typedef adjacency_list<vecS, vecS, directedS> Graph;
272  typedef graph_traits<Graph>::vertex_descriptor Vertex;
273 
274  int n = dof_handler.n_dofs();
275 
276  Graph G(n);
277 
278  std::vector<::types::global_dof_index> dofs_on_this_cell;
279 
280  typename DoFHandlerType::active_cell_iterator cell = dof_handler.begin_active(),
281  endc = dof_handler.end();
282 
283  for (; cell!=endc; ++cell)
284  {
285 
286  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
287 
288  dofs_on_this_cell.resize (dofs_per_cell);
289 
290  cell->get_active_or_mg_dof_indices (dofs_on_this_cell);
291  for (unsigned int i=0; i<dofs_per_cell; ++i)
292  for (unsigned int j=0; j<dofs_per_cell; ++j)
293  if (dofs_on_this_cell[i] > dofs_on_this_cell[j])
294  {
295  add_edge (dofs_on_this_cell[i], dofs_on_this_cell[j], G);
296  add_edge (dofs_on_this_cell[j], dofs_on_this_cell[i], G);
297  }
298  }
299 
300 
301  typedef std::vector<int> Vector;
302 
303 
304  Vector inverse_perm(n, 0);
305 
306  Vector perm(n, 0);
307 
308 
309  Vector supernode_sizes(n, 1);
310  // init has to be 1
311 
312  ::boost::property_map<Graph, vertex_index_t>::type
313  id = get(vertex_index, G);
314 
315 
316  Vector degree(n, 0);
317 
318 
319  minimum_degree_ordering
320  (G,
321  make_iterator_property_map(&degree[0], id, degree[0]),
322  &inverse_perm[0],
323  &perm[0],
324  make_iterator_property_map(&supernode_sizes[0], id, supernode_sizes[0]),
325  delta, id);
326 
327 
328  for (int i=0; i<n; ++i)
329  {
330  Assert (std::find (perm.begin(), perm.end(), i)
331  != perm.end(),
332  ExcInternalError());
333  Assert (std::find (inverse_perm.begin(), inverse_perm.end(), i)
334  != inverse_perm.end(),
335  ExcInternalError());
336  Assert (inverse_perm[perm[i]] == i, ExcInternalError());
337  }
338 
339  if (reversed_numbering == true)
340  std::copy (perm.begin(), perm.end(),
341  new_dof_indices.begin());
342  else
343  std::copy (inverse_perm.begin(), inverse_perm.end(),
344  new_dof_indices.begin());
345  }
346 
347  } // namespace boost
348 
349 
350 
351  template <typename DoFHandlerType>
352  void
353  Cuthill_McKee (DoFHandlerType &dof_handler,
354  const bool reversed_numbering,
355  const bool use_constraints,
356  const std::vector<types::global_dof_index> &starting_indices)
357  {
358  std::vector<types::global_dof_index> renumbering(dof_handler.locally_owned_dofs().n_elements(),
359  DoFHandlerType::invalid_dof_index);
360  compute_Cuthill_McKee(renumbering, dof_handler, reversed_numbering,
361  use_constraints, starting_indices);
362 
363  // actually perform renumbering;
364  // this is dimension specific and
365  // thus needs an own function
366  dof_handler.renumber_dofs (renumbering);
367  }
368 
369 
370 
371  template <typename DoFHandlerType>
372  void
373  compute_Cuthill_McKee (std::vector<types::global_dof_index> &new_indices,
374  const DoFHandlerType &dof_handler,
375  const bool reversed_numbering,
376  const bool use_constraints,
377  const std::vector<types::global_dof_index> &starting_indices)
378  {
379  // see if there is anything to do at all or whether we can skip the work on this processor
380  if (dof_handler.locally_owned_dofs().n_elements() == 0)
381  {
382  Assert (new_indices.size() == 0, ExcInternalError());
383  return;
384  }
385 
386  // make the connection graph. in 2d/3d use an intermediate compressed
387  // sparsity pattern since the we don't have very good estimates for
388  // max_couplings_between_dofs() in 3d and this then leads to excessive
389  // memory consumption
390  //
391  // note that if constraints are not requested, then the 'constraints'
392  // object will be empty and nothing happens
393  ConstraintMatrix constraints;
394  if (use_constraints)
395  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
396  constraints.close ();
397 
398  const IndexSet locally_owned = dof_handler.locally_owned_dofs();
399 
400  // otherwise compute the Cuthill-McKee permutation
401  DynamicSparsityPattern dsp (dof_handler.n_dofs(),
402  dof_handler.n_dofs(),
403  locally_owned);
404  DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints);
405 
406  // constraints are not needed anymore
407  constraints.clear ();
408 
409  // If the index set is not complete, need to get indices in local index
410  // space.
411  if (locally_owned.n_elements() !=
412  locally_owned.size())
413  {
414  // Create sparsity pattern from dsp by transferring its indices to
415  // processor-local index space and doing Cuthill-McKee there
416  DynamicSparsityPattern sparsity(locally_owned.n_elements(),
417  locally_owned.n_elements());
418  std::vector<types::global_dof_index> row_entries;
419  for (unsigned int i=0; i<locally_owned.n_elements(); ++i)
420  {
421  const types::global_dof_index row = locally_owned.nth_index_in_set(i);
422  row_entries.clear();
424  dsp.begin(row); it != dsp.end(row); ++it)
425  if (it->column() != row && locally_owned.is_element(it->column()))
426  row_entries.push_back(locally_owned.index_within_set(it->column()));
427  sparsity.add_entries(i, row_entries.begin(), row_entries.end(),
428  true);
429  }
430 
431  // translate starting indices from global to local indices
432  std::vector<types::global_dof_index> local_starting_indices (starting_indices.size());
433  for (unsigned int i=0; i<starting_indices.size(); ++i)
434  {
435  Assert (locally_owned.is_element (starting_indices[i]),
436  ExcMessage ("You specified global degree of freedom "
437  + Utilities::to_string(starting_indices[i]) +
438  " as a starting index, but this index is not among the "
439  "locally owned ones on this processor."));
440  local_starting_indices[i] = locally_owned.index_within_set(starting_indices[i]);
441  }
442 
443  // then do the renumbering on the locally owned portion
444  AssertDimension(new_indices.size(), locally_owned.n_elements());
445  SparsityTools::reorder_Cuthill_McKee (sparsity, new_indices,
446  local_starting_indices);
447  if (reversed_numbering)
448  new_indices = Utilities::reverse_permutation (new_indices);
449 
450  // convert indices back to global index space
451  for (std::size_t i=0; i<new_indices.size(); ++i)
452  new_indices[i] = locally_owned.nth_index_in_set(new_indices[i]);
453  }
454  else
455  {
456  AssertDimension(new_indices.size(), dsp.n_rows());
457  SparsityTools::reorder_Cuthill_McKee (dsp, new_indices,
458  starting_indices);
459  if (reversed_numbering)
460  new_indices = Utilities::reverse_permutation (new_indices);
461  }
462  }
463 
464 
465 
466  template <typename DoFHandlerType>
467  void Cuthill_McKee (DoFHandlerType &dof_handler,
468  const unsigned int level,
469  const bool reversed_numbering,
470  const std::vector<types::global_dof_index> &starting_indices)
471  {
472  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
473  ExcNotInitialized());
474 
475  // make the connection graph
476  DynamicSparsityPattern dsp (dof_handler.n_dofs(level),
477  dof_handler.n_dofs(level));
478  MGTools::make_sparsity_pattern (dof_handler, dsp, level);
479 
480  std::vector<types::global_dof_index> new_indices(dsp.n_rows());
481  SparsityTools::reorder_Cuthill_McKee (dsp, new_indices,
482  starting_indices);
483 
484  if (reversed_numbering)
485  new_indices = Utilities::reverse_permutation (new_indices);
486 
487  // actually perform renumbering;
488  // this is dimension specific and
489  // thus needs an own function
490  dof_handler.renumber_dofs (level, new_indices);
491  }
492 
493 
494 
495  template <int dim, int spacedim>
496  void
498  const std::vector<unsigned int> &component_order_arg)
499  {
500  std::vector<types::global_dof_index> renumbering (dof_handler.n_locally_owned_dofs(),
502 
504  start = dof_handler.begin_active();
505  const typename DoFHandler<dim,spacedim>::level_cell_iterator
506  end = dof_handler.end();
507 
508  const types::global_dof_index result =
509  compute_component_wise<dim,spacedim,
511  typename DoFHandler<dim,spacedim>::level_cell_iterator>
512  (renumbering, start, end, component_order_arg, false);
513  if (result == 0)
514  return;
515 
516  // verify that the last numbered
517  // degree of freedom is either
518  // equal to the number of degrees
519  // of freedom in total (the
520  // sequential case) or in the
521  // distributed case at least
522  // makes sense
523  Assert ((result == dof_handler.n_locally_owned_dofs())
524  ||
525  ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs())
526  &&
527  (result <= dof_handler.n_dofs())),
528  ExcRenumberingIncomplete());
529 
530  dof_handler.renumber_dofs (renumbering);
531 
532  // for (unsigned int level=0;level<dof_handler.get_triangulation().n_levels();++level)
533  // if (dof_handler.n_dofs(level) != numbers::invalid_dof_index)
534  // component_wise(dof_handler, level, component_order_arg);
535  }
536 
537 
538 
539  template <int dim>
540  void
542  const std::vector<unsigned int> &component_order_arg)
543  {
544 //TODO: Merge with previous function
545  std::vector<types::global_dof_index> renumbering (dof_handler.n_dofs(),
547 
549  start = dof_handler.begin_active();
550  const typename hp::DoFHandler<dim>::level_cell_iterator
551  end = dof_handler.end();
552 
553  const types::global_dof_index result =
554  compute_component_wise<hp::DoFHandler<dim>::dimension,hp::DoFHandler<dim>::space_dimension,
556  typename hp::DoFHandler<dim>::level_cell_iterator>
557  (renumbering, start, end, component_order_arg, false);
558 
559  if (result == 0) return;
560 
561  Assert (result == dof_handler.n_dofs(),
562  ExcRenumberingIncomplete());
563 
564  dof_handler.renumber_dofs (renumbering);
565  }
566 
567 
568 
569  template <typename DoFHandlerType>
570  void
571  component_wise (DoFHandlerType &dof_handler,
572  const unsigned int level,
573  const std::vector<unsigned int> &component_order_arg)
574  {
575  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
576  ExcNotInitialized());
577 
578  std::vector<types::global_dof_index> renumbering (dof_handler.n_dofs(level),
579  DoFHandlerType::invalid_dof_index);
580 
581  typename DoFHandlerType::level_cell_iterator start =dof_handler.begin(level);
582  typename DoFHandlerType::level_cell_iterator end = dof_handler.end(level);
583 
584  const types::global_dof_index result =
585  compute_component_wise<DoFHandlerType::dimension, DoFHandlerType::space_dimension,
586  typename DoFHandlerType::level_cell_iterator, typename DoFHandlerType::level_cell_iterator>
587  (renumbering, start, end, component_order_arg, true);
588 
589  if (result == 0) return;
590 
591  Assert (result == dof_handler.n_dofs(level),
592  ExcRenumberingIncomplete());
593 
594  if (renumbering.size()!=0)
595  dof_handler.renumber_dofs (level, renumbering);
596  }
597 
598 
599 
600  template <int dim, int spacedim, class ITERATOR, class ENDITERATOR>
602  compute_component_wise (std::vector<types::global_dof_index> &new_indices,
603  const ITERATOR &start,
604  const ENDITERATOR &end,
605  const std::vector<unsigned int> &component_order_arg,
606  bool is_level_operation)
607  {
609  fe_collection (start->get_dof_handler().get_fe ());
610 
611  // do nothing if the FE has only
612  // one component
613  if (fe_collection.n_components() == 1)
614  {
615  new_indices.resize(0);
616  return 0;
617  }
618 
619  // Copy last argument into a
620  // writable vector.
621  std::vector<unsigned int> component_order (component_order_arg);
622  // If the last argument was an
623  // empty vector, set up things to
624  // store components in the order
625  // found in the system.
626  if (component_order.size() == 0)
627  for (unsigned int i=0; i<fe_collection.n_components(); ++i)
628  component_order.push_back (i);
629 
630  Assert (component_order.size() == fe_collection.n_components(),
631  ExcDimensionMismatch(component_order.size(), fe_collection.n_components()));
632 
633  for (unsigned int i=0; i<component_order.size(); ++i)
634  Assert(component_order[i] < fe_collection.n_components(),
635  ExcIndexRange(component_order[i], 0, fe_collection.n_components()));
636 
637  // vector to hold the dof indices on
638  // the cell we visit at a time
639  std::vector<types::global_dof_index> local_dof_indices;
640 
641  // prebuilt list to which component
642  // a given dof on a cell
643  // should go. note that we get into
644  // trouble here if the shape
645  // function is not primitive, since
646  // then there is no single vector
647  // component to which it
648  // belongs. in this case, assign it
649  // to the first vector component to
650  // which it belongs
651  std::vector<std::vector<unsigned int> > component_list (fe_collection.size());
652  for (unsigned int f=0; f<fe_collection.size(); ++f)
653  {
654  const FiniteElement<dim,spacedim> &fe = fe_collection[f];
655  const unsigned int dofs_per_cell = fe.dofs_per_cell;
656  component_list[f].resize(dofs_per_cell);
657  for (unsigned int i=0; i<dofs_per_cell; ++i)
658  if (fe.is_primitive(i))
659  component_list[f][i]
660  = component_order[fe.system_to_component_index(i).first];
661  else
662  {
663  const unsigned int comp
665 
666  // then associate this degree
667  // of freedom with this
668  // component
669  component_list[f][i] = component_order[comp];
670  }
671  }
672 
673  // set up a map where for each
674  // component the respective degrees
675  // of freedom are collected.
676  //
677  // note that this map is sorted by
678  // component but that within each
679  // component it is NOT sorted by
680  // dof index. note also that some
681  // dof indices are entered
682  // multiply, so we will have to
683  // take care of that
684  std::vector<std::vector<types::global_dof_index> >
685  component_to_dof_map (fe_collection.n_components());
686  for (ITERATOR cell=start; cell!=end; ++cell)
687  {
688  if (is_level_operation)
689  {
690  //we are dealing with mg dofs, skip foreign level cells:
691  if (!cell->is_locally_owned_on_level())
692  continue;
693  }
694  else
695  {
696  //we are dealing with active dofs, skip the loop if not locally
697  // owned:
698  if (!cell->is_locally_owned())
699  continue;
700  }
701  // on each cell: get dof indices
702  // and insert them into the global
703  // list using their component
704  const unsigned int fe_index = cell->active_fe_index();
705  const unsigned int dofs_per_cell = fe_collection[fe_index].dofs_per_cell;
706  local_dof_indices.resize (dofs_per_cell);
707  cell->get_active_or_mg_dof_indices (local_dof_indices);
708  for (unsigned int i=0; i<dofs_per_cell; ++i)
709  if (start->get_dof_handler().locally_owned_dofs().is_element(local_dof_indices[i]))
710  component_to_dof_map[component_list[fe_index][i]].
711  push_back (local_dof_indices[i]);
712  }
713 
714  // now we've got all indices sorted
715  // into buckets labeled by their
716  // target component number. we've
717  // only got to traverse this list
718  // and assign the new indices
719  //
720  // however, we first want to sort
721  // the indices entered into the
722  // buckets to preserve the order
723  // within each component and during
724  // this also remove duplicate
725  // entries
726  //
727  // note that we no longer have to
728  // care about non-primitive shape
729  // functions since the buckets
730  // corresponding to the second and
731  // following vector components of a
732  // non-primitive FE will simply be
733  // empty, everything being shoved
734  // into the first one. The same
735  // holds if several components were
736  // joined into a single target.
737  for (unsigned int component=0; component<fe_collection.n_components();
738  ++component)
739  {
740  std::sort (component_to_dof_map[component].begin(),
741  component_to_dof_map[component].end());
742  component_to_dof_map[component]
743  .erase (std::unique (component_to_dof_map[component].begin(),
744  component_to_dof_map[component].end()),
745  component_to_dof_map[component].end());
746  }
747 
748  // calculate the number of locally owned
749  // DoFs per bucket
750  const unsigned int n_buckets = fe_collection.n_components();
751  std::vector<types::global_dof_index> shifts(n_buckets);
752 
754  = (dynamic_cast<const parallel::Triangulation<dim,spacedim>*>
755  (&start->get_dof_handler().get_triangulation())))
756  {
757 #ifdef DEAL_II_WITH_MPI
758  std::vector<types::global_dof_index> local_dof_count(n_buckets);
759 
760  for (unsigned int c=0; c<n_buckets; ++c)
761  local_dof_count[c] = component_to_dof_map[c].size();
762 
763 
764  // gather information from all CPUs
765  std::vector<types::global_dof_index>
766  all_dof_counts(fe_collection.n_components() *
767  Utilities::MPI::n_mpi_processes (tria->get_communicator()));
768 
769  MPI_Allgather ( &local_dof_count[0],
770  n_buckets, DEAL_II_DOF_INDEX_MPI_TYPE,
771  &all_dof_counts[0],
772  n_buckets, DEAL_II_DOF_INDEX_MPI_TYPE,
773  tria->get_communicator());
774 
775  for (unsigned int i=0; i<n_buckets; ++i)
776  Assert (all_dof_counts[n_buckets*tria->locally_owned_subdomain()+i]
777  ==
778  local_dof_count[i],
779  ExcInternalError());
780 
781  //calculate shifts
782  unsigned int cumulated = 0;
783  for (unsigned int c=0; c<n_buckets; ++c)
784  {
785  shifts[c]=cumulated;
786  for (types::subdomain_id i=0; i<tria->locally_owned_subdomain(); ++i)
787  shifts[c] += all_dof_counts[c+n_buckets*i];
788  for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes (tria->get_communicator()); ++i)
789  cumulated += all_dof_counts[c+n_buckets*i];
790  }
791 #else
792  (void)tria;
793  Assert (false, ExcInternalError());
794 #endif
795  }
796  else
797  {
798  shifts[0] = 0;
799  for (unsigned int c=1; c<fe_collection.n_components(); ++c)
800  shifts[c] = shifts[c-1] + component_to_dof_map[c-1].size();
801  }
802 
803 
804 
805 
806  // now concatenate all the
807  // components in the order the user
808  // desired to see
809  types::global_dof_index next_free_index = 0;
810  for (unsigned int component=0; component<fe_collection.n_components(); ++component)
811  {
812  const typename std::vector<types::global_dof_index>::const_iterator
813  begin_of_component = component_to_dof_map[component].begin(),
814  end_of_component = component_to_dof_map[component].end();
815 
816  next_free_index = shifts[component];
817 
818  for (typename std::vector<types::global_dof_index>::const_iterator
819  dof_index = begin_of_component;
820  dof_index != end_of_component; ++dof_index)
821  {
822  Assert (start->get_dof_handler().locally_owned_dofs()
823  .index_within_set(*dof_index)
824  <
825  new_indices.size(),
826  ExcInternalError());
827  new_indices[start->get_dof_handler().locally_owned_dofs()
828  .index_within_set(*dof_index)]
829  = next_free_index++;
830  }
831  }
832 
833  return next_free_index;
834  }
835 
836 
837 
838  template <int dim, int spacedim>
839  void
841  {
842  std::vector<types::global_dof_index> renumbering (dof_handler.n_locally_owned_dofs(),
844 
846  start = dof_handler.begin_active();
847  const typename DoFHandler<dim,spacedim>::level_cell_iterator
848  end = dof_handler.end();
849 
850  const types::global_dof_index result =
851  compute_block_wise<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,
852  typename DoFHandler<dim,spacedim>::level_cell_iterator>
853  (renumbering, start, end, false);
854  if (result == 0)
855  return;
856 
857  // verify that the last numbered
858  // degree of freedom is either
859  // equal to the number of degrees
860  // of freedom in total (the
861  // sequential case) or in the
862  // distributed case at least
863  // makes sense
864  Assert ((result == dof_handler.n_locally_owned_dofs())
865  ||
866  ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs())
867  &&
868  (result <= dof_handler.n_dofs())),
869  ExcRenumberingIncomplete());
870 
871  dof_handler.renumber_dofs (renumbering);
872  }
873 
874 
875 
876  template <int dim, int spacedim>
877  void
879  {
880  std::vector<types::global_dof_index> renumbering (dof_handler.n_dofs(),
882 
884  start = dof_handler.begin_active();
885  const typename hp::DoFHandler<dim,spacedim>::level_cell_iterator
886  end = dof_handler.end();
887 
888  const types::global_dof_index result =
889  compute_block_wise<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,
890  typename hp::DoFHandler<dim,spacedim>::level_cell_iterator>(renumbering,
891  start, end, false);
892 
893  if (result == 0)
894  return;
895 
896  Assert (result == dof_handler.n_dofs(),
897  ExcRenumberingIncomplete());
898 
899  dof_handler.renumber_dofs (renumbering);
900  }
901 
902 
903 
904  template <int dim, int spacedim>
905  void
906  block_wise (DoFHandler<dim,spacedim> &dof_handler, const unsigned int level)
907  {
908  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
909  ExcNotInitialized());
910 
911  std::vector<types::global_dof_index> renumbering (dof_handler.n_dofs(level),
913 
914  typename DoFHandler<dim, spacedim>::level_cell_iterator
915  start =dof_handler.begin(level);
916  typename DoFHandler<dim, spacedim>::level_cell_iterator
917  end = dof_handler.end(level);
918 
919  const types::global_dof_index result =
920  compute_block_wise<dim, spacedim, typename DoFHandler<dim, spacedim>::level_cell_iterator,
921  typename DoFHandler<dim, spacedim>::level_cell_iterator>(
922  renumbering, start, end, true);
923 
924  if (result == 0) return;
925 
926  Assert (result == dof_handler.n_dofs(level),
927  ExcRenumberingIncomplete());
928 
929  if (renumbering.size()!=0)
930  dof_handler.renumber_dofs (level, renumbering);
931  }
932 
933 
934 
935  template <int dim, int spacedim, class ITERATOR, class ENDITERATOR>
937  compute_block_wise (std::vector<types::global_dof_index> &new_indices,
938  const ITERATOR &start,
939  const ENDITERATOR &end,
940  const bool is_level_operation)
941  {
943  fe_collection (start->get_dof_handler().get_fe ());
944 
945  // do nothing if the FE has only
946  // one component
947  if (fe_collection.n_blocks() == 1)
948  {
949  new_indices.resize(0);
950  return 0;
951  }
952 
953  // vector to hold the dof indices on
954  // the cell we visit at a time
955  std::vector<types::global_dof_index> local_dof_indices;
956 
957  // prebuilt list to which block
958  // a given dof on a cell
959  // should go.
960  std::vector<std::vector<types::global_dof_index> > block_list (fe_collection.size());
961  for (unsigned int f=0; f<fe_collection.size(); ++f)
962  {
963  const FiniteElement<dim,spacedim> &fe = fe_collection[f];
964  block_list[f].resize(fe.dofs_per_cell);
965  for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
966  block_list[f][i]
967  = fe.system_to_block_index(i).first;
968  }
969 
970  // set up a map where for each
971  // block the respective degrees
972  // of freedom are collected.
973  //
974  // note that this map is sorted by
975  // block but that within each
976  // block it is NOT sorted by
977  // dof index. note also that some
978  // dof indices are entered
979  // multiply, so we will have to
980  // take care of that
981  std::vector<std::vector<types::global_dof_index> >
982  block_to_dof_map (fe_collection.n_blocks());
983  for (ITERATOR cell=start; cell!=end; ++cell)
984  {
985  if (is_level_operation)
986  {
987  //we are dealing with mg dofs, skip foreign level cells:
988  if (!cell->is_locally_owned_on_level())
989  continue;
990  }
991  else
992  {
993  //we are dealing with active dofs, skip the loop if not locally
994  // owned:
995  if (!cell->is_locally_owned())
996  continue;
997  }
998 
999  // on each cell: get dof indices
1000  // and insert them into the global
1001  // list using their component
1002  const unsigned int fe_index = cell->active_fe_index();
1003  const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell;
1004  local_dof_indices.resize (dofs_per_cell);
1005  cell->get_active_or_mg_dof_indices (local_dof_indices);
1006  for (unsigned int i=0; i<dofs_per_cell; ++i)
1007  if (start->get_dof_handler().locally_owned_dofs().is_element(local_dof_indices[i]))
1008  block_to_dof_map[block_list[fe_index][i]].
1009  push_back (local_dof_indices[i]);
1010  }
1011 
1012  // now we've got all indices sorted
1013  // into buckets labeled by their
1014  // target block number. we've
1015  // only got to traverse this list
1016  // and assign the new indices
1017  //
1018  // however, we first want to sort
1019  // the indices entered into the
1020  // buckets to preserve the order
1021  // within each component and during
1022  // this also remove duplicate
1023  // entries
1024  for (unsigned int block=0; block<fe_collection.n_blocks();
1025  ++block)
1026  {
1027  std::sort (block_to_dof_map[block].begin(),
1028  block_to_dof_map[block].end());
1029  block_to_dof_map[block]
1030  .erase (std::unique (block_to_dof_map[block].begin(),
1031  block_to_dof_map[block].end()),
1032  block_to_dof_map[block].end());
1033  }
1034 
1035  // calculate the number of locally owned
1036  // DoFs per bucket
1037  const unsigned int n_buckets = fe_collection.n_blocks();
1038  std::vector<types::global_dof_index> shifts(n_buckets);
1039 
1041  = (dynamic_cast<const parallel::Triangulation<dim,spacedim>*>
1042  (&start->get_dof_handler().get_triangulation())))
1043  {
1044 #ifdef DEAL_II_WITH_MPI
1045  std::vector<types::global_dof_index> local_dof_count(n_buckets);
1046 
1047  for (unsigned int c=0; c<n_buckets; ++c)
1048  local_dof_count[c] = block_to_dof_map[c].size();
1049 
1050 
1051  // gather information from all CPUs
1052  std::vector<types::global_dof_index>
1053  all_dof_counts(fe_collection.n_components() *
1054  Utilities::MPI::n_mpi_processes (tria->get_communicator()));
1055 
1056  MPI_Allgather ( &local_dof_count[0],
1057  n_buckets, DEAL_II_DOF_INDEX_MPI_TYPE,
1058  &all_dof_counts[0],
1059  n_buckets, DEAL_II_DOF_INDEX_MPI_TYPE,
1060  tria->get_communicator());
1061 
1062  for (unsigned int i=0; i<n_buckets; ++i)
1063  Assert (all_dof_counts[n_buckets*tria->locally_owned_subdomain()+i]
1064  ==
1065  local_dof_count[i],
1066  ExcInternalError());
1067 
1068  //calculate shifts
1069  types::global_dof_index cumulated = 0;
1070  for (unsigned int c=0; c<n_buckets; ++c)
1071  {
1072  shifts[c]=cumulated;
1073  for (types::subdomain_id i=0; i<tria->locally_owned_subdomain(); ++i)
1074  shifts[c] += all_dof_counts[c+n_buckets*i];
1075  for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes (tria->get_communicator()); ++i)
1076  cumulated += all_dof_counts[c+n_buckets*i];
1077  }
1078 #else
1079  (void)tria;
1080  Assert (false, ExcInternalError());
1081 #endif
1082  }
1083  else
1084  {
1085  shifts[0] = 0;
1086  for (unsigned int c=1; c<fe_collection.n_blocks(); ++c)
1087  shifts[c] = shifts[c-1] + block_to_dof_map[c-1].size();
1088  }
1089 
1090 
1091 
1092 
1093  // now concatenate all the
1094  // components in the order the user
1095  // desired to see
1096  types::global_dof_index next_free_index = 0;
1097  for (unsigned int block=0; block<fe_collection.n_blocks(); ++block)
1098  {
1099  const typename std::vector<types::global_dof_index>::const_iterator
1100  begin_of_component = block_to_dof_map[block].begin(),
1101  end_of_component = block_to_dof_map[block].end();
1102 
1103  next_free_index = shifts[block];
1104 
1105  for (typename std::vector<types::global_dof_index>::const_iterator
1106  dof_index = begin_of_component;
1107  dof_index != end_of_component; ++dof_index)
1108  {
1109  Assert (start->get_dof_handler().locally_owned_dofs()
1110  .index_within_set(*dof_index)
1111  <
1112  new_indices.size(),
1113  ExcInternalError());
1114  new_indices[start->get_dof_handler().locally_owned_dofs()
1115  .index_within_set(*dof_index)]
1116  = next_free_index++;
1117  }
1118  }
1119 
1120  return next_free_index;
1121  }
1122 
1123 
1124 
1125  namespace
1126  {
1127  // helper function for hierarchical()
1128 
1129 // Note that this function only works for active dofs.
1130  template <int dim, class iterator>
1132  compute_hierarchical_recursive (
1133  types::global_dof_index next_free,
1134  std::vector<types::global_dof_index> &new_indices,
1135  const iterator &cell,
1136  const IndexSet &locally_owned)
1137  {
1138  if (cell->has_children())
1139  {
1140  //recursion
1141  for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; ++c)
1142  next_free = compute_hierarchical_recursive<dim> (
1143  next_free,
1144  new_indices,
1145  cell->child (c),
1146  locally_owned);
1147  }
1148  else
1149  {
1150  if (cell->is_locally_owned())
1151  {
1152  const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
1153  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1154  cell->get_dof_indices (local_dof_indices);
1155 
1156  for (unsigned int i = 0; i < dofs_per_cell; ++i)
1157  {
1158  if (locally_owned.is_element (local_dof_indices[i]))
1159  {
1160  // this is a locally owned DoF, assign new number if not assigned a number yet
1161  unsigned int idx = locally_owned.index_within_set (local_dof_indices[i]);
1162  if (new_indices[idx] == DoFHandler<dim>::invalid_dof_index)
1163  {
1164  new_indices[idx] = locally_owned.nth_index_in_set (next_free);
1165  next_free++;
1166  }
1167  }
1168  }
1169  }
1170  }
1171  return next_free;
1172  }
1173  }
1174 
1175 
1176 
1177  template <int dim>
1178  void
1180  {
1181  std::vector<types::global_dof_index> renumbering (dof_handler.n_locally_owned_dofs(),
1183 
1185 
1186  types::global_dof_index next_free = 0;
1187  const IndexSet locally_owned = dof_handler.locally_owned_dofs();
1188 
1190  = dynamic_cast<const parallel::distributed::Triangulation<dim>*>
1191  (&dof_handler.get_triangulation());
1192 
1193  if (tria)
1194  {
1195 #ifdef DEAL_II_WITH_P4EST
1196  // this is a distributed Triangulation. We need to traverse the coarse
1197  // cells in the order p4est does
1198  for (unsigned int c = 0; c < tria->n_cells (0); ++c)
1199  {
1200  const unsigned int coarse_cell_index =
1202 
1204  this_cell (tria, 0, coarse_cell_index, &dof_handler);
1205 
1206  next_free = compute_hierarchical_recursive<dim> (next_free,
1207  renumbering,
1208  this_cell,
1209  locally_owned);
1210  }
1211 #else
1212  Assert (false, ExcNotImplemented());
1213 #endif
1214  }
1215  else
1216  {
1217  //this is not a distributed Triangulation. Traverse coarse cells in the
1218  //normal order
1219  for (cell = dof_handler.begin (0); cell != dof_handler.end (0); ++cell)
1220  next_free = compute_hierarchical_recursive<dim> (next_free,
1221  renumbering,
1222  cell,
1223  locally_owned);
1224  }
1225 
1226  // verify that the last numbered
1227  // degree of freedom is either
1228  // equal to the number of degrees
1229  // of freedom in total (the
1230  // sequential case) or in the
1231  // distributed case at least
1232  // makes sense
1233  Assert ((next_free == dof_handler.n_locally_owned_dofs())
1234  ||
1235  ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs())
1236  &&
1237  (next_free <= dof_handler.n_dofs())),
1238  ExcRenumberingIncomplete());
1239 
1240  // make sure that all local DoFs got new numbers assigned
1241  Assert (std::find (renumbering.begin(), renumbering.end(),
1243  == renumbering.end(),
1244  ExcInternalError());
1245 
1246  dof_handler.renumber_dofs(renumbering);
1247  }
1248 
1249 
1250 
1251  template <typename DoFHandlerType>
1252  void
1253  sort_selected_dofs_back (DoFHandlerType &dof_handler,
1254  const std::vector<bool> &selected_dofs)
1255  {
1256  std::vector<types::global_dof_index> renumbering(dof_handler.n_dofs(),
1257  DoFHandlerType::invalid_dof_index);
1258  compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs);
1259 
1260  dof_handler.renumber_dofs(renumbering);
1261  }
1262 
1263 
1264 
1265  template <typename DoFHandlerType>
1266  void
1267  sort_selected_dofs_back (DoFHandlerType &dof_handler,
1268  const std::vector<bool> &selected_dofs,
1269  const unsigned int level)
1270  {
1271  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
1272  ExcNotInitialized());
1273 
1274  std::vector<types::global_dof_index> renumbering(dof_handler.n_dofs(level),
1275  DoFHandlerType::invalid_dof_index);
1276  compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs, level);
1277 
1278  dof_handler.renumber_dofs(level, renumbering);
1279  }
1280 
1281 
1282 
1283  template <typename DoFHandlerType>
1284  void
1285  compute_sort_selected_dofs_back (std::vector<types::global_dof_index> &new_indices,
1286  const DoFHandlerType &dof_handler,
1287  const std::vector<bool> &selected_dofs)
1288  {
1289  const types::global_dof_index n_dofs = dof_handler.n_dofs();
1290  Assert (selected_dofs.size() == n_dofs,
1291  ExcDimensionMismatch (selected_dofs.size(), n_dofs));
1292 
1293  // re-sort the dofs according to
1294  // their selection state
1295  Assert (new_indices.size() == n_dofs,
1296  ExcDimensionMismatch(new_indices.size(), n_dofs));
1297 
1298  const types::global_dof_index n_selected_dofs = std::count (selected_dofs.begin(),
1299  selected_dofs.end(),
1300  false);
1301 
1302  types::global_dof_index next_unselected = 0;
1303  types::global_dof_index next_selected = n_selected_dofs;
1304  for (types::global_dof_index i=0; i<n_dofs; ++i)
1305  if (selected_dofs[i] == false)
1306  {
1307  new_indices[i] = next_unselected;
1308  ++next_unselected;
1309  }
1310  else
1311  {
1312  new_indices[i] = next_selected;
1313  ++next_selected;
1314  };
1315  Assert (next_unselected == n_selected_dofs, ExcInternalError());
1316  Assert (next_selected == n_dofs, ExcInternalError());
1317  }
1318 
1319 
1320 
1321  template <typename DoFHandlerType>
1322  void
1323  compute_sort_selected_dofs_back (std::vector<types::global_dof_index> &new_indices,
1324  const DoFHandlerType &dof_handler,
1325  const std::vector<bool> &selected_dofs,
1326  const unsigned int level)
1327  {
1328  Assert(dof_handler.n_dofs(level) != numbers::invalid_dof_index,
1329  ExcNotInitialized());
1330 
1331  const unsigned int n_dofs = dof_handler.n_dofs(level);
1332  Assert (selected_dofs.size() == n_dofs,
1333  ExcDimensionMismatch (selected_dofs.size(), n_dofs));
1334 
1335  // re-sort the dofs according to
1336  // their selection state
1337  Assert (new_indices.size() == n_dofs,
1338  ExcDimensionMismatch(new_indices.size(), n_dofs));
1339 
1340  const unsigned int n_selected_dofs = std::count (selected_dofs.begin(),
1341  selected_dofs.end(),
1342  false);
1343 
1344  unsigned int next_unselected = 0;
1345  unsigned int next_selected = n_selected_dofs;
1346  for (unsigned int i=0; i<n_dofs; ++i)
1347  if (selected_dofs[i] == false)
1348  {
1349  new_indices[i] = next_unselected;
1350  ++next_unselected;
1351  }
1352  else
1353  {
1354  new_indices[i] = next_selected;
1355  ++next_selected;
1356  };
1357  Assert (next_unselected == n_selected_dofs, ExcInternalError());
1358  Assert (next_selected == n_dofs, ExcInternalError());
1359  }
1360 
1361 
1362 
1363  template <typename DoFHandlerType>
1364  void
1365  cell_wise (DoFHandlerType &dof,
1366  const std::vector<typename DoFHandlerType::active_cell_iterator> &cells)
1367  {
1368  std::vector<types::global_dof_index> renumbering(dof.n_dofs());
1369  std::vector<types::global_dof_index> reverse(dof.n_dofs());
1370  compute_cell_wise(renumbering, reverse, dof, cells);
1371 
1372  dof.renumber_dofs(renumbering);
1373  }
1374 
1375 
1376  template <typename DoFHandlerType>
1377  void
1379  (std::vector<types::global_dof_index> &new_indices,
1380  std::vector<types::global_dof_index> &reverse,
1381  const DoFHandlerType &dof,
1382  const typename std::vector<typename DoFHandlerType::active_cell_iterator> &cells)
1383  {
1384  Assert(cells.size() == dof.get_triangulation().n_active_cells(),
1385  ExcDimensionMismatch(cells.size(),
1386  dof.get_triangulation().n_active_cells()));
1387 
1388  types::global_dof_index n_global_dofs = dof.n_dofs();
1389 
1390  // Actually, we compute the
1391  // inverse of the reordering
1392  // vector, called reverse here.
1393  // Later, irs inverse is computed
1394  // into new_indices, which is the
1395  // return argument.
1396 
1397  Assert(new_indices.size() == n_global_dofs,
1398  ExcDimensionMismatch(new_indices.size(), n_global_dofs));
1399  Assert(reverse.size() == n_global_dofs,
1400  ExcDimensionMismatch(reverse.size(), n_global_dofs));
1401 
1402  // For continuous elements, we must
1403  // make sure, that each dof is
1404  // reordered only once.
1405  std::vector<bool> already_sorted(n_global_dofs, false);
1406  std::vector<types::global_dof_index> cell_dofs;
1407 
1408  unsigned int global_index = 0;
1409 
1410  typename std::vector<typename DoFHandlerType::active_cell_iterator>::const_iterator cell;
1411 
1412  for (cell = cells.begin(); cell != cells.end(); ++cell)
1413  {
1414  // Determine the number of dofs
1415  // on this cell and reinit the
1416  // vector storing these
1417  // numbers.
1418  unsigned int n_cell_dofs = (*cell)->get_fe().n_dofs_per_cell();
1419  cell_dofs.resize(n_cell_dofs);
1420 
1421  (*cell)->get_active_or_mg_dof_indices(cell_dofs);
1422 
1423  // Sort here to make sure that
1424  // degrees of freedom inside a
1425  // single cell are in the same
1426  // order after renumbering.
1427  std::sort(cell_dofs.begin(), cell_dofs.end());
1428 
1429  for (unsigned int i=0; i<n_cell_dofs; ++i)
1430  {
1431  if (!already_sorted[cell_dofs[i]])
1432  {
1433  already_sorted[cell_dofs[i]] = true;
1434  reverse[global_index++] = cell_dofs[i];
1435  }
1436  }
1437  }
1438  Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
1439 
1440  for (types::global_dof_index i=0; i<reverse.size(); ++i)
1441  new_indices[reverse[i]] = i;
1442  }
1443 
1444 
1445 
1446  template <typename DoFHandlerType>
1447  void cell_wise
1448  (DoFHandlerType &dof,
1449  const unsigned int level,
1450  const typename std::vector<typename DoFHandlerType::level_cell_iterator> &cells)
1451  {
1452  Assert(dof.n_dofs(level) != numbers::invalid_dof_index,
1453  ExcNotInitialized());
1454 
1455  std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1456  std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1457 
1458  compute_cell_wise(renumbering, reverse, dof, level, cells);
1459  dof.renumber_dofs(level, renumbering);
1460  }
1461 
1462 
1463 
1464  template <typename DoFHandlerType>
1465  void compute_cell_wise
1466  (std::vector<types::global_dof_index> &new_order,
1467  std::vector<types::global_dof_index> &reverse,
1468  const DoFHandlerType &dof,
1469  const unsigned int level,
1470  const typename std::vector<typename DoFHandlerType::level_cell_iterator> &cells)
1471  {
1472  Assert(cells.size() == dof.get_triangulation().n_cells(level),
1473  ExcDimensionMismatch(cells.size(),
1474  dof.get_triangulation().n_cells(level)));
1475  Assert (new_order.size() == dof.n_dofs(level),
1476  ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
1477  Assert (reverse.size() == dof.n_dofs(level),
1478  ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
1479 
1480  unsigned int n_global_dofs = dof.n_dofs(level);
1481  unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
1482 
1483  std::vector<bool> already_sorted(n_global_dofs, false);
1484  std::vector<types::global_dof_index> cell_dofs(n_cell_dofs);
1485 
1486  unsigned int global_index = 0;
1487 
1488  typename std::vector<typename DoFHandlerType::level_cell_iterator>::const_iterator cell;
1489 
1490  for (cell = cells.begin(); cell != cells.end(); ++cell)
1491  {
1492  Assert ((*cell)->level() == (int) level, ExcInternalError());
1493 
1494  (*cell)->get_active_or_mg_dof_indices(cell_dofs);
1495  std::sort(cell_dofs.begin(), cell_dofs.end());
1496 
1497  for (unsigned int i=0; i<n_cell_dofs; ++i)
1498  {
1499  if (!already_sorted[cell_dofs[i]])
1500  {
1501  already_sorted[cell_dofs[i]] = true;
1502  reverse[global_index++] = cell_dofs[i];
1503  }
1504  }
1505  }
1506  Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
1507 
1508  for (types::global_dof_index i=0; i<new_order.size(); ++i)
1509  new_order[reverse[i]] = i;
1510  }
1511 
1512 
1513 
1514 
1515 
1516 
1517 
1518  template <typename DoFHandlerType>
1519  void
1521  (std::vector<types::global_dof_index> &new_indices,
1522  std::vector<types::global_dof_index> &reverse,
1523  const DoFHandlerType &dof,
1524  const Point<DoFHandlerType::space_dimension> &direction,
1525  const bool dof_wise_renumbering)
1526  {
1527  if (dof_wise_renumbering == false)
1528  {
1529  std::vector<typename DoFHandlerType::active_cell_iterator> ordered_cells;
1530  ordered_cells.reserve(dof.get_triangulation().n_active_cells());
1531  const CompareDownstream<typename DoFHandlerType::active_cell_iterator,
1532  DoFHandlerType::space_dimension> comparator(direction);
1533 
1534  typename DoFHandlerType::active_cell_iterator p = dof.begin_active();
1535  typename DoFHandlerType::active_cell_iterator end = dof.end();
1536 
1537  while (p!=end)
1538  {
1539  ordered_cells.push_back(p);
1540  ++p;
1541  }
1542  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
1543 
1544  compute_cell_wise(new_indices, reverse, dof, ordered_cells);
1545  }
1546  else
1547  {
1548  // similar code as for
1549  // DoFTools::map_dofs_to_support_points, but
1550  // need to do this for general DoFHandlerType classes and
1551  // want to be able to sort the result
1552  // (otherwise, could use something like
1553  // DoFTools::map_support_points_to_dofs)
1554  const unsigned int n_dofs = dof.n_dofs();
1555  std::vector<std::pair<Point<DoFHandlerType::space_dimension>,unsigned int> > support_point_list
1556  (n_dofs);
1557 
1558  const hp::FECollection<DoFHandlerType::dimension> fe_collection (dof.get_fe ());
1559  Assert (fe_collection[0].has_support_points(),
1561  hp::QCollection<DoFHandlerType::dimension> quadrature_collection;
1562  for (unsigned int comp=0; comp<fe_collection.size(); ++comp)
1563  {
1564  Assert (fe_collection[comp].has_support_points(),
1566  quadrature_collection.push_back
1567  (Quadrature<DoFHandlerType::dimension> (fe_collection[comp].
1568  get_unit_support_points()));
1569  }
1571  hp_fe_values (fe_collection, quadrature_collection,
1573 
1574  std::vector<bool> already_touched (n_dofs, false);
1575 
1576  std::vector<types::global_dof_index> local_dof_indices;
1577  typename DoFHandlerType::active_cell_iterator begin = dof.begin_active();
1578  typename DoFHandlerType::active_cell_iterator end = dof.end();
1579  for ( ; begin != end; ++begin)
1580  {
1581  const unsigned int dofs_per_cell = begin->get_fe().dofs_per_cell;
1582  local_dof_indices.resize (dofs_per_cell);
1583  hp_fe_values.reinit (begin);
1584  const FEValues<DoFHandlerType::dimension> &fe_values =
1585  hp_fe_values.get_present_fe_values ();
1586  begin->get_active_or_mg_dof_indices(local_dof_indices);
1587  const std::vector<Point<DoFHandlerType::space_dimension> > &points
1588  = fe_values.get_quadrature_points ();
1589  for (unsigned int i=0; i<dofs_per_cell; ++i)
1590  if (!already_touched[local_dof_indices[i]])
1591  {
1592  support_point_list[local_dof_indices[i]].first = points[i];
1593  support_point_list[local_dof_indices[i]].second =
1594  local_dof_indices[i];
1595  already_touched[local_dof_indices[i]] = true;
1596  }
1597  }
1598 
1600  std::sort (support_point_list.begin(), support_point_list.end(),
1601  comparator);
1602  for (types::global_dof_index i=0; i<n_dofs; ++i)
1603  new_indices[support_point_list[i].second] = i;
1604  }
1605  }
1606 
1607 
1608 
1609  template <typename DoFHandlerType>
1610  void downstream (DoFHandlerType &dof,
1611  const unsigned int level,
1612  const Point<DoFHandlerType::space_dimension> &direction,
1613  const bool dof_wise_renumbering)
1614  {
1615  std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1616  std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1617  compute_downstream(renumbering, reverse, dof, level, direction,
1618  dof_wise_renumbering);
1619 
1620  dof.renumber_dofs(level, renumbering);
1621  }
1622 
1623 
1624 
1625  template <typename DoFHandlerType>
1626  void
1628  (std::vector<types::global_dof_index> &new_indices,
1629  std::vector<types::global_dof_index> &reverse,
1630  const DoFHandlerType &dof,
1631  const unsigned int level,
1632  const Point<DoFHandlerType::space_dimension> &direction,
1633  const bool dof_wise_renumbering)
1634  {
1635  if (dof_wise_renumbering == false)
1636  {
1637  std::vector<typename DoFHandlerType::level_cell_iterator> ordered_cells;
1638  ordered_cells.reserve (dof.get_triangulation().n_cells(level));
1639  const CompareDownstream<typename DoFHandlerType::level_cell_iterator,
1640  DoFHandlerType::space_dimension> comparator(direction);
1641 
1642  typename DoFHandlerType::level_cell_iterator p = dof.begin(level);
1643  typename DoFHandlerType::level_cell_iterator end = dof.end(level);
1644 
1645  while (p!=end)
1646  {
1647  ordered_cells.push_back(p);
1648  ++p;
1649  }
1650  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
1651 
1652  compute_cell_wise(new_indices, reverse, dof, level, ordered_cells);
1653  }
1654  else
1655  {
1656  Assert (dof.get_fe().has_support_points(),
1658  const unsigned int n_dofs = dof.n_dofs(level);
1659  std::vector<std::pair<Point<DoFHandlerType::space_dimension>,unsigned int> > support_point_list
1660  (n_dofs);
1661 
1662  Quadrature<DoFHandlerType::dimension> q_dummy(dof.get_fe().get_unit_support_points());
1665 
1666  std::vector<bool> already_touched (dof.n_dofs(), false);
1667 
1668  const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
1669  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1670  typename DoFHandlerType::level_cell_iterator begin = dof.begin(level);
1671  typename DoFHandlerType::level_cell_iterator end = dof.end(level);
1672  for ( ; begin != end; ++begin)
1673  {
1674  const typename Triangulation<DoFHandlerType::dimension,
1675  DoFHandlerType::space_dimension>::cell_iterator &begin_tria = begin;
1676  begin->get_active_or_mg_dof_indices(local_dof_indices);
1677  fe_values.reinit (begin_tria);
1678  const std::vector<Point<DoFHandlerType::space_dimension> > &points
1679  = fe_values.get_quadrature_points ();
1680  for (unsigned int i=0; i<dofs_per_cell; ++i)
1681  if (!already_touched[local_dof_indices[i]])
1682  {
1683  support_point_list[local_dof_indices[i]].first = points[i];
1684  support_point_list[local_dof_indices[i]].second =
1685  local_dof_indices[i];
1686  already_touched[local_dof_indices[i]] = true;
1687  }
1688  }
1689 
1691  std::sort (support_point_list.begin(), support_point_list.end(),
1692  comparator);
1693  for (types::global_dof_index i=0; i<n_dofs; ++i)
1694  new_indices[support_point_list[i].second] = i;
1695  }
1696  }
1697 
1698 
1699 
1703  namespace internal
1704  {
1705  template <int dim>
1706  struct ClockCells
1707  {
1711  const Point<dim> &center;
1715  bool counter;
1716 
1720  ClockCells (const Point<dim> &center, bool counter) :
1721  center(center),
1722  counter(counter)
1723  {}
1724 
1728  template <class DHCellIterator>
1729  bool operator () (const DHCellIterator &c1,
1730  const DHCellIterator &c2) const
1731  {
1732  // dispatch to
1733  // dimension-dependent functions
1734  return compare (c1, c2, ::internal::int2type<dim>());
1735  }
1736 
1737  private:
1741  template <class DHCellIterator, int xdim>
1742  bool compare (const DHCellIterator &c1,
1743  const DHCellIterator &c2,
1744  ::internal::int2type<xdim>) const
1745  {
1746  const Tensor<1,dim> v1 = c1->center() - center;
1747  const Tensor<1,dim> v2 = c2->center() - center;
1748  const double s1 = std::atan2(v1[0], v1[1]);
1749  const double s2 = std::atan2(v2[0], v2[1]);
1750  return ( counter ? (s1>s2) : (s2>s1));
1751  }
1752 
1753 
1758  template <class DHCellIterator>
1759  bool compare (const DHCellIterator &,
1760  const DHCellIterator &,
1761  ::internal::int2type<1>) const
1762  {
1763  Assert (dim >= 2,
1764  ExcMessage ("This operation only makes sense for dim>=2."));
1765  return false;
1766  }
1767 
1768  };
1769  }
1770 
1771 
1772 
1773  template <typename DoFHandlerType>
1774  void
1776  DoFHandlerType &dof,
1778  const bool counter)
1779  {
1780  std::vector<types::global_dof_index> renumbering(dof.n_dofs());
1781  compute_clockwise_dg(renumbering, dof, center, counter);
1782 
1783  dof.renumber_dofs(renumbering);
1784  }
1785 
1786 
1787 
1788  template <typename DoFHandlerType>
1789  void
1791  (std::vector<types::global_dof_index> &new_indices,
1792  const DoFHandlerType &dof,
1794  const bool counter)
1795  {
1796  std::vector<typename DoFHandlerType::active_cell_iterator> ordered_cells;
1797  ordered_cells.reserve (dof.get_triangulation().n_active_cells());
1798  internal::ClockCells<DoFHandlerType::space_dimension> comparator(center, counter);
1799 
1800  typename DoFHandlerType::active_cell_iterator p = dof.begin_active();
1801  typename DoFHandlerType::active_cell_iterator end = dof.end();
1802 
1803  while (p!=end)
1804  {
1805  ordered_cells.push_back(p);
1806  ++p;
1807  }
1808  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
1809 
1810  std::vector<types::global_dof_index> reverse(new_indices.size());
1811  compute_cell_wise(new_indices, reverse, dof, ordered_cells);
1812  }
1813 
1814 
1815 
1816  template <typename DoFHandlerType>
1817  void clockwise_dg (DoFHandlerType &dof,
1818  const unsigned int level,
1820  const bool counter)
1821  {
1822  std::vector<typename DoFHandlerType::level_cell_iterator> ordered_cells;
1823  ordered_cells.reserve(dof.get_triangulation().n_active_cells());
1824  internal::ClockCells<DoFHandlerType::space_dimension> comparator(center, counter);
1825 
1826  typename DoFHandlerType::level_cell_iterator p = dof.begin(level);
1827  typename DoFHandlerType::level_cell_iterator end = dof.end(level);
1828 
1829  while (p!=end)
1830  {
1831  ordered_cells.push_back(p);
1832  ++p;
1833  }
1834  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
1835 
1836  cell_wise(dof, level, ordered_cells);
1837  }
1838 
1839 
1840 
1841  template <typename DoFHandlerType>
1842  void
1843  random (DoFHandlerType &dof_handler)
1844  {
1845  std::vector<types::global_dof_index> renumbering(dof_handler.n_dofs(),
1846  DoFHandlerType::invalid_dof_index);
1847  compute_random(renumbering, dof_handler);
1848 
1849  dof_handler.renumber_dofs(renumbering);
1850  }
1851 
1852 
1853 
1854  template <typename DoFHandlerType>
1855  void
1857  std::vector<types::global_dof_index> &new_indices,
1858  const DoFHandlerType &dof_handler)
1859  {
1860  const types::global_dof_index n_dofs = dof_handler.n_dofs();
1861  Assert(new_indices.size() == n_dofs,
1862  ExcDimensionMismatch(new_indices.size(), n_dofs));
1863 
1864  for (unsigned int i=0; i<n_dofs; ++i)
1865  new_indices[i] = i;
1866 
1867  // shuffle the elements; the following is essentially the
1868  // std::random_shuffle algorithm but uses a predictable
1869  // random number generator
1870  ::boost::mt19937 random_number_generator;
1871  for (unsigned int i=1; i<n_dofs; ++i)
1872  {
1873  // get a random number between 0 and i (inclusive)
1874  const unsigned int j
1875  = ::boost::random::uniform_int_distribution<>(0, i)(random_number_generator);
1876 
1877  // if possible, swap the elements
1878  if (i != j)
1879  std::swap (new_indices[i], new_indices[j]);
1880  }
1881  }
1882 
1883 
1884 
1885  template <typename DoFHandlerType>
1886  void
1887  subdomain_wise (DoFHandlerType &dof_handler)
1888  {
1889  std::vector<types::global_dof_index> renumbering(dof_handler.n_dofs(),
1890  DoFHandlerType::invalid_dof_index);
1891  compute_subdomain_wise(renumbering, dof_handler);
1892 
1893  dof_handler.renumber_dofs(renumbering);
1894  }
1895 
1896 
1897 
1898  template <typename DoFHandlerType>
1899  void
1900  compute_subdomain_wise (std::vector<types::global_dof_index> &new_dof_indices,
1901  const DoFHandlerType &dof_handler)
1902  {
1903  const types::global_dof_index n_dofs = dof_handler.n_dofs();
1904  Assert (new_dof_indices.size() == n_dofs,
1905  ExcDimensionMismatch (new_dof_indices.size(), n_dofs));
1906 
1907  // first get the association of each dof
1908  // with a subdomain and determine the total
1909  // number of subdomain ids used
1910  std::vector<types::subdomain_id> subdomain_association (n_dofs);
1912  subdomain_association);
1913  const unsigned int n_subdomains
1914  = *std::max_element (subdomain_association.begin(),
1915  subdomain_association.end()) + 1;
1916 
1917  // then renumber the subdomains by first
1918  // looking at those belonging to subdomain
1919  // 0, then those of subdomain 1, etc. note
1920  // that the algorithm is stable, i.e. if
1921  // two dofs i,j have i<j and belong to the
1922  // same subdomain, then they will be in
1923  // this order also after reordering
1924  std::fill (new_dof_indices.begin(), new_dof_indices.end(),
1926  types::global_dof_index next_free_index = 0;
1927  for (types::subdomain_id subdomain=0; subdomain<n_subdomains; ++subdomain)
1928  for (types::global_dof_index i=0; i<n_dofs; ++i)
1929  if (subdomain_association[i] == subdomain)
1930  {
1931  Assert (new_dof_indices[i] == numbers::invalid_dof_index,
1932  ExcInternalError());
1933  new_dof_indices[i] = next_free_index;
1934  ++next_free_index;
1935  }
1936 
1937  // we should have numbered all dofs
1938  Assert (next_free_index == n_dofs, ExcInternalError());
1939  Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
1941  == new_dof_indices.end(),
1942  ExcInternalError());
1943  }
1944 
1945 } // namespace DoFRenumbering
1946 
1947 
1948 
1949 /*-------------- Explicit Instantiations -------------------------------*/
1950 #include "dof_renumbering.inst"
1951 
1952 
1953 DEAL_II_NAMESPACE_CLOSE
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: dof_handler.cc:845
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
Definition: fe.h:2885
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1052
types::global_dof_index n_dofs() const
void minimum_degree(DoFHandlerType &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
void compute_cell_wise(std::vector< types::global_dof_index > &renumbering, std::vector< types::global_dof_index > &inverse_renumbering, const DoFHandlerType &dof_handler, const std::vector< typename DoFHandlerType::active_cell_iterator > &cell_order)
void compute_king_ordering(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &, const bool reversed_numbering=false, const bool use_constraints=false)
void sort_selected_dofs_back(DoFHandlerType &dof_handler, const std::vector< bool > &selected_dofs)
void get_subdomain_association(const DoFHandlerType &dof_handler, std::vector< types::subdomain_id > &subdomain)
Definition: dof_tools.cc:1252
::ExceptionBase & ExcMessage(std::string arg1)
void hierarchical(DoFHandler< dim > &dof_handler)
size_type n_elements() const
Definition: index_set.h:1380
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity, const unsigned int level)
Definition: mg_tools.cc:539
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 >())
const Triangulation< dim, spacedim > & get_triangulation() const
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:221
STL namespace.
Transformed quadrature points.
size_type size() const
Definition: index_set.h:1247
iterator end()
bool is_primitive(const unsigned int i) const
Definition: fe.h:2937
Definition: point.h:89
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:89
ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:157
active_cell_iterator begin_active(const unsigned int level=0) const
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType, lda > > cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
Definition: fe_values.cc:144
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, ConstraintMatrix &constraints)
void random(DoFHandlerType &dof_handler)
const std::vector< Point< spacedim > > & get_quadrature_points() const
void compute_minimum_degree(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &, const bool reversed_numbering=false, const bool use_constraints=false)
unsigned int n_cells() const
Definition: tria.cc:10966
void Cuthill_McKee(DoFHandlerType &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
unsigned int global_dof_index
Definition: types.h:88
const std::vector< types::global_dof_index > & get_p4est_tree_to_coarse_cell_permutation() const
Definition: tria.cc:4214
#define Assert(cond, exc)
Definition: exceptions.h:294
types::global_dof_index n_dofs() const
void compute_subdomain_wise(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler)
unsigned int n_locally_owned_dofs() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:2915
void compute_random(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler)
iterator begin()
unsigned int subdomain_id
Definition: types.h:42
cell_iterator end() const
Definition: dof_handler.cc:861
void downstream(DoFHandlerType &dof_handler, const Point< DoFHandlerType::space_dimension > &direction, const bool dof_wise_renumbering=false)
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:95
void compute_sort_selected_dofs_back(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler, const std::vector< bool > &selected_dofs)
const unsigned int dofs_per_cell
Definition: fe_base.h:283
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
cell_iterator begin(const unsigned int level=0) const
Definition: dof_handler.cc:833
cell_iterator end() const
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const ConstraintMatrix &constraints=ConstraintMatrix(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void compute_clockwise_dg(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler, const Point< DoFHandlerType::space_dimension > &center, const bool counter)
void king_ordering(DoFHandlerType &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:2742
void block_wise(DoFHandler< dim, spacedim > &dof_handler)
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1475
Definition: fe.h:31
types::global_dof_index compute_component_wise(std::vector< types::global_dof_index > &new_dof_indices, const ITERATOR &start, const ENDITERATOR &end, const std::vector< unsigned int > &target_component, bool is_level_operation)
void compute_Cuthill_McKee(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &, const bool reversed_numbering=false, const bool use_constraints=false)
std::vector< unsigned int > reverse_permutation(const std::vector< unsigned int > &permutation)
Definition: utilities.cc:469
bool is_element(const size_type index) const
Definition: index_set.h:1317
const ::FEValues< dim, spacedim > & get_present_fe_values() const
const FiniteElement< dim, spacedim > & get_fe() const
void subdomain_wise(DoFHandlerType &dof_handler)
void clockwise_dg(DoFHandlerType &dof_handler, const Point< DoFHandlerType::space_dimension > &center, const bool counter=false)
const types::global_dof_index invalid_dof_index
Definition: types.h:178
const IndexSet & locally_owned_dofs() const
size_type nth_index_in_set(const unsigned int local_index) const
Definition: index_set.h:1432
void compute_downstream(std::vector< types::global_dof_index > &new_dof_indices, std::vector< types::global_dof_index > &reverse, const DoFHandlerType &dof_handler, const Point< DoFHandlerType::space_dimension > &direction, const bool dof_wise_renumbering)
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 push_back(const Quadrature< dim > &new_quadrature)
Definition: q_collection.h:219
void cell_wise(DoFHandlerType &dof_handler, const std::vector< typename DoFHandlerType::active_cell_iterator > &cell_order)