16 #include <deal.II/base/quadrature_lib.h> 17 #include <deal.II/base/template_constraints.h> 18 #include <deal.II/base/thread_management.h> 19 #include <deal.II/base/types.h> 20 #include <deal.II/base/utilities.h> 22 #include <deal.II/distributed/tria.h> 24 #include <deal.II/dofs/dof_accessor.h> 25 #include <deal.II/dofs/dof_handler.h> 26 #include <deal.II/dofs/dof_renumbering.h> 27 #include <deal.II/dofs/dof_tools.h> 29 #include <deal.II/fe/fe.h> 31 #include <deal.II/grid/tria.h> 32 #include <deal.II/grid/tria_iterator.h> 34 #include <deal.II/hp/dof_handler.h> 35 #include <deal.II/hp/fe_collection.h> 36 #include <deal.II/hp/fe_values.h> 38 #include <deal.II/lac/affine_constraints.h> 39 #include <deal.II/lac/dynamic_sparsity_pattern.h> 40 #include <deal.II/lac/sparsity_pattern.h> 41 #include <deal.II/lac/sparsity_tools.h> 43 #include <deal.II/multigrid/mg_tools.h> 45 #include <boost/config.hpp> 46 #include <boost/graph/adjacency_list.hpp> 47 #include <boost/graph/bandwidth.hpp> 48 #include <boost/graph/cuthill_mckee_ordering.hpp> 49 #include <boost/graph/king_ordering.hpp> 50 #include <boost/graph/minimum_degree_ordering.hpp> 51 #include <boost/graph/properties.hpp> 52 #include <boost/random.hpp> 53 #include <boost/random/uniform_int_distribution.hpp> 62 DEAL_II_NAMESPACE_OPEN
70 using namespace ::
boost;
73 using Graph = adjacency_list<vecS,
76 property<vertex_color_t,
78 property<vertex_degree_t, int>>>;
79 using Vertex = graph_traits<Graph>::vertex_descriptor;
80 using size_type = graph_traits<Graph>::vertices_size_type;
82 using Pair = std::pair<size_type, size_type>;
88 template <
typename DoFHandlerType>
90 create_graph(
const DoFHandlerType &dof_handler,
91 const bool use_constraints,
92 boosttypes::Graph & graph,
93 boosttypes::property_map<boosttypes::Graph,
94 boosttypes::vertex_degree_t>::type
105 dof_handler.n_dofs());
109 for (
unsigned int row = 0; row < dsp.n_rows(); ++row)
110 for (
unsigned int col = 0; col < dsp.row_length(row); ++col)
111 add_edge(row, dsp.column_number(row, col), graph);
114 boosttypes::graph_traits<boosttypes::Graph>::vertex_iterator ui, ui_end;
116 graph_degree =
get(::boost::vertex_degree, graph);
117 for (::boost::tie(ui, ui_end) = vertices(graph); ui != ui_end; ++ui)
118 graph_degree[*ui] = degree(*ui, graph);
123 template <
typename DoFHandlerType>
126 const bool reversed_numbering,
127 const bool use_constraints)
129 std::vector<types::global_dof_index> renumbering(
139 dof_handler.renumber_dofs(renumbering);
143 template <
typename DoFHandlerType>
146 const DoFHandlerType & dof_handler,
147 const bool reversed_numbering,
148 const bool use_constraints)
150 boosttypes::Graph graph(dof_handler.n_dofs());
151 boosttypes::property_map<boosttypes::Graph,
152 boosttypes::vertex_degree_t>::type graph_degree;
154 internal::create_graph(dof_handler, use_constraints, graph, graph_degree);
156 boosttypes::property_map<boosttypes::Graph,
157 boosttypes::vertex_index_t>::type index_map =
158 get(::boost::vertex_index, graph);
161 std::vector<boosttypes::Vertex> inv_perm(num_vertices(graph));
163 if (reversed_numbering ==
false)
164 ::boost::cuthill_mckee_ordering(graph,
166 get(::boost::vertex_color, graph),
167 make_degree_map(graph));
169 ::boost::cuthill_mckee_ordering(graph,
171 get(::boost::vertex_color, graph),
172 make_degree_map(graph));
174 for (boosttypes::size_type c = 0; c != inv_perm.size(); ++c)
175 new_dof_indices[index_map[inv_perm[c]]] = c;
177 Assert(std::find(new_dof_indices.begin(),
178 new_dof_indices.end(),
185 template <
typename DoFHandlerType>
188 const bool reversed_numbering,
189 const bool use_constraints)
191 std::vector<types::global_dof_index> renumbering(
201 dof_handler.renumber_dofs(renumbering);
205 template <
typename DoFHandlerType>
208 const DoFHandlerType & dof_handler,
209 const bool reversed_numbering,
210 const bool use_constraints)
212 boosttypes::Graph graph(dof_handler.n_dofs());
213 boosttypes::property_map<boosttypes::Graph,
214 boosttypes::vertex_degree_t>::type graph_degree;
216 internal::create_graph(dof_handler, use_constraints, graph, graph_degree);
218 boosttypes::property_map<boosttypes::Graph,
219 boosttypes::vertex_index_t>::type index_map =
220 get(::boost::vertex_index, graph);
223 std::vector<boosttypes::Vertex> inv_perm(num_vertices(graph));
225 if (reversed_numbering ==
false)
226 ::boost::king_ordering(graph, inv_perm.rbegin());
228 ::boost::king_ordering(graph, inv_perm.begin());
230 for (boosttypes::size_type c = 0; c != inv_perm.size(); ++c)
231 new_dof_indices[index_map[inv_perm[c]]] = c;
233 Assert(std::find(new_dof_indices.begin(),
234 new_dof_indices.end(),
241 template <
typename DoFHandlerType>
244 const bool reversed_numbering,
245 const bool use_constraints)
247 std::vector<types::global_dof_index> renumbering(
257 dof_handler.renumber_dofs(renumbering);
261 template <
typename DoFHandlerType>
264 std::vector<types::global_dof_index> &new_dof_indices,
265 const DoFHandlerType & dof_handler,
266 const bool reversed_numbering,
267 const bool use_constraints)
269 (void)use_constraints;
278 using namespace ::
boost;
283 using Graph = adjacency_list<vecS, vecS, directedS>;
285 int n = dof_handler.n_dofs();
289 std::vector<::types::global_dof_index> dofs_on_this_cell;
291 typename DoFHandlerType::active_cell_iterator cell = dof_handler
293 endc = dof_handler.end();
295 for (; cell != endc; ++cell)
297 const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
299 dofs_on_this_cell.resize(dofs_per_cell);
301 cell->get_active_or_mg_dof_indices(dofs_on_this_cell);
302 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
303 for (
unsigned int j = 0; j < dofs_per_cell; ++j)
304 if (dofs_on_this_cell[i] > dofs_on_this_cell[j])
306 add_edge(dofs_on_this_cell[i], dofs_on_this_cell[j], G);
307 add_edge(dofs_on_this_cell[j], dofs_on_this_cell[i], G);
312 using Vector = std::vector<int>;
315 Vector inverse_perm(n, 0);
320 Vector supernode_sizes(n, 1);
323 ::boost::property_map<Graph, vertex_index_t>::type
id =
324 get(vertex_index, G);
330 minimum_degree_ordering(
332 make_iterator_property_map(degree.begin(), id, degree[0]),
335 make_iterator_property_map(supernode_sizes.begin(),
342 for (
int i = 0; i < n; ++i)
344 Assert(std::find(perm.begin(), perm.end(), i) != perm.end(),
346 Assert(std::find(inverse_perm.begin(), inverse_perm.end(), i) !=
352 if (reversed_numbering ==
true)
353 std::copy(perm.begin(), perm.end(), new_dof_indices.begin());
355 std::copy(inverse_perm.begin(),
357 new_dof_indices.begin());
364 template <
typename DoFHandlerType>
367 const bool reversed_numbering,
368 const bool use_constraints,
369 const std::vector<types::global_dof_index> &starting_indices)
371 std::vector<types::global_dof_index> renumbering(
372 dof_handler.locally_owned_dofs().n_elements(),
383 dof_handler.renumber_dofs(renumbering);
388 template <
typename DoFHandlerType>
391 std::vector<types::global_dof_index> & new_indices,
392 const DoFHandlerType & dof_handler,
393 const bool reversed_numbering,
394 const bool use_constraints,
395 const std::vector<types::global_dof_index> &starting_indices)
399 if (dof_handler.locally_owned_dofs().n_elements() == 0)
415 constraints.
reinit(locally_relevant_dofs);
420 const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
433 if (reversed_numbering)
449 bool needs_locally_active =
false;
450 for (
unsigned int i = 0; i < starting_indices.size(); ++i)
452 if ((needs_locally_active ==
454 (locally_owned_dofs.
is_element(starting_indices[i]) ==
false))
457 locally_active_dofs.
is_element(starting_indices[i]),
459 "You specified global degree of freedom " +
461 " as a starting index, but this index is not among the " 462 "locally active ones on this processor, as required " 463 "for this function."));
464 needs_locally_active =
true;
469 (needs_locally_active ? locally_active_dofs : locally_owned_dofs);
475 dof_handler.n_dofs(),
481 std::vector<types::global_dof_index> row_entries;
482 for (
unsigned int i = 0; i < index_set_to_use.
n_elements(); ++i)
486 const unsigned int row_length = dsp.row_length(row);
488 for (
unsigned int j = 0; j < row_length; ++j)
490 const unsigned int col = dsp.column_number(row, j);
491 if (col != row && index_set_to_use.
is_element(col))
494 local_sparsity.add_entries(i,
501 std::vector<types::global_dof_index> local_starting_indices(
502 starting_indices.size());
503 for (
unsigned int i = 0; i < starting_indices.size(); ++i)
504 local_starting_indices[i] =
509 std::vector<types::global_dof_index> my_new_indices(
513 local_starting_indices);
514 if (reversed_numbering)
523 if (needs_locally_active ==
true)
526 IndexSet active_but_not_owned_dofs = locally_active_dofs;
527 active_but_not_owned_dofs.
subtract_set(locally_owned_dofs);
529 std::set<types::global_dof_index> erase_these_indices;
530 for (
const auto p : active_but_not_owned_dofs)
535 erase_these_indices.insert(my_new_indices[index]);
538 Assert(erase_these_indices.size() ==
539 active_but_not_owned_dofs.n_elements(),
541 Assert(static_cast<unsigned int>(
542 std::count(my_new_indices.begin(),
543 my_new_indices.end(),
545 active_but_not_owned_dofs.n_elements(),
549 std::vector<types::global_dof_index> translate_indices(
550 my_new_indices.size());
552 std::set<types::global_dof_index>::const_iterator
553 next_erased_index = erase_these_indices.begin();
555 for (
unsigned int i = 0; i < translate_indices.size(); ++i)
556 if ((next_erased_index != erase_these_indices.end()) &&
557 (*next_erased_index == i))
564 translate_indices[i] = next_new_index;
574 new_indices.reserve(locally_owned_dofs.
n_elements());
575 for (
const auto &p : my_new_indices)
580 new_indices.push_back(translate_indices[p]);
586 new_indices = std::move(my_new_indices);
592 for (std::size_t i = 0; i < new_indices.size(); ++i)
599 template <
typename DoFHandlerType>
602 const unsigned int level,
603 const bool reversed_numbering,
604 const std::vector<types::global_dof_index> &starting_indices)
611 dof_handler.n_dofs(level));
614 std::vector<types::global_dof_index> new_indices(dsp.n_rows());
617 if (reversed_numbering)
623 dof_handler.renumber_dofs(level, new_indices);
628 template <
typename DoFHandlerType>
631 const std::vector<unsigned int> &component_order_arg)
633 std::vector<types::global_dof_index> renumbering(
636 typename DoFHandlerType::active_cell_iterator start =
637 dof_handler.begin_active();
638 const typename DoFHandlerType::level_cell_iterator end = dof_handler.end();
642 DoFHandlerType::space_dimension>(
643 renumbering, start, end, component_order_arg,
false);
654 Assert((result == dof_handler.n_locally_owned_dofs()) ||
655 ((dof_handler.n_locally_owned_dofs() < dof_handler.n_dofs()) &&
656 (result <= dof_handler.n_dofs())),
659 dof_handler.renumber_dofs(renumbering);
664 template <
typename DoFHandlerType>
667 const unsigned int level,
668 const std::vector<unsigned int> &component_order_arg)
673 std::vector<types::global_dof_index> renumbering(
676 typename DoFHandlerType::level_cell_iterator start =
677 dof_handler.begin(level);
678 typename DoFHandlerType::level_cell_iterator end = dof_handler.end(level);
682 DoFHandlerType::space_dimension>(
683 renumbering, start, end, component_order_arg,
true);
690 if (renumbering.size() != 0)
691 dof_handler.renumber_dofs(level, renumbering);
696 template <
int dim,
int spacedim,
typename CellIterator>
699 const CellIterator & start,
700 const typename identity<CellIterator>::type &end,
701 const std::vector<unsigned int> &component_order_arg,
702 const bool is_level_operation)
705 start->get_dof_handler().get_fe_collection());
709 if (fe_collection.n_components() == 1)
711 new_indices.resize(0);
717 std::vector<unsigned int> component_order(component_order_arg);
722 if (component_order.size() == 0)
723 for (
unsigned int i = 0; i < fe_collection.n_components(); ++i)
724 component_order.push_back(i);
726 Assert(component_order.size() == fe_collection.n_components(),
728 fe_collection.n_components()));
730 for (
unsigned int i = 0; i < component_order.size(); ++i)
731 Assert(component_order[i] < fe_collection.n_components(),
734 fe_collection.n_components()));
738 std::vector<types::global_dof_index> local_dof_indices;
750 std::vector<std::vector<unsigned int>> component_list(fe_collection.size());
751 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
755 component_list[f].resize(dofs_per_cell);
756 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
758 component_list[f][i] =
762 const unsigned int comp =
768 component_list[f][i] = component_order[comp];
783 std::vector<std::vector<types::global_dof_index>> component_to_dof_map(
784 fe_collection.n_components());
785 for (CellIterator cell = start; cell != end; ++cell)
787 if (is_level_operation)
790 if (!cell->is_locally_owned_on_level())
797 if (!cell->is_locally_owned())
803 const unsigned int fe_index = cell->active_fe_index();
804 const unsigned int dofs_per_cell =
805 fe_collection[fe_index].dofs_per_cell;
806 local_dof_indices.resize(dofs_per_cell);
807 cell->get_active_or_mg_dof_indices(local_dof_indices);
808 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
809 if (start->get_dof_handler().locally_owned_dofs().is_element(
810 local_dof_indices[i]))
811 component_to_dof_map[component_list[fe_index][i]].push_back(
812 local_dof_indices[i]);
838 for (
unsigned int component = 0; component < fe_collection.n_components();
841 std::sort(component_to_dof_map[component].begin(),
842 component_to_dof_map[component].end());
843 component_to_dof_map[component].erase(
844 std::unique(component_to_dof_map[component].begin(),
845 component_to_dof_map[component].end()),
846 component_to_dof_map[component].end());
851 const unsigned int n_buckets = fe_collection.n_components();
852 std::vector<types::global_dof_index> shifts(n_buckets);
856 &start->get_dof_handler().get_triangulation())))
858 #ifdef DEAL_II_WITH_MPI 859 std::vector<types::global_dof_index> local_dof_count(n_buckets);
861 for (
unsigned int c = 0; c < n_buckets; ++c)
862 local_dof_count[c] = component_to_dof_map[c].size();
866 std::vector<types::global_dof_index> all_dof_counts(
867 fe_collection.n_components() *
870 const int ierr = MPI_Allgather(local_dof_count.data(),
872 DEAL_II_DOF_INDEX_MPI_TYPE,
873 all_dof_counts.data(),
875 DEAL_II_DOF_INDEX_MPI_TYPE,
876 tria->get_communicator());
879 for (
unsigned int i = 0; i < n_buckets; ++i)
881 all_dof_counts[n_buckets * tria->locally_owned_subdomain() + i] ==
886 unsigned int cumulated = 0;
887 for (
unsigned int c = 0; c < n_buckets; ++c)
889 shifts[c] = cumulated;
892 shifts[c] += all_dof_counts[c + n_buckets * i];
893 for (
unsigned int i = 0;
896 cumulated += all_dof_counts[c + n_buckets * i];
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();
916 for (
unsigned int component = 0; component < fe_collection.n_components();
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();
923 next_free_index = shifts[component];
925 for (
typename std::vector<types::global_dof_index>::const_iterator
926 dof_index = begin_of_component;
927 dof_index != end_of_component;
931 start->get_dof_handler().locally_owned_dofs().index_within_set(
932 *dof_index) < new_indices.size(),
934 new_indices[start->get_dof_handler()
935 .locally_owned_dofs()
936 .index_within_set(*dof_index)] = next_free_index++;
940 return next_free_index;
945 template <
int dim,
int spacedim>
949 std::vector<types::global_dof_index> renumbering(
954 const typename DoFHandler<dim, spacedim>::level_cell_iterator end =
961 typename DoFHandler<dim, spacedim>::level_cell_iterator>(renumbering,
977 (result <= dof_handler.
n_dofs())),
985 template <
int dim,
int spacedim>
989 std::vector<types::global_dof_index> renumbering(
994 const typename hp::DoFHandler<dim, spacedim>::level_cell_iterator end =
1001 typename hp::DoFHandler<dim, spacedim>::level_cell_iterator>(renumbering,
1016 template <
int dim,
int spacedim>
1023 std::vector<types::global_dof_index> renumbering(
1026 typename DoFHandler<dim, spacedim>::level_cell_iterator start =
1027 dof_handler.
begin(level);
1028 typename DoFHandler<dim, spacedim>::level_cell_iterator end =
1029 dof_handler.
end(level);
1034 typename DoFHandler<dim, spacedim>::level_cell_iterator,
1035 typename DoFHandler<dim, spacedim>::level_cell_iterator>(renumbering,
1045 if (renumbering.size() != 0)
1051 template <
int dim,
int spacedim,
class ITERATOR,
class ENDITERATOR>
1054 const ITERATOR & start,
1055 const ENDITERATOR & end,
1056 const bool is_level_operation)
1059 start->get_dof_handler().get_fe_collection());
1063 if (fe_collection.n_blocks() == 1)
1065 new_indices.resize(0);
1071 std::vector<types::global_dof_index> local_dof_indices;
1076 std::vector<std::vector<types::global_dof_index>> block_list(
1077 fe_collection.size());
1078 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
1097 std::vector<std::vector<types::global_dof_index>> block_to_dof_map(
1098 fe_collection.n_blocks());
1099 for (ITERATOR cell = start; cell != end; ++cell)
1101 if (is_level_operation)
1104 if (!cell->is_locally_owned_on_level())
1111 if (!cell->is_locally_owned())
1118 const unsigned int fe_index = cell->active_fe_index();
1119 const unsigned int dofs_per_cell =
1120 fe_collection[fe_index].dofs_per_cell;
1121 local_dof_indices.resize(dofs_per_cell);
1122 cell->get_active_or_mg_dof_indices(local_dof_indices);
1123 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1124 if (start->get_dof_handler().locally_owned_dofs().is_element(
1125 local_dof_indices[i]))
1126 block_to_dof_map[block_list[fe_index][i]].push_back(
1127 local_dof_indices[i]);
1142 for (
unsigned int block = 0; block < fe_collection.n_blocks(); ++block)
1144 std::sort(block_to_dof_map[block].begin(),
1145 block_to_dof_map[block].end());
1146 block_to_dof_map[block].erase(
1147 std::unique(block_to_dof_map[block].begin(),
1148 block_to_dof_map[block].end()),
1149 block_to_dof_map[block].end());
1154 const unsigned int n_buckets = fe_collection.n_blocks();
1155 std::vector<types::global_dof_index> shifts(n_buckets);
1159 &start->get_dof_handler().get_triangulation())))
1161 #ifdef DEAL_II_WITH_MPI 1162 std::vector<types::global_dof_index> local_dof_count(n_buckets);
1164 for (
unsigned int c = 0; c < n_buckets; ++c)
1165 local_dof_count[c] = block_to_dof_map[c].size();
1169 std::vector<types::global_dof_index> all_dof_counts(
1170 fe_collection.n_components() *
1173 const int ierr = MPI_Allgather(local_dof_count.data(),
1175 DEAL_II_DOF_INDEX_MPI_TYPE,
1176 all_dof_counts.data(),
1178 DEAL_II_DOF_INDEX_MPI_TYPE,
1179 tria->get_communicator());
1182 for (
unsigned int i = 0; i < n_buckets; ++i)
1184 all_dof_counts[n_buckets * tria->locally_owned_subdomain() + i] ==
1190 for (
unsigned int c = 0; c < n_buckets; ++c)
1192 shifts[c] = cumulated;
1195 shifts[c] += all_dof_counts[c + n_buckets * i];
1196 for (
unsigned int i = 0;
1199 cumulated += all_dof_counts[c + n_buckets * i];
1209 for (
unsigned int c = 1; c < fe_collection.n_blocks(); ++c)
1210 shifts[c] = shifts[c - 1] + block_to_dof_map[c - 1].size();
1219 for (
unsigned int block = 0; block < fe_collection.n_blocks(); ++block)
1221 const typename std::vector<types::global_dof_index>::const_iterator
1222 begin_of_component = block_to_dof_map[block].begin(),
1223 end_of_component = block_to_dof_map[block].end();
1225 next_free_index = shifts[block];
1227 for (
typename std::vector<types::global_dof_index>::const_iterator
1228 dof_index = begin_of_component;
1229 dof_index != end_of_component;
1233 start->get_dof_handler().locally_owned_dofs().index_within_set(
1234 *dof_index) < new_indices.size(),
1236 new_indices[start->get_dof_handler()
1237 .locally_owned_dofs()
1238 .index_within_set(*dof_index)] = next_free_index++;
1242 return next_free_index;
1254 template <
int dim,
class CellIteratorType>
1256 compute_hierarchical_recursive(
1259 const CellIteratorType & cell,
1260 const IndexSet & locally_owned_dof_indices,
1261 std::vector<types::global_dof_index> &new_indices)
1264 next_free_dof_offset;
1266 if (cell->has_children())
1269 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
1271 current_next_free_dof_offset =
1272 compute_hierarchical_recursive<dim>(current_next_free_dof_offset,
1275 locally_owned_dof_indices,
1295 if (cell->is_locally_owned())
1298 const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
1299 std::vector<types::global_dof_index> local_dof_indices(
1301 cell->get_dof_indices(local_dof_indices);
1320 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1321 if (locally_owned_dof_indices.
is_element(local_dof_indices[i]))
1325 const unsigned int idx =
1327 local_dof_indices[i]);
1331 my_starting_index + current_next_free_dof_offset;
1332 ++current_next_free_dof_offset;
1338 return current_next_free_dof_offset;
1348 std::vector<types::global_dof_index> renumbering(
1371 const std::vector<types::global_dof_index>
1372 &n_locally_owned_dofs_per_processor =
1375 std::accumulate(n_locally_owned_dofs_per_processor.begin(),
1376 n_locally_owned_dofs_per_processor.begin() +
1377 tria->locally_owned_subdomain(),
1385 #ifdef DEAL_II_WITH_P4EST 1390 for (
unsigned int c = 0; c < tria->n_cells(0); ++c)
1392 const unsigned int coarse_cell_index =
1393 tria->get_p4est_tree_to_coarse_cell_permutation()[c];
1396 tria, 0, coarse_cell_index, &dof_handler);
1398 next_free_dof_offset =
1399 compute_hierarchical_recursive<dim>(next_free_dof_offset,
1414 dof_handler.
begin(0);
1415 cell != dof_handler.
end(0);
1417 next_free_dof_offset =
1418 compute_hierarchical_recursive<dim>(next_free_dof_offset,
1431 (next_free_dof_offset <= dof_handler.
n_dofs())),
1435 Assert(std::find(renumbering.begin(),
1445 template <
typename DoFHandlerType>
1448 const std::vector<bool> &selected_dofs)
1450 std::vector<types::global_dof_index> renumbering(
1454 dof_handler.renumber_dofs(renumbering);
1459 template <
typename DoFHandlerType>
1462 const std::vector<bool> &selected_dofs,
1463 const unsigned int level)
1468 std::vector<types::global_dof_index> renumbering(
1475 dof_handler.renumber_dofs(level, renumbering);
1480 template <
typename DoFHandlerType>
1483 std::vector<types::global_dof_index> &new_indices,
1484 const DoFHandlerType & dof_handler,
1485 const std::vector<bool> & selected_dofs)
1488 Assert(selected_dofs.size() == n_dofs,
1493 Assert(new_indices.size() == n_dofs,
1497 std::count(selected_dofs.begin(), selected_dofs.end(),
false);
1502 if (selected_dofs[i] ==
false)
1504 new_indices[i] = next_unselected;
1509 new_indices[i] = next_selected;
1518 template <
typename DoFHandlerType>
1521 std::vector<types::global_dof_index> &new_indices,
1522 const DoFHandlerType & dof_handler,
1523 const std::vector<bool> & selected_dofs,
1524 const unsigned int level)
1529 const unsigned int n_dofs = dof_handler.n_dofs(level);
1530 Assert(selected_dofs.size() == n_dofs,
1535 Assert(new_indices.size() == n_dofs,
1538 const unsigned int n_selected_dofs =
1539 std::count(selected_dofs.begin(), selected_dofs.end(),
false);
1541 unsigned int next_unselected = 0;
1542 unsigned int next_selected = n_selected_dofs;
1543 for (
unsigned int i = 0; i < n_dofs; ++i)
1544 if (selected_dofs[i] ==
false)
1546 new_indices[i] = next_unselected;
1551 new_indices[i] = next_selected;
1560 template <
typename DoFHandlerType>
1563 DoFHandlerType & dof,
1564 const std::vector<typename DoFHandlerType::active_cell_iterator> &cells)
1566 std::vector<types::global_dof_index> renumbering(dof.n_dofs());
1567 std::vector<types::global_dof_index> reverse(dof.n_dofs());
1570 dof.renumber_dofs(renumbering);
1574 template <
typename DoFHandlerType>
1577 std::vector<types::global_dof_index> &new_indices,
1578 std::vector<types::global_dof_index> &reverse,
1579 const DoFHandlerType & dof,
1580 const typename std::vector<typename DoFHandlerType::active_cell_iterator>
1583 Assert(cells.size() == dof.get_triangulation().n_active_cells(),
1585 dof.get_triangulation().n_active_cells()));
1596 Assert(new_indices.size() == n_global_dofs,
1598 Assert(reverse.size() == n_global_dofs,
1604 std::vector<bool> already_sorted(n_global_dofs,
false);
1605 std::vector<types::global_dof_index> cell_dofs;
1607 unsigned int global_index = 0;
1609 typename std::vector<
1610 typename DoFHandlerType::active_cell_iterator>::const_iterator cell;
1612 for (cell = cells.begin(); cell != cells.end(); ++cell)
1618 unsigned int n_cell_dofs = (*cell)->get_fe().n_dofs_per_cell();
1619 cell_dofs.resize(n_cell_dofs);
1621 (*cell)->get_active_or_mg_dof_indices(cell_dofs);
1627 std::sort(cell_dofs.begin(), cell_dofs.end());
1629 for (
unsigned int i = 0; i < n_cell_dofs; ++i)
1631 if (!already_sorted[cell_dofs[i]])
1633 already_sorted[cell_dofs[i]] =
true;
1634 reverse[global_index++] = cell_dofs[i];
1638 Assert(global_index == n_global_dofs,
1640 "Traversing over the given set of cells did not cover all " 1641 "degrees of freedom in the DoFHandler. Does the set of cells " 1642 "not include all active cells?"));
1645 new_indices[reverse[i]] = i;
1650 template <
typename DoFHandlerType>
1653 DoFHandlerType & dof,
1654 const unsigned int level,
1655 const typename std::vector<typename DoFHandlerType::level_cell_iterator>
1661 std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1662 std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1665 dof.renumber_dofs(level, renumbering);
1670 template <
typename DoFHandlerType>
1673 std::vector<types::global_dof_index> &new_order,
1674 std::vector<types::global_dof_index> &reverse,
1675 const DoFHandlerType & dof,
1676 const unsigned int level,
1677 const typename std::vector<typename DoFHandlerType::level_cell_iterator>
1680 Assert(cells.size() == dof.get_triangulation().n_cells(level),
1682 dof.get_triangulation().n_cells(level)));
1683 Assert(new_order.size() == dof.n_dofs(level),
1685 Assert(reverse.size() == dof.n_dofs(level),
1688 unsigned int n_global_dofs = dof.n_dofs(level);
1689 unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
1691 std::vector<bool> already_sorted(n_global_dofs,
false);
1692 std::vector<types::global_dof_index> cell_dofs(n_cell_dofs);
1694 unsigned int global_index = 0;
1696 typename std::vector<
1697 typename DoFHandlerType::level_cell_iterator>::const_iterator cell;
1699 for (cell = cells.begin(); cell != cells.end(); ++cell)
1703 (*cell)->get_active_or_mg_dof_indices(cell_dofs);
1704 std::sort(cell_dofs.begin(), cell_dofs.end());
1706 for (
unsigned int i = 0; i < n_cell_dofs; ++i)
1708 if (!already_sorted[cell_dofs[i]])
1710 already_sorted[cell_dofs[i]] =
true;
1711 reverse[global_index++] = cell_dofs[i];
1715 Assert(global_index == n_global_dofs,
1717 "Traversing over the given set of cells did not cover all " 1718 "degrees of freedom in the DoFHandler. Does the set of cells " 1719 "not include all cells of the specified level?"));
1722 new_order[reverse[i]] = i;
1727 template <
typename DoFHandlerType>
1731 const bool dof_wise_renumbering)
1733 std::vector<types::global_dof_index> renumbering(dof.n_dofs());
1734 std::vector<types::global_dof_index> reverse(dof.n_dofs());
1736 renumbering, reverse, dof, direction, dof_wise_renumbering);
1738 dof.renumber_dofs(renumbering);
1743 template <
typename DoFHandlerType>
1746 std::vector<types::global_dof_index> & new_indices,
1747 std::vector<types::global_dof_index> & reverse,
1748 const DoFHandlerType & dof,
1750 const bool dof_wise_renumbering)
1754 DoFHandlerType::space_dimension
> *>(
1755 &dof.get_triangulation()) ==
nullptr),
1758 if (dof_wise_renumbering ==
false)
1760 std::vector<typename DoFHandlerType::active_cell_iterator>
1762 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
1764 DoFHandlerType::space_dimension>
1765 comparator(direction);
1767 typename DoFHandlerType::active_cell_iterator p = dof.begin_active();
1768 typename DoFHandlerType::active_cell_iterator end = dof.end();
1772 ordered_cells.push_back(p);
1775 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
1787 const unsigned int n_dofs = dof.n_dofs();
1789 std::pair<Point<DoFHandlerType::space_dimension>,
unsigned int>>
1790 support_point_list(n_dofs);
1793 dof.get_fe_collection();
1794 Assert(fe_collection[0].has_support_points(),
1796 DoFHandlerType::dimension>::ExcFEHasNoSupportPoints());
1798 for (
unsigned int comp = 0; comp < fe_collection.
size(); ++comp)
1800 Assert(fe_collection[comp].has_support_points(),
1802 DoFHandlerType::dimension>::ExcFEHasNoSupportPoints());
1805 fe_collection[comp].get_unit_support_points()));
1808 hp_fe_values(fe_collection,
1809 quadrature_collection,
1812 std::vector<bool> already_touched(n_dofs,
false);
1814 std::vector<types::global_dof_index> local_dof_indices;
1815 typename DoFHandlerType::active_cell_iterator begin =
1817 typename DoFHandlerType::active_cell_iterator end = dof.end();
1818 for (; begin != end; ++begin)
1820 const unsigned int dofs_per_cell = begin->get_fe().dofs_per_cell;
1821 local_dof_indices.resize(dofs_per_cell);
1822 hp_fe_values.
reinit(begin);
1825 begin->get_active_or_mg_dof_indices(local_dof_indices);
1826 const std::vector<Point<DoFHandlerType::space_dimension>> &points =
1828 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1829 if (!already_touched[local_dof_indices[i]])
1831 support_point_list[local_dof_indices[i]].first = points[i];
1832 support_point_list[local_dof_indices[i]].second =
1833 local_dof_indices[i];
1834 already_touched[local_dof_indices[i]] =
true;
1840 std::sort(support_point_list.begin(),
1841 support_point_list.end(),
1844 new_indices[support_point_list[i].second] = i;
1850 template <
typename DoFHandlerType>
1853 const unsigned int level,
1855 const bool dof_wise_renumbering)
1857 std::vector<types::global_dof_index> renumbering(dof.n_dofs(level));
1858 std::vector<types::global_dof_index> reverse(dof.n_dofs(level));
1860 renumbering, reverse, dof, level, direction, dof_wise_renumbering);
1862 dof.renumber_dofs(level, renumbering);
1867 template <
typename DoFHandlerType>
1870 std::vector<types::global_dof_index> & new_indices,
1871 std::vector<types::global_dof_index> & reverse,
1872 const DoFHandlerType & dof,
1873 const unsigned int level,
1875 const bool dof_wise_renumbering)
1877 if (dof_wise_renumbering ==
false)
1879 std::vector<typename DoFHandlerType::level_cell_iterator> ordered_cells;
1880 ordered_cells.reserve(dof.get_triangulation().n_cells(level));
1882 DoFHandlerType::space_dimension>
1883 comparator(direction);
1885 typename DoFHandlerType::level_cell_iterator p = dof.begin(level);
1886 typename DoFHandlerType::level_cell_iterator end = dof.end(level);
1890 ordered_cells.push_back(p);
1893 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
1899 Assert(dof.get_fe().has_support_points(),
1901 DoFHandlerType::dimension>::ExcFEHasNoSupportPoints());
1902 const unsigned int n_dofs = dof.n_dofs(level);
1904 std::pair<Point<DoFHandlerType::space_dimension>,
unsigned int>>
1905 support_point_list(n_dofs);
1908 dof.get_fe().get_unit_support_points());
1912 std::vector<bool> already_touched(dof.n_dofs(),
false);
1914 const unsigned int dofs_per_cell = dof.
get_fe().dofs_per_cell;
1915 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1916 typename DoFHandlerType::level_cell_iterator begin = dof.begin(level);
1917 typename DoFHandlerType::level_cell_iterator end = dof.end(level);
1918 for (; begin != end; ++begin)
1921 DoFHandlerType::dimension,
1922 DoFHandlerType::space_dimension>::cell_iterator &begin_tria =
1924 begin->get_active_or_mg_dof_indices(local_dof_indices);
1925 fe_values.reinit(begin_tria);
1926 const std::vector<Point<DoFHandlerType::space_dimension>> &points =
1927 fe_values.get_quadrature_points();
1928 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1929 if (!already_touched[local_dof_indices[i]])
1931 support_point_list[local_dof_indices[i]].first = points[i];
1932 support_point_list[local_dof_indices[i]].second =
1933 local_dof_indices[i];
1934 already_touched[local_dof_indices[i]] =
true;
1940 std::sort(support_point_list.begin(),
1941 support_point_list.end(),
1944 new_indices[support_point_list[i].second] = i;
1970 ClockCells(
const Point<dim> ¢er,
bool counter)
1978 template <
class DHCellIterator>
1980 operator()(
const DHCellIterator &c1,
const DHCellIterator &c2)
const 1984 return compare(c1, c2, std::integral_constant<int, dim>());
1991 template <
class DHCellIterator,
int xdim>
1993 compare(
const DHCellIterator &c1,
1994 const DHCellIterator &c2,
1995 std::integral_constant<int, xdim>)
const 1999 const double s1 = std::atan2(v1[0], v1[1]);
2000 const double s2 = std::atan2(v2[0], v2[1]);
2001 return (counter ? (s1 > s2) : (s2 > s1));
2009 template <
class DHCellIterator>
2011 compare(
const DHCellIterator &,
2012 const DHCellIterator &,
2013 std::integral_constant<int, 1>)
const 2016 ExcMessage(
"This operation only makes sense for dim>=2."));
2024 template <
typename DoFHandlerType>
2030 std::vector<types::global_dof_index> renumbering(dof.n_dofs());
2033 dof.renumber_dofs(renumbering);
2038 template <
typename DoFHandlerType>
2041 const DoFHandlerType & dof,
2045 std::vector<typename DoFHandlerType::active_cell_iterator> ordered_cells;
2046 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
2047 internal::ClockCells<DoFHandlerType::space_dimension> comparator(center,
2050 typename DoFHandlerType::active_cell_iterator p = dof.begin_active();
2051 typename DoFHandlerType::active_cell_iterator end = dof.end();
2055 ordered_cells.push_back(p);
2058 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
2060 std::vector<types::global_dof_index> reverse(new_indices.size());
2066 template <
typename DoFHandlerType>
2069 const unsigned int level,
2073 std::vector<typename DoFHandlerType::level_cell_iterator> ordered_cells;
2074 ordered_cells.reserve(dof.get_triangulation().n_active_cells());
2075 internal::ClockCells<DoFHandlerType::space_dimension> comparator(center,
2078 typename DoFHandlerType::level_cell_iterator p = dof.begin(level);
2079 typename DoFHandlerType::level_cell_iterator end = dof.end(level);
2083 ordered_cells.push_back(p);
2086 std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
2093 template <
typename DoFHandlerType>
2097 std::vector<types::global_dof_index> renumbering(
2101 dof_handler.renumber_dofs(renumbering);
2106 template <
typename DoFHandlerType>
2109 const DoFHandlerType & dof_handler)
2112 Assert(new_indices.size() == n_dofs,
2115 for (
unsigned int i = 0; i < n_dofs; ++i)
2120 ::boost::mt19937 random_number_generator;
2121 for (
unsigned int i = 1; i < n_dofs; ++i)
2124 const unsigned int j =
2125 ::boost::random::uniform_int_distribution<>(0, i)(
2126 random_number_generator);
2130 std::swap(new_indices[i], new_indices[j]);
2136 template <
typename DoFHandlerType>
2143 DoFHandlerType::space_dimension
> *>(
2144 &dof_handler.get_triangulation())),
2146 "Parallel triangulations are already enumerated according to their MPI process id."));
2148 std::vector<types::global_dof_index> renumbering(
2152 dof_handler.renumber_dofs(renumbering);
2157 template <
typename DoFHandlerType>
2160 const DoFHandlerType & dof_handler)
2163 Assert(new_dof_indices.size() == n_dofs,
2169 std::vector<types::subdomain_id> subdomain_association(n_dofs);
2171 const unsigned int n_subdomains =
2172 *std::max_element(subdomain_association.begin(),
2173 subdomain_association.end()) +
2183 std::fill(new_dof_indices.begin(),
2184 new_dof_indices.end(),
2190 if (subdomain_association[i] == subdomain)
2194 new_dof_indices[i] = next_free_index;
2200 Assert(std::find(new_dof_indices.begin(),
2201 new_dof_indices.end(),
2211 #include "dof_renumbering.inst" 2214 DEAL_II_NAMESPACE_CLOSE
void compute_downstream(std::vector< types::global_dof_index > &new_dof_indices, std::vector< types::global_dof_index > &reverse, const DoFHandlerType &dof_handler, const Tensor< 1, DoFHandlerType::space_dimension > &direction, const bool dof_wise_renumbering)
active_cell_iterator begin_active(const unsigned int level=0) const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
#define AssertDimension(dim1, dim2)
void downstream(DoFHandlerType &dof_handler, const Tensor< 1, DoFHandlerType::space_dimension > &direction, const bool dof_wise_renumbering=false)
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 hierarchical(DoFHandler< dim > &dof_handler)
size_type n_elements() const
void component_wise(DoFHandlerType &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
unsigned int size() const
Transformed quadrature points.
const Triangulation< dim, spacedim > & get_triangulation() const
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)
const std::vector< types::global_dof_index > & n_locally_owned_dofs_per_processor() const
bool is_primitive() const
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
unsigned long long int global_dof_index
static::ExceptionBase & ExcDoFHandlerNotInitialized()
active_cell_iterator begin_active(const unsigned int level=0) const
typename ActiveSelector::active_cell_iterator active_cell_iterator
void random(DoFHandlerType &dof_handler)
static::ExceptionBase & ExcMessage(std::string arg1)
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 subdomain_id
void Cuthill_McKee(DoFHandlerType &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false)
void subtract_set(const IndexSet &other)
#define Assert(cond, exc)
types::global_dof_index n_dofs() const
void compute_subdomain_wise(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
unsigned int n_locally_owned_dofs() const
const FiniteElement< dim, spacedim > & get_fe() const
const ComponentMask & get_nonzero_components(const unsigned int i) 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)
void compute_random(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler)
cell_iterator end() const
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)
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
types::global_dof_index size_type
void renumber_dofs(const std::vector< types::global_dof_index > &new_numbers)
void swap(Vector< Number > &u, Vector< Number > &v)
void reinit(const IndexSet &local_constraints=IndexSet())
cell_iterator begin(const unsigned int level=0) const
cell_iterator end() const
#define AssertThrowMPI(error_code)
void compute_clockwise_dg(std::vector< types::global_dof_index > &new_dof_indices, const DoFHandlerType &dof_handler, const Point< DoFHandlerType::space_dimension > ¢er, 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
void block_wise(DoFHandler< dim, spacedim > &dof_handler)
size_type index_within_set(const size_type global_index) const
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)
bool is_element(const size_type index) const
static::ExceptionBase & ExcNotImplemented()
const ::FEValues< dim, spacedim > & get_present_fe_values() const
void subdomain_wise(DoFHandlerType &dof_handler)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void clockwise_dg(DoFHandlerType &dof_handler, const Point< DoFHandlerType::space_dimension > ¢er, const bool counter=false)
const types::global_dof_index invalid_dof_index
const IndexSet & locally_owned_dofs() const
size_type nth_index_in_set(const unsigned int local_index) const
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
types::global_dof_index compute_block_wise(std::vector< types::global_dof_index > &new_dof_indices, const ITERATOR &start, const ENDITERATOR &end, bool is_level_operation)
void push_back(const Quadrature< dim > &new_quadrature)
static::ExceptionBase & ExcInternalError()
void cell_wise(DoFHandlerType &dof_handler, const std::vector< typename DoFHandlerType::active_cell_iterator > &cell_order)