71 template <
int dim,
typename Number =
double>
83 double downstream_size = 0;
85 for (
unsigned int d = 0; d < dim; ++d)
87 downstream_size += (rhs[d] - lhs[d]) * weight;
90 if (downstream_size < 0)
92 else if (downstream_size > 0)
96 for (
unsigned int d = 0; d < dim; ++d)
100 return lhs[d] < rhs[d];
120 template <
int dim,
int spacedim>
121 std::vector<unsigned char>
125 std::vector<unsigned char> local_component_association(
134 local_component_association[i] =
145 const unsigned int first_comp =
150 local_component_association[i] = first_comp;
157 for (
unsigned int c = first_comp; c < fe.
n_components(); ++c)
158 if (component_mask[c] ==
true)
160 local_component_association[i] = c;
165 Assert(std::find(local_component_association.begin(),
166 local_component_association.end(),
167 static_cast<unsigned char>(-1)) ==
168 local_component_association.end(),
171 return local_component_association;
191 template <
int dim,
int spacedim>
195 std::vector<unsigned char> &dofs_by_component)
197 const ::hp::FECollection<dim, spacedim> &fe_collection =
211 std::vector<std::vector<unsigned char>> local_component_association(
212 fe_collection.size());
213 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
216 local_component_association[f] =
221 std::vector<types::global_dof_index> indices;
226 const unsigned int dofs_per_cell = c->get_fe().n_dofs_per_cell();
227 indices.resize(dofs_per_cell);
228 c->get_dof_indices(indices);
229 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
232 indices[i])] = local_component_association[fe_index][i];
244 template <
int dim,
int spacedim>
247 std::vector<unsigned char> &dofs_by_block)
249 const ::hp::FECollection<dim, spacedim> &fe_collection =
263 std::vector<std::vector<unsigned char>> local_block_association(
264 fe_collection.size());
265 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
269 static_cast<unsigned char>(-1));
273 Assert(std::find(local_block_association[f].begin(),
274 local_block_association[f].end(),
275 static_cast<unsigned char>(-1)) ==
276 local_block_association[f].end(),
281 std::vector<types::global_dof_index> indices;
283 if (cell->is_locally_owned())
286 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
287 indices.resize(dofs_per_cell);
288 cell->get_dof_indices(indices);
289 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
292 indices[i])] = local_block_association[fe_index][i];
299 template <
int dim,
int spacedim,
typename Number>
304 const unsigned int component)
313 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
315 Assert(fe_collection[i].is_primitive() ==
true,
322 const bool consider_components =
326 if (consider_components ==
false)
330 std::vector<unsigned char> component_dofs(
337 for (
unsigned int i = 0; i < dof_data.
size(); ++i)
338 if (component_dofs[i] ==
static_cast<unsigned char>(component))
343 std::vector<unsigned char> touch_count(dof_handler.
n_dofs(), 0);
345 std::vector<types::global_dof_index> dof_indices;
346 dof_indices.reserve(fe_collection.max_dofs_per_cell());
350 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
351 dof_indices.resize(dofs_per_cell);
352 cell->get_dof_indices(dof_indices);
354 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
357 if (!consider_components ||
358 (cell->get_fe().system_to_component_index(i).first == component))
361 dof_data(dof_indices[i]) += cell_data(cell->active_cell_index());
364 ++touch_count[dof_indices[i]];
374 Assert(consider_components || (touch_count[i] != 0),
376 if (touch_count[i] != 0)
377 dof_data(i) /= touch_count[i];
383 template <
int dim,
int spacedim>
391 "The given component mask is not sized correctly to represent the "
392 "components of the given finite element."));
410 std::vector<types::global_dof_index> selected_dofs;
413 if (component_mask[dofs_by_component[i]] ==
true)
418 result.
add_indices(selected_dofs.begin(), selected_dofs.end());
424 template <
int dim,
int spacedim>
430 return extract_dofs<dim, spacedim>(
436 template <
int dim,
int spacedim>
437 std::vector<IndexSet>
444 "The given component mask is not sized correctly to represent the "
445 "components of the given finite element."));
454 std::vector<IndexSet> index_per_comp(n_comps,
IndexSet(dof.
n_dofs()));
458 const auto &comp_i = dofs_by_component[i];
459 if (component_mask[comp_i])
460 index_per_comp[comp_i].add_index(
461 locally_owned_dofs.nth_index_in_set(i));
463 for (
auto &c : index_per_comp)
465 return index_per_comp;
470 template <
int dim,
int spacedim>
475 std::vector<bool> &selected_dofs)
482 "The given component mask is not sized correctly to represent the "
483 "components of the given finite element."));
492 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
false);
499 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
true);
504 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
false);
508 std::vector<unsigned char> local_component_association =
512 local_selected_dofs[i] = component_mask[local_component_association[i]];
518 cell->get_mg_dof_indices(indices);
520 selected_dofs[indices[i]] = local_selected_dofs[i];
526 template <
int dim,
int spacedim>
531 std::vector<bool> &selected_dofs)
542 template <
int dim,
int spacedim>
546 std::vector<bool> &selected_dofs,
547 const std::set<types::boundary_id> &boundary_ids)
553 "This function can not be used with distributed triangulations. "
554 "See the documentation for more information."));
560 selected_dofs.clear();
561 selected_dofs.resize(dof_handler.
n_dofs(),
false);
569 template <
int dim,
int spacedim>
574 const std::set<types::boundary_id> &boundary_ids)
583 template <
int dim,
int spacedim>
587 const std::set<types::boundary_id> &boundary_ids)
591 ExcMessage(
"Component mask has invalid size."));
600 const bool check_boundary_id = (boundary_ids.size() != 0);
604 const bool check_vector_component =
610 std::vector<types::global_dof_index> face_dof_indices;
611 face_dof_indices.reserve(
623 if (cell->is_artificial() ==
false)
624 for (
const unsigned int face : cell->face_indices())
625 if (cell->at_boundary(face))
626 if (!check_boundary_id ||
627 (boundary_ids.find(cell->face(face)->boundary_id()) !=
632 const auto reference_cell = cell->reference_cell();
634 const unsigned int n_vertices_per_cell =
635 reference_cell.n_vertices();
636 const unsigned int n_lines_per_cell = reference_cell.n_lines();
637 const unsigned int n_vertices_per_face =
638 reference_cell.face_reference_cell(face).n_vertices();
639 const unsigned int n_lines_per_face =
640 reference_cell.face_reference_cell(face).n_lines();
643 face_dof_indices.resize(dofs_per_face);
644 cell->face(face)->get_dof_indices(face_dof_indices,
645 cell->active_fe_index());
648 if (!check_vector_component)
649 selected_dofs.
add_index(face_dof_indices[i]);
664 (dim == 3 ? (i < n_vertices_per_face *
667 (i < n_vertices_per_face *
671 (i - n_vertices_per_face *
673 n_vertices_per_cell *
676 n_vertices_per_face *
680 n_vertices_per_cell *
690 selected_dofs.
add_index(face_dof_indices[i]);
694 const unsigned int first_nonzero_comp =
700 if (component_mask[first_nonzero_comp] ==
true)
701 selected_dofs.
add_index(face_dof_indices[i]);
706 return selected_dofs;
711 template <
int dim,
int spacedim>
716 std::vector<bool> &selected_dofs,
717 const std::set<types::boundary_id> &boundary_ids)
721 ExcMessage(
"This component mask has the wrong size."));
728 const bool check_boundary_id = (boundary_ids.size() != 0);
732 const bool check_vector_component =
736 selected_dofs.clear();
737 selected_dofs.resize(dof_handler.
n_dofs(),
false);
738 std::vector<types::global_dof_index> cell_dof_indices;
739 cell_dof_indices.reserve(
749 for (
const unsigned int face : cell->face_indices())
750 if (cell->at_boundary(face))
751 if (!check_boundary_id ||
752 (boundary_ids.find(cell->face(face)->boundary_id()) !=
758 cell_dof_indices.resize(dofs_per_cell);
759 cell->get_dof_indices(cell_dof_indices);
764 if (!check_vector_component)
765 selected_dofs[cell_dof_indices[i]] =
true;
772 selected_dofs[cell_dof_indices[i]] =
777 const unsigned int first_nonzero_comp =
783 selected_dofs[cell_dof_indices[i]] =
784 (component_mask[first_nonzero_comp] ==
true);
793 template <
int dim,
int spacedim,
typename number>
802 const std::function<
bool(
807 ->
bool {
return cell->is_locally_owned() && predicate(cell); };
809 std::vector<types::global_dof_index> local_dof_indices;
810 local_dof_indices.reserve(
814 std::set<types::global_dof_index> predicate_dofs;
817 if (!cell->is_artificial() && predicate(cell))
819 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
820 cell->get_dof_indices(local_dof_indices);
821 predicate_dofs.insert(local_dof_indices.begin(),
822 local_dof_indices.end());
826 std::set<types::global_dof_index> dofs_with_support_on_halo_cells;
828 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
832 active_cell_iterator>::const_iterator it =
834 it != halo_cells.end();
846 const unsigned int dofs_per_cell = (*it)->get_fe().n_dofs_per_cell();
847 local_dof_indices.resize(dofs_per_cell);
848 (*it)->get_dof_indices(local_dof_indices);
849 dofs_with_support_on_halo_cells.insert(local_dof_indices.begin(),
850 local_dof_indices.end());
862 std::set<types::global_dof_index> extra_halo;
863 for (
const auto dof : dofs_with_support_on_halo_cells)
869 const unsigned int line_size = line_ptr->size();
870 for (
unsigned int j = 0; j < line_size; ++j)
871 extra_halo.insert((*line_ptr)[j].first);
874 dofs_with_support_on_halo_cells.insert(extra_halo.begin(),
881 support_set.
add_indices(predicate_dofs.begin(), predicate_dofs.end());
885 halo_set.
add_indices(dofs_with_support_on_halo_cells.begin(),
886 dofs_with_support_on_halo_cells.end());
901 template <
int spacedim>
910 template <
int spacedim>
914 const unsigned int dim = 2;
922 for (
const auto &cell : dof_handler.active_cell_iterators())
923 if (!cell->is_artificial())
925 for (
const unsigned int face : cell->face_indices())
926 if (cell->face(face)->has_children())
929 line = cell->face(face);
933 selected_dofs.add_index(
934 line->child(0)->vertex_dof_index(1, dof));
936 for (
unsigned int child = 0; child < 2; ++child)
938 if (cell->neighbor_child_on_subface(face, child)
943 selected_dofs.add_index(
944 line->child(child)->dof_index(dof));
949 selected_dofs.compress();
950 return selected_dofs;
954 template <
int spacedim>
958 const unsigned int dim = 3;
965 for (
const auto &cell : dof_handler.active_cell_iterators())
966 if (!cell->is_artificial())
967 for (auto f : cell->face_indices())
971 if (cell->face(f)->has_children())
973 for (
unsigned int child = 0; child < 4; ++child)
974 if (!cell->neighbor_child_on_subface(f, child)
978 std::vector<types::global_dof_index> ldi(
980 face->child(child)->get_dof_indices(ldi);
981 selected_dofs.add_indices(ldi.begin(), ldi.end());
986 for (
unsigned int vertex = 0; vertex < 4; ++vertex)
989 unconstrained_dofs.add_index(
990 face->vertex_dof_index(vertex, dof));
993 selected_dofs.subtract_set(unconstrained_dofs);
994 return selected_dofs;
1001 template <
int dim,
int spacedim>
1005 return internal::extract_hanging_node_dofs(dof_handler);
1010 template <
int dim,
int spacedim>
1014 std::vector<bool> &selected_dofs)
1020 std::fill_n(selected_dofs.begin(), dof_handler.
n_dofs(),
false);
1022 std::vector<types::global_dof_index> local_dof_indices;
1023 local_dof_indices.reserve(
1029 if (cell->subdomain_id() == subdomain_id)
1031 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1032 local_dof_indices.resize(dofs_per_cell);
1033 cell->get_dof_indices(local_dof_indices);
1034 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1035 selected_dofs[local_dof_indices[i]] =
true;
1041 template <
int dim,
int spacedim>
1051 std::vector<types::global_dof_index> dof_indices;
1052 std::set<types::global_dof_index> global_dof_indices;
1057 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1058 cell->get_dof_indices(dof_indices);
1062 global_dof_indices.insert(dof_index);
1065 dof_set.
add_indices(global_dof_indices.begin(), global_dof_indices.end());
1074 template <
int dim,
int spacedim>
1084 template <
int dim,
int spacedim>
1088 const unsigned int level)
1096 std::vector<types::global_dof_index> dof_indices;
1097 std::set<types::global_dof_index> global_dof_indices;
1099 const auto filtered_iterators_range =
1102 for (
const auto &cell : filtered_iterators_range)
1104 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1105 cell->get_mg_dof_indices(dof_indices);
1109 global_dof_indices.insert(dof_index);
1112 dof_set.
add_indices(global_dof_indices.begin(), global_dof_indices.end());
1121 template <
int dim,
int spacedim>
1126 const unsigned int level)
1133 template <
int dim,
int spacedim>
1148 std::vector<types::global_dof_index> dof_indices;
1149 std::vector<types::global_dof_index> dofs_on_ghosts;
1152 if (cell->is_ghost())
1154 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1155 cell->get_dof_indices(dof_indices);
1156 for (
const auto dof_index : dof_indices)
1158 dofs_on_ghosts.push_back(dof_index);
1162 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1163 dof_set.
add_indices(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1171 template <
int dim,
int spacedim>
1181 template <
int dim,
int spacedim>
1185 const unsigned int level)
1198 std::vector<types::global_dof_index> dof_indices;
1199 std::vector<types::global_dof_index> dofs_on_ghosts;
1210 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1211 cell->get_mg_dof_indices(dof_indices);
1212 for (
const auto dof_index : dof_indices)
1214 dofs_on_ghosts.push_back(dof_index);
1218 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1219 dof_set.
add_indices(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1227 template <
int dim,
int spacedim>
1231 const unsigned int level,
1239 template <
int dim,
int spacedim>
1243 std::vector<std::vector<bool>> &constant_modes)
1249 constant_modes = std::vector<std::vector<bool>>(0);
1257 std::vector<unsigned char> dofs_by_component(
1262 unsigned int n_selected_dofs = 0;
1263 for (
unsigned int i = 0; i < n_components; ++i)
1264 if (component_mask[i] ==
true)
1266 std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
1270 std::vector<unsigned int> component_numbering(
1272 for (
unsigned int i = 0, count = 0; i < locally_owned_dofs.
n_elements();
1274 if (component_mask[dofs_by_component[i]])
1275 component_numbering[i] = count++;
1282 const ::hp::FECollection<dim, spacedim> &fe_collection =
1284 std::vector<Table<2, bool>> element_constant_modes;
1285 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
1286 constant_mode_to_component_translation(n_components);
1288 unsigned int n_constant_modes = 0;
1289 int first_non_empty_constant_mode = -1;
1290 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
1292 std::pair<Table<2, bool>, std::vector<unsigned int>> data =
1293 fe_collection[f].get_constant_modes();
1297 if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
1299 first_non_empty_constant_mode = f;
1303 for (
unsigned int i = 0; i < data.second.size(); ++i)
1304 if (component_mask[data.second[i]])
1305 constant_mode_to_component_translation[data.second[i]]
1306 .emplace_back(n_constant_modes++, i);
1312 element_constant_modes.push_back(data.first);
1314 element_constant_modes.back().n_rows() == 0 ||
1315 element_constant_modes.back().n_rows() ==
1316 element_constant_modes[first_non_empty_constant_mode].n_rows(),
1323 constant_modes.clear();
1324 constant_modes.resize(n_constant_modes,
1325 std::vector<bool>(n_selected_dofs,
false));
1329 std::vector<types::global_dof_index> dof_indices;
1333 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1334 cell->get_dof_indices(dof_indices);
1336 for (
unsigned int i = 0; i < dof_indices.size(); ++i)
1337 if (locally_owned_dofs.
is_element(dof_indices[i]))
1339 const unsigned int loc_index =
1341 const unsigned int comp = dofs_by_component[loc_index];
1342 if (component_mask[comp])
1343 for (
auto &indices :
1344 constant_mode_to_component_translation[comp])
1345 constant_modes[indices
1346 .
first][component_numbering[loc_index]] =
1347 element_constant_modes[cell->active_fe_index()](
1355 template <
int dim,
int spacedim>
1356 std::map<
typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
1357 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1368 std::map<
typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
1369 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1371 c1_to_c0_cell_and_face;
1374 if (c1_to_c0_face.empty())
1375 return c1_to_c0_cell_and_face;
1378 std::map<typename Triangulation<dim, spacedim>::face_iterator,
1379 typename DoFHandler<dim - 1, spacedim>::active_cell_iterator>
1383 for (
const auto &[c1_cell, c0_cell] : c1_to_c0_face)
1384 if (!c1_cell->has_children())
1385 c0_to_c1[c0_cell] = c1_cell->as_dof_handler_iterator(c1_dh);
1389 for (
const auto &cell :
1391 for (
const auto f : cell->face_indices())
1392 if (cell->face(f)->at_boundary())
1394 const auto &it = c0_to_c1.find(cell->face(f));
1395 if (it != c0_to_c1.end())
1397 const auto &c1_cell = it->second;
1398 c1_to_c0_cell_and_face[c1_cell] = {cell, f};
1404 return c1_to_c0_cell_and_face;
1409 template <
int dim,
int spacedim>
1412 std::vector<unsigned int> &active_fe_indices)
1419 active_fe_indices.assign(indices.begin(), indices.end());
1422 template <
int dim,
int spacedim>
1423 std::vector<IndexSet>
1427 ExcMessage(
"The given DoFHandler has no DoFs."));
1435 "For parallel::distributed::Triangulation objects and "
1436 "associated DoF handler objects, asking for any information "
1437 "related to a subdomain other than the locally owned one does "
1438 "not make sense."));
1442 std::vector<::types::subdomain_id> subdomain_association(
1445 subdomain_association);
1459 const unsigned int n_subdomains =
1463 unsigned int max_subdomain_id = 0;
1466 std::max(max_subdomain_id, cell->subdomain_id());
1467 return max_subdomain_id + 1;
1473 Assert(n_subdomains > *std::max_element(subdomain_association.begin(),
1474 subdomain_association.end()),
1477 std::vector<::IndexSet> index_sets(
1486 index < subdomain_association.size();
1490 if (subdomain_association[index] != this_subdomain)
1492 index_sets[this_subdomain].add_range(i_min, index);
1494 this_subdomain = subdomain_association[index];
1499 if (i_min == subdomain_association.size() - 1)
1501 index_sets[this_subdomain].add_index(i_min);
1507 index_sets[this_subdomain].add_range(i_min,
1508 subdomain_association.size());
1511 for (
unsigned int i = 0; i < n_subdomains; ++i)
1512 index_sets[i].compress();
1517 template <
int dim,
int spacedim>
1518 std::vector<IndexSet>
1528 "For parallel::distributed::Triangulation objects and "
1529 "associated DoF handler objects, asking for any information "
1530 "related to a subdomain other than the locally owned one does "
1531 "not make sense."));
1539 std::vector<IndexSet> dof_set =
1541 const ::types::subdomain_id n_subdomains = dof_set.size();
1547 subdomain_id < n_subdomains;
1560 std::vector<types::global_dof_index> local_dof_indices;
1561 std::set<types::global_dof_index> subdomain_halo_global_dof_indices;
1562 for (
typename std::vector<
1564 const_iterator it_cell = active_halo_layer.begin();
1565 it_cell != active_halo_layer.end();
1571 cell->subdomain_id() != subdomain_id,
1573 "The subdomain ID of the halo cell should not match that of the vector entry."));
1576 cell->get_dof_indices(local_dof_indices);
1580 subdomain_halo_global_dof_indices.insert(local_dof_index);
1583 dof_set[subdomain_id].add_indices(
1584 subdomain_halo_global_dof_indices.begin(),
1585 subdomain_halo_global_dof_indices.end());
1587 dof_set[subdomain_id].compress();
1593 template <
int dim,
int spacedim>
1597 std::vector<types::subdomain_id> &subdomain_association)
1605 "For parallel::distributed::Triangulation objects and "
1606 "associated DoF handler objects, asking for any subdomain other "
1607 "than the locally owned one does not make sense."));
1609 Assert(subdomain_association.size() == dof_handler.
n_dofs(),
1623 std::vector<types::subdomain_id> cell_owners(
1629 cell_owners = tr->get_true_subdomain_ids_of_cells();
1630 Assert(tr->get_true_subdomain_ids_of_cells().size() ==
1631 tr->n_active_cells(),
1638 cell_owners[cell->active_cell_index()] = cell->subdomain_id();
1642 std::fill_n(subdomain_association.begin(),
1646 std::vector<types::global_dof_index> local_dof_indices;
1647 local_dof_indices.reserve(
1655 cell_owners[cell->active_cell_index()];
1656 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1657 local_dof_indices.resize(dofs_per_cell);
1658 cell->get_dof_indices(local_dof_indices);
1664 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1665 if (subdomain_association[local_dof_indices[i]] ==
1667 subdomain_association[local_dof_indices[i]] = subdomain_id;
1668 else if (subdomain_association[local_dof_indices[i]] > subdomain_id)
1670 subdomain_association[local_dof_indices[i]] = subdomain_id;
1674 Assert(std::find(subdomain_association.begin(),
1675 subdomain_association.end(),
1677 subdomain_association.end(),
1683 template <
int dim,
int spacedim>
1689 std::vector<types::subdomain_id> subdomain_association(
1693 return std::count(subdomain_association.begin(),
1694 subdomain_association.end(),
1700 template <
int dim,
int spacedim>
1713 "For parallel::distributed::Triangulation objects and "
1714 "associated DoF handler objects, asking for any subdomain other "
1715 "than the locally owned one does not make sense."));
1719 std::vector<types::global_dof_index> local_dof_indices;
1720 local_dof_indices.reserve(
1728 std::vector<types::global_dof_index> subdomain_indices;
1731 if ((cell->is_artificial() ==
false) &&
1732 (cell->subdomain_id() == subdomain))
1734 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1735 local_dof_indices.resize(dofs_per_cell);
1736 cell->get_dof_indices(local_dof_indices);
1737 subdomain_indices.insert(subdomain_indices.end(),
1738 local_dof_indices.begin(),
1739 local_dof_indices.end());
1742 std::sort(subdomain_indices.begin(), subdomain_indices.end());
1743 index_set.
add_indices(subdomain_indices.begin(), subdomain_indices.end());
1751 template <
int dim,
int spacedim>
1756 std::vector<unsigned int> &n_dofs_on_subdomain)
1761 std::fill(n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
1770 return cell.subdomain_id() == subdomain;
1772 ExcMessage(
"There are no cells for the given subdomain!"));
1774 std::vector<types::subdomain_id> subdomain_association(
1778 std::vector<unsigned char> component_association(dof_handler.
n_dofs());
1781 component_association);
1783 for (
unsigned int c = 0; c < dof_handler.
get_fe(0).n_components(); ++c)
1786 if ((subdomain_association[i] == subdomain) &&
1787 (component_association[i] ==
static_cast<unsigned char>(c)))
1788 ++n_dofs_on_subdomain[c];
1800 template <
int dim,
int spacedim>
1803 const std::vector<unsigned char> &dofs_by_component,
1804 const std::vector<unsigned int> &target_component,
1805 const bool only_once,
1806 std::vector<types::global_dof_index> &dofs_per_component,
1807 unsigned int &component)
1826 for (
unsigned int dd = 0; dd < d; ++dd, ++component)
1827 dofs_per_component[target_component[component]] +=
1828 std::count(dofs_by_component.begin(),
1829 dofs_by_component.end(),
1836 for (
unsigned int dd = 1; dd < d; ++dd)
1837 dofs_per_component[target_component[component - d + dd]] =
1838 dofs_per_component[target_component[component - d]];
1845 template <
int dim,
int spacedim>
1848 const std::vector<unsigned char> &dofs_by_component,
1849 const std::vector<unsigned int> &target_component,
1850 const bool only_once,
1851 std::vector<types::global_dof_index> &dofs_per_component,
1852 unsigned int &component)
1857 for (
unsigned int fe = 1; fe < fe_collection.
size(); ++fe)
1859 Assert(fe_collection[fe].n_components() ==
1860 fe_collection[0].n_components(),
1862 Assert(fe_collection[fe].n_base_elements() ==
1863 fe_collection[0].n_base_elements(),
1865 for (
unsigned int b = 0; b < fe_collection[0].n_base_elements(); ++b)
1867 Assert(fe_collection[fe].base_element(b).n_components() ==
1868 fe_collection[0].base_element(b).n_components(),
1870 Assert(fe_collection[fe].base_element(b).n_base_elements() ==
1871 fe_collection[0].base_element(b).n_base_elements(),
1894 template <
int dim,
int spacedim>
1906 template <
int dim,
int spacedim>
1908 all_elements_are_primitive(
1909 const ::hp::FECollection<dim, spacedim> &fe_collection)
1911 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
1912 if (fe_collection[i].is_primitive() ==
false)
1922 template <
int dim,
int spacedim>
1923 std::vector<types::global_dof_index>
1926 const bool only_once,
1927 const std::vector<unsigned int> &target_component_)
1933 std::vector<unsigned int> target_component = target_component_;
1934 if (target_component.empty())
1936 target_component.resize(n_components);
1937 for (
unsigned int i = 0; i < n_components; ++i)
1938 target_component[i] = i;
1941 Assert(target_component.size() == n_components,
1945 const unsigned int max_component =
1946 *std::max_element(target_component.begin(), target_component.end());
1947 const unsigned int n_target_components = max_component + 1;
1949 std::vector<types::global_dof_index> dofs_per_component(
1954 if (n_components == 1)
1957 return dofs_per_component;
1963 std::vector<unsigned char> dofs_by_component(
1970 unsigned int component = 0;
1981 Assert((internal::all_elements_are_primitive(
1983 (std::accumulate(dofs_per_component.begin(),
1984 dofs_per_component.end(),
1990#ifdef DEAL_II_WITH_MPI
1996 std::vector<types::global_dof_index> local_dof_count =
1999 const int ierr = MPI_Allreduce(local_dof_count.data(),
2000 dofs_per_component.data(),
2001 n_target_components,
2009 return dofs_per_component;
2014 template <
int dim,
int spacedim>
2015 std::vector<types::global_dof_index>
2017 const std::vector<unsigned int> &target_block_)
2019 const ::hp::FECollection<dim, spacedim> &fe_collection =
2028 const unsigned int n_blocks = fe_collection[0].n_blocks();
2030 std::vector<unsigned int> target_block = target_block_;
2031 if (target_block.empty())
2033 target_block.resize(fe_collection[0].n_blocks());
2034 for (
unsigned int i = 0; i < n_blocks; ++i)
2035 target_block[i] = i;
2038 Assert(target_block.size() == n_blocks,
2040 for (
unsigned int f = 1; f < fe_collection.size(); ++f)
2041 Assert(fe_collection[0].n_blocks() == fe_collection[f].n_blocks(),
2042 ExcMessage(
"This function can only work if all elements in a "
2043 "collection have the same number of blocks."));
2049 std::vector<types::global_dof_index> dofs_per_block(1);
2050 dofs_per_block[0] = dof_handler.
n_dofs();
2051 return dofs_per_block;
2055 const unsigned int max_block =
2056 *std::max_element(target_block.begin(), target_block.end());
2057 const unsigned int n_target_blocks = max_block + 1;
2059 std::vector<types::global_dof_index> dofs_per_block(n_target_blocks);
2063 for (
unsigned int this_fe = fe_collection.size() - 1;
2064 this_fe < fe_collection.size();
2069 std::vector<unsigned char> dofs_by_block(
2074 for (
unsigned int block = 0; block < fe.
n_blocks(); ++block)
2075 dofs_per_block[target_block[block]] +=
2076 std::count(dofs_by_block.begin(), dofs_by_block.end(), block);
2078#ifdef DEAL_II_WITH_MPI
2085 std::vector<types::global_dof_index> local_dof_count =
2087 const int ierr = MPI_Allreduce(local_dof_count.data(),
2088 dofs_per_block.data(),
2098 return dofs_per_block;
2103 template <
int dim,
int spacedim>
2106 std::vector<types::global_dof_index> &mapping)
2109 mapping.insert(mapping.end(),
2113 std::vector<types::global_dof_index> dofs_on_face;
2124 for (
const unsigned int f : cell->face_indices())
2125 if (cell->at_boundary(f))
2127 const unsigned int dofs_per_face =
2128 cell->get_fe().n_dofs_per_face(f);
2129 dofs_on_face.resize(dofs_per_face);
2130 cell->face(f)->get_dof_indices(dofs_on_face,
2131 cell->active_fe_index());
2132 for (
unsigned int i = 0; i < dofs_per_face; ++i)
2134 mapping[dofs_on_face[i]] = next_boundary_index++;
2142 template <
int dim,
int spacedim>
2145 const std::set<types::boundary_id> &boundary_ids,
2146 std::vector<types::global_dof_index> &mapping)
2153 mapping.insert(mapping.end(),
2158 if (boundary_ids.empty())
2161 std::vector<types::global_dof_index> dofs_on_face;
2166 for (
const unsigned int f : cell->face_indices())
2167 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
2170 const unsigned int dofs_per_face =
2171 cell->get_fe().n_dofs_per_face(f);
2172 dofs_on_face.resize(dofs_per_face);
2173 cell->face(f)->get_dof_indices(dofs_on_face,
2174 cell->active_fe_index());
2175 for (
unsigned int i = 0; i < dofs_per_face; ++i)
2177 mapping[dofs_on_face[i]] = next_boundary_index++;
2188 template <
int dim,
int spacedim>
2189 std::map<types::global_dof_index, Point<spacedim>>
2195 std::map<types::global_dof_index, Point<spacedim>> support_points;
2201 for (
unsigned int fe_index = 0; fe_index < fe_collection.
size();
2205 Assert((fe_collection[fe_index].n_dofs_per_cell() == 0) ||
2206 (fe_collection[fe_index].has_support_points()),
2209 fe_collection[fe_index].get_unit_support_points()));
2214 (in_mask.
size() == 0 ?
2230 std::vector<types::global_dof_index> local_dof_indices;
2231 for (
const auto &cell : dof_handler.active_cell_iterators())
2233 if (cell->is_artificial() == false)
2235 hp_fe_values.reinit(cell);
2239 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2240 cell->get_dof_indices(local_dof_indices);
2242 const std::vector<Point<spacedim>> &points =
2244 for (
unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell();
2247 const unsigned int dof_comp =
2248 cell->get_fe().system_to_component_index(i).first;
2252 support_points[local_dof_indices[i]] = points[i];
2256 return support_points;
2260 template <
int dim,
int spacedim>
2261 std::vector<Point<spacedim>>
2262 map_dofs_to_support_points_vector(
2267 std::vector<Point<spacedim>> support_points(dof_handler.
n_dofs());
2270 const std::map<types::global_dof_index, Point<spacedim>>
2278 Assert(x_support_points.find(i) != x_support_points.end(),
2281 support_points[i] = x_support_points.find(i)->second;
2284 return support_points;
2290 template <
int dim,
int spacedim>
2302 "This function can not be used with distributed triangulations. "
2303 "See the documentation for more information."));
2310 internal::map_dofs_to_support_points_vector(mapping_collection,
2316 template <
int dim,
int spacedim>
2329 "This function can not be used with distributed triangulations. "
2330 "See the documentation for more information."));
2335 internal::map_dofs_to_support_points_vector(mapping, dof_handler, mask);
2340 template <
int dim,
int spacedim>
2348 support_points.clear();
2354 support_points = internal::map_dofs_to_support_points(mapping_collection,
2361 template <
int dim,
int spacedim>
2369 support_points.clear();
2374 internal::map_dofs_to_support_points(mapping, dof_handler, mask);
2378 template <
int dim,
int spacedim>
2379 std::map<types::global_dof_index, Point<spacedim>>
2388 return internal::map_dofs_to_support_points(mapping_collection,
2394 template <
int dim,
int spacedim>
2395 std::map<types::global_dof_index, Point<spacedim>>
2401 return internal::map_dofs_to_support_points(mapping, dof_handler, mask);
2405 template <
int spacedim>
2413 std::map<Point<spacedim>,
2414 std::vector<types::global_dof_index>,
2419 for (
const auto &it : support_points)
2421 std::vector<types::global_dof_index> &v = point_map[it.second];
2422 v.push_back(it.first);
2426 for (
const auto &it : point_map)
2428 out << it.first <<
" \"";
2429 const std::vector<types::global_dof_index> &v = it.second;
2430 for (
unsigned int i = 0; i < v.size(); ++i)
2444 template <
int dim,
int spacedim>
2453 const unsigned int nb = fe.
n_blocks();
2455 tables_by_block.resize(1);
2456 tables_by_block[0].reinit(nb, nb);
2457 tables_by_block[0].fill(
none);
2465 tables_by_block[0](ib, jb) |= table(i, j);
2473 tables_by_block.resize(fe_collection.
size());
2475 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
2479 const unsigned int nb = fe.
n_blocks();
2480 tables_by_block[f].reinit(nb, nb);
2481 tables_by_block[f].fill(
none);
2488 tables_by_block[f](ib, jb) |= table(i, j);
2497 template <
int dim,
int spacedim>
2501 const unsigned int level,
2502 const std::vector<bool> &selected_dofs,
2505 std::vector<types::global_dof_index> indices;
2510 if (cell->is_locally_owned_on_level())
2517 if (cell->is_locally_owned_on_level())
2519 indices.resize(cell->get_fe().n_dofs_per_cell());
2520 cell->get_mg_dof_indices(indices);
2522 if (selected_dofs.size() != 0)
2527 if (selected_dofs.empty())
2528 block_list.
add(i, indices[j] - offset);
2531 if (selected_dofs[j])
2532 block_list.
add(i, indices[j] - offset);
2540 template <
int dim,
int spacedim>
2544 const unsigned int level,
2545 const bool interior_only)
2550 std::vector<types::global_dof_index> indices;
2551 std::vector<bool> exclude;
2555 indices.resize(cell->get_fe().n_dofs_per_cell());
2556 cell->get_mg_dof_indices(indices);
2562 std::fill(exclude.begin(), exclude.end(),
false);
2564 for (
const unsigned int face : cell->face_indices())
2565 if (cell->at_boundary(face) ||
2566 cell->neighbor(face)->level() != cell->level())
2571 block_list.
add(0, indices[j]);
2575 for (
const auto index : indices)
2576 block_list.
add(0, index);
2582 template <
int dim,
int spacedim>
2586 const unsigned int level,
2587 const bool interior_dofs_only,
2588 const bool boundary_dofs)
2593 std::vector<types::global_dof_index> indices;
2594 std::vector<bool> exclude;
2596 unsigned int block = 0;
2599 if (pcell->is_active())
2602 for (
unsigned int child = 0; child < pcell->n_children(); ++child)
2604 const auto cell = pcell->child(child);
2609 indices.resize(n_dofs);
2610 exclude.resize(n_dofs);
2611 std::fill(exclude.begin(), exclude.end(),
false);
2612 cell->get_mg_dof_indices(indices);
2614 if (interior_dofs_only)
2618 for (
unsigned int d = 0; d < dim; ++d)
2620 const unsigned int face =
2629 for (
const unsigned int face :
2631 if (cell->at_boundary(face))
2637 for (
unsigned int i = 0; i < n_dofs; ++i)
2639 block_list.
add(block, indices[i]);
2645 template <
int dim,
int spacedim>
2646 std::vector<unsigned int>
2649 const unsigned int level,
2650 const bool interior_only,
2651 const bool boundary_patches,
2652 const bool level_boundary_patches,
2653 const bool single_cell_patches,
2654 const bool invert_vertex_mapping)
2661 exclude_boundary_dofs,
2663 level_boundary_patches,
2664 single_cell_patches,
2665 invert_vertex_mapping);
2668 template <
int dim,
int spacedim>
2669 std::vector<unsigned int>
2672 const unsigned int level,
2674 const bool boundary_patches,
2675 const bool level_boundary_patches,
2676 const bool single_cell_patches,
2677 const bool invert_vertex_mapping)
2681 std::vector<unsigned int> vertex_cell_count(
2685 std::vector<bool> vertex_boundary(
2688 std::vector<unsigned int> vertex_mapping(
2693 std::vector<unsigned int> vertex_dof_count(
2699 for (
const unsigned int v : cell->vertex_indices())
2701 const unsigned int vg = cell->vertex_index(v);
2702 vertex_dof_count[vg] += cell->get_fe().n_dofs_per_cell();
2703 ++vertex_cell_count[vg];
2704 for (
unsigned int d = 0; d < dim; ++d)
2707 if (cell->at_boundary(face))
2708 vertex_boundary[vg] =
true;
2709 else if ((!level_boundary_patches) &&
2710 (cell->neighbor(face)->level() !=
2711 static_cast<int>(
level)))
2712 vertex_boundary[vg] =
true;
2718 for (
unsigned int vg = 0; vg < vertex_dof_count.size(); ++vg)
2719 if ((!single_cell_patches && vertex_cell_count[vg] < 2) ||
2720 (!boundary_patches && vertex_boundary[vg]))
2721 vertex_dof_count[vg] = 0;
2724 unsigned int n_vertex_count = 0;
2725 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2726 if (vertex_dof_count[vg] != 0)
2727 vertex_mapping[vg] = n_vertex_count++;
2730 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2731 if (vertex_dof_count[vg] != 0)
2732 vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
2736 vertex_dof_count.resize(n_vertex_count);
2740 block_list.
reinit(vertex_dof_count.size(),
2744 std::vector<types::global_dof_index> indices;
2745 std::vector<bool> exclude;
2751 cell->get_mg_dof_indices(indices);
2753 for (
const unsigned int v : cell->vertex_indices())
2755 const unsigned int vg = cell->vertex_index(v);
2756 const unsigned int block = vertex_mapping[vg];
2762 if (exclude_boundary_dofs.
size() == 0 ||
2768 std::fill(exclude.begin(), exclude.end(),
false);
2770 for (
unsigned int d = 0; d < dim; ++d)
2772 const unsigned int a_face =
2774 const unsigned int face =
2787 for (
unsigned int j = 0; j < indices.size(); ++j)
2789 block_list.
add(block, indices[j]);
2793 for (
const auto index : indices)
2794 block_list.
add(block, index);
2799 if (invert_vertex_mapping)
2802 unsigned int n_vertex_count = 0;
2803 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2805 vertex_mapping[n_vertex_count++] = vg;
2808 vertex_mapping.resize(n_vertex_count);
2811 return vertex_mapping;
2815 template <
int dim,
int spacedim>
2821 std::set<types::global_dof_index> dofs_on_patch;
2822 std::vector<types::global_dof_index> local_dof_indices;
2827 for (
unsigned int i = 0; i < patch.size(); ++i)
2831 Assert(cell->is_artificial() ==
false,
2832 ExcMessage(
"This function can not be called with cells that are "
2833 "not either locally owned or ghost cells."));
2835 cell->get_dof_indices(local_dof_indices);
2836 dofs_on_patch.insert(local_dof_indices.begin(),
2837 local_dof_indices.end());
2841 return dofs_on_patch.size();
2846 template <
int dim,
int spacedim>
2847 std::vector<types::global_dof_index>
2852 std::set<types::global_dof_index> dofs_on_patch;
2853 std::vector<types::global_dof_index> local_dof_indices;
2858 for (
unsigned int i = 0; i < patch.size(); ++i)
2862 Assert(cell->is_artificial() ==
false,
2863 ExcMessage(
"This function can not be called with cells that are "
2864 "not either locally owned or ghost cells."));
2866 cell->get_dof_indices(local_dof_indices);
2867 dofs_on_patch.insert(local_dof_indices.begin(),
2868 local_dof_indices.end());
2871 Assert((dofs_on_patch.size() == count_dofs_on_patch<dim, spacedim>(patch)),
2878 return std::vector<types::global_dof_index>(dofs_on_patch.begin(),
2879 dofs_on_patch.end());
2889#include "dof_tools.inst"
const std::vector< std::pair< size_type, number > > * get_constraint_entries(const size_type line_n) const
size_type n_constraints() const
unsigned int size() const
unsigned int n_selected_blocks(const unsigned int overall_number_of_blocks=numbers::invalid_unsigned_int) const
bool represents_n_components(const unsigned int n) const
bool represents_the_all_selected_mask() const
unsigned int size() const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
const IndexSet & locally_owned_dofs() const
active_cell_iterator begin_active(const unsigned int level=0) const
std::vector< types::fe_index > get_active_fe_indices() const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() const
types::global_dof_index n_locally_owned_dofs() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
const FEValues< dim, spacedim > & get_present_fe_values() const
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_blocks() const
unsigned int n_components() const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
const ComponentMask & get_nonzero_components(const unsigned int i) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const unsigned char combined_orientation=ReferenceCell::default_combined_face_orientation()) const
bool is_primitive() const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
unsigned int component_to_block_index(const unsigned int component) const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
unsigned int element_multiplicity(const unsigned int index) const
unsigned int n_base_elements() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
size_type index_within_set(const size_type global_index) const
size_type n_elements() const
bool is_element(const size_type index) const
void add_index(const size_type index)
void fill_binary_vector(VectorType &vector) const
void subtract_set(const IndexSet &other)
size_type nth_index_in_set(const size_type local_index) const
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Abstract base class for mapping classes.
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
void add(const size_type i, const size_type j)
virtual MPI_Comm get_communicator() const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_cells() const
unsigned int n_levels() const
unsigned int n_vertices() const
virtual size_type size() const override
unsigned int size() const
unsigned int max_dofs_per_face() const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
unsigned int n_components() const
unsigned int max_dofs_per_cell() const
void push_back(const Quadrature< dim_in > &new_quadrature)
virtual MPI_Comm get_communicator() const override
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< FilteredIterator< BaseIterator > > filter_iterators(IteratorRange< BaseIterator > i, const Predicate &p)
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::line_iterator line_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
typename ActiveSelector::CellAccessor cell_accessor
@ update_quadrature_points
Transformed quadrature points.
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
const types::boundary_id internal_face_boundary_id
const types::subdomain_id artificial_subdomain_id
const types::subdomain_id invalid_subdomain_id
static const unsigned int invalid_unsigned_int
const types::global_dof_index invalid_dof_index
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
unsigned short int fe_index
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
#define DEAL_II_DOF_INDEX_MPI_TYPE