75 template <
int dim,
typename Number =
double>
87 double downstream_size = 0;
89 for (
unsigned int d = 0;
d < dim; ++
d)
91 downstream_size += (rhs[
d] - lhs[
d]) * weight;
94 if (downstream_size < 0)
96 else if (downstream_size > 0)
100 for (
unsigned int d = 0;
d < dim; ++
d)
102 if (lhs[
d] == rhs[
d])
104 return lhs[
d] < rhs[
d];
124 template <
int dim,
int spacedim>
125 std::vector<unsigned char>
129 std::vector<unsigned char> local_component_association(
138 local_component_association[i] =
149 const unsigned int first_comp =
154 local_component_association[i] = first_comp;
161 for (
unsigned int c = first_comp; c < fe.
n_components(); ++c)
162 if (component_mask[c] ==
true)
164 local_component_association[i] = c;
169 Assert(std::find(local_component_association.begin(),
170 local_component_association.end(),
171 static_cast<unsigned char>(-1)) ==
172 local_component_association.end(),
175 return local_component_association;
195 template <
int dim,
int spacedim>
199 std::vector<unsigned char> &dofs_by_component)
201 const ::hp::FECollection<dim, spacedim> &fe_collection =
215 std::vector<std::vector<unsigned char>> local_component_association(
216 fe_collection.size());
217 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
220 local_component_association[f] =
225 std::vector<types::global_dof_index> indices;
229 const unsigned int fe_index = c->active_fe_index();
230 const unsigned int dofs_per_cell = c->get_fe().n_dofs_per_cell();
231 indices.resize(dofs_per_cell);
232 c->get_dof_indices(indices);
233 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
236 indices[i])] = local_component_association[fe_index][i];
248 template <
int dim,
int spacedim>
251 std::vector<unsigned char> & dofs_by_block)
253 const ::hp::FECollection<dim, spacedim> &fe_collection =
267 std::vector<std::vector<unsigned char>> local_block_association(
268 fe_collection.size());
269 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
273 static_cast<unsigned char>(-1));
277 Assert(std::find(local_block_association[f].
begin(),
278 local_block_association[f].
end(),
279 static_cast<unsigned char>(-1)) ==
280 local_block_association[f].
end(),
285 std::vector<types::global_dof_index> indices;
287 if (cell->is_locally_owned())
289 const unsigned int fe_index = cell->active_fe_index();
290 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
291 indices.resize(dofs_per_cell);
292 cell->get_dof_indices(indices);
293 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
296 indices[i])] = local_block_association[fe_index][i];
303 template <
int dim,
int spacedim,
typename Number>
308 const unsigned int component)
317 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
319 Assert(fe_collection[i].is_primitive() ==
true,
326 const bool consider_components =
330 if (consider_components ==
false)
334 std::vector<unsigned char> component_dofs(
341 for (
unsigned int i = 0; i < dof_data.
size(); ++i)
342 if (component_dofs[i] ==
static_cast<unsigned char>(component))
347 std::vector<unsigned char> touch_count(dof_handler.
n_dofs(), 0);
349 std::vector<types::global_dof_index> dof_indices;
350 dof_indices.reserve(fe_collection.max_dofs_per_cell());
354 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
355 dof_indices.resize(dofs_per_cell);
356 cell->get_dof_indices(dof_indices);
358 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
361 if (!consider_components ||
362 (cell->get_fe().system_to_component_index(i).first == component))
365 dof_data(dof_indices[i]) += cell_data(cell->active_cell_index());
368 ++touch_count[dof_indices[i]];
378 Assert(consider_components || (touch_count[i] != 0),
380 if (touch_count[i] != 0)
381 dof_data(i) /= touch_count[i];
387 template <
int dim,
int spacedim>
395 "The given component mask is not sized correctly to represent the "
396 "components of the given finite element."));
414 std::vector<types::global_dof_index> selected_dofs;
417 if (component_mask[dofs_by_component[i]] ==
true)
422 result.
add_indices(selected_dofs.begin(), selected_dofs.end());
428 template <
int dim,
int spacedim>
434 return extract_dofs<dim, spacedim>(
440 template <
int dim,
int spacedim>
441 std::vector<IndexSet>
448 "The given component mask is not sized correctly to represent the "
449 "components of the given finite element."));
458 std::vector<IndexSet> index_per_comp(n_comps,
IndexSet(dof.
n_dofs()));
462 const auto &comp_i = dofs_by_component[i];
463 if (component_mask[comp_i])
464 index_per_comp[comp_i].add_index(
465 locally_owned_dofs.nth_index_in_set(i));
467 for (
auto &c : index_per_comp)
469 return index_per_comp;
474 template <
int dim,
int spacedim>
479 std::vector<bool> & selected_dofs)
486 "The given component mask is not sized correctly to represent the "
487 "components of the given finite element."));
496 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
false);
503 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
true);
508 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
false);
512 std::vector<unsigned char> local_component_asssociation =
516 local_selected_dofs[i] = component_mask[local_component_asssociation[i]];
522 cell->get_mg_dof_indices(indices);
524 selected_dofs[indices[i]] = local_selected_dofs[i];
530 template <
int dim,
int spacedim>
535 std::vector<bool> & selected_dofs)
546 template <
int dim,
int spacedim>
550 std::vector<bool> & selected_dofs,
551 const std::set<types::boundary_id> &boundary_ids)
557 "This function can not be used with distributed triangulations. "
558 "See the documentation for more information."));
564 selected_dofs.clear();
565 selected_dofs.resize(dof_handler.
n_dofs(),
false);
573 template <
int dim,
int spacedim>
578 const std::set<types::boundary_id> &boundary_ids)
587 template <
int dim,
int spacedim>
591 const std::set<types::boundary_id> &boundary_ids)
595 ExcMessage(
"Component mask has invalid size."));
604 const bool check_boundary_id = (boundary_ids.size() != 0);
608 const bool check_vector_component =
614 std::vector<types::global_dof_index> face_dof_indices;
615 face_dof_indices.reserve(
627 if (cell->is_artificial() ==
false)
628 for (
const unsigned int face : cell->face_indices())
629 if (cell->at_boundary(face))
630 if (!check_boundary_id ||
631 (boundary_ids.find(cell->face(face)->boundary_id()) !=
638 const unsigned int n_vertices_per_cell =
641 const unsigned int n_vertices_per_face =
643 const unsigned int n_lines_per_face =
647 face_dof_indices.resize(dofs_per_face);
648 cell->face(face)->get_dof_indices(face_dof_indices,
649 cell->active_fe_index());
652 if (!check_vector_component)
653 selected_dofs.
add_index(face_dof_indices[i]);
668 (dim == 3 ? (i < n_vertices_per_face *
671 (i < n_vertices_per_face *
675 (i - n_vertices_per_face *
677 n_vertices_per_cell *
680 n_vertices_per_face *
684 n_vertices_per_cell *
694 selected_dofs.
add_index(face_dof_indices[i]);
698 const unsigned int first_nonzero_comp =
704 if (component_mask[first_nonzero_comp] ==
true)
705 selected_dofs.
add_index(face_dof_indices[i]);
710 return selected_dofs;
715 template <
int dim,
int spacedim>
720 std::vector<bool> & selected_dofs,
721 const std::set<types::boundary_id> &boundary_ids)
725 ExcMessage(
"This component mask has the wrong size."));
732 const bool check_boundary_id = (boundary_ids.size() != 0);
736 const bool check_vector_component =
740 selected_dofs.clear();
741 selected_dofs.resize(dof_handler.
n_dofs(),
false);
742 std::vector<types::global_dof_index> cell_dof_indices;
743 cell_dof_indices.reserve(
753 for (
const unsigned int face : cell->face_indices())
754 if (cell->at_boundary(face))
755 if (!check_boundary_id ||
756 (boundary_ids.find(cell->face(face)->boundary_id()) !=
762 cell_dof_indices.resize(dofs_per_cell);
763 cell->get_dof_indices(cell_dof_indices);
768 if (!check_vector_component)
769 selected_dofs[cell_dof_indices[i]] =
true;
776 selected_dofs[cell_dof_indices[i]] =
781 const unsigned int first_nonzero_comp =
787 selected_dofs[cell_dof_indices[i]] =
788 (component_mask[first_nonzero_comp] ==
true);
797 template <
int dim,
int spacedim,
typename number>
806 const std::function<
bool(
811 ->
bool {
return cell->is_locally_owned() && predicate(cell); };
813 std::vector<types::global_dof_index> local_dof_indices;
814 local_dof_indices.reserve(
818 std::set<types::global_dof_index> predicate_dofs;
821 if (!cell->is_artificial() && predicate(cell))
823 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
824 cell->get_dof_indices(local_dof_indices);
825 predicate_dofs.insert(local_dof_indices.begin(),
826 local_dof_indices.end());
830 std::set<types::global_dof_index> dofs_with_support_on_halo_cells;
832 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
836 active_cell_iterator>::const_iterator it =
838 it != halo_cells.end();
850 const unsigned int dofs_per_cell = (*it)->get_fe().n_dofs_per_cell();
851 local_dof_indices.resize(dofs_per_cell);
852 (*it)->get_dof_indices(local_dof_indices);
853 dofs_with_support_on_halo_cells.insert(local_dof_indices.begin(),
854 local_dof_indices.end());
866 std::set<types::global_dof_index> extra_halo;
867 for (
const auto dof : dofs_with_support_on_halo_cells)
873 const unsigned int line_size = line_ptr->size();
874 for (
unsigned int j = 0; j < line_size; ++j)
875 extra_halo.insert((*line_ptr)[j].first);
878 dofs_with_support_on_halo_cells.insert(extra_halo.begin(),
885 support_set.
add_indices(predicate_dofs.begin(), predicate_dofs.end());
889 halo_set.
add_indices(dofs_with_support_on_halo_cells.begin(),
890 dofs_with_support_on_halo_cells.end());
905 template <
int spacedim>
908 const ::DoFHandler<1, spacedim> &dof_handler)
911 return IndexSet(dof_handler.n_dofs());
915 template <
int spacedim>
918 const ::DoFHandler<2, spacedim> &dof_handler)
920 const unsigned int dim = 2;
922 IndexSet selected_dofs(dof_handler.n_dofs());
928 for (
const auto &cell : dof_handler.active_cell_iterators())
929 if (!cell->is_artificial())
931 for (
const unsigned int face : cell->face_indices())
932 if (cell->face(face)->has_children())
934 const typename ::DoFHandler<dim,
935 spacedim>::line_iterator
936 line = cell->face(face);
940 selected_dofs.add_index(
941 line->child(0)->vertex_dof_index(1, dof));
943 for (
unsigned int child = 0; child < 2; ++child)
945 if (cell->neighbor_child_on_subface(face, child)
950 selected_dofs.add_index(
951 line->child(child)->dof_index(dof));
956 selected_dofs.compress();
957 return selected_dofs;
961 template <
int spacedim>
964 const ::DoFHandler<3, spacedim> &dof_handler)
966 const unsigned int dim = 3;
968 IndexSet selected_dofs(dof_handler.n_dofs());
969 IndexSet unconstrained_dofs(dof_handler.n_dofs());
973 for (
const auto &cell : dof_handler.active_cell_iterators())
974 if (!cell->is_artificial())
975 for (
auto f : cell->face_indices())
977 const typename ::DoFHandler<dim, spacedim>::face_iterator
978 face = cell->face(f);
979 if (cell->face(f)->has_children())
981 for (
unsigned int child = 0; child < 4; ++child)
982 if (!cell->neighbor_child_on_subface(f, child)
986 std::vector<types::global_dof_index> ldi(
988 face->child(child)->get_dof_indices(ldi);
989 selected_dofs.add_indices(ldi.begin(), ldi.end());
994 for (
unsigned int vertex = 0; vertex < 4; ++vertex)
997 unconstrained_dofs.add_index(
998 face->vertex_dof_index(vertex, dof));
1001 selected_dofs.subtract_set(unconstrained_dofs);
1002 return selected_dofs;
1009 template <
int dim,
int spacedim>
1018 template <
int dim,
int spacedim>
1022 std::vector<bool> & selected_dofs)
1028 std::fill_n(selected_dofs.begin(), dof_handler.
n_dofs(),
false);
1030 std::vector<types::global_dof_index> local_dof_indices;
1031 local_dof_indices.reserve(
1039 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1040 local_dof_indices.resize(dofs_per_cell);
1041 cell->get_dof_indices(local_dof_indices);
1042 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1043 selected_dofs[local_dof_indices[i]] =
true;
1049 template <
int dim,
int spacedim>
1059 std::vector<types::global_dof_index> dof_indices;
1060 std::set<types::global_dof_index> global_dof_indices;
1065 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1066 cell->get_dof_indices(dof_indices);
1070 global_dof_indices.insert(dof_index);
1073 dof_set.
add_indices(global_dof_indices.begin(), global_dof_indices.end());
1082 template <
int dim,
int spacedim>
1092 template <
int dim,
int spacedim>
1096 const unsigned int level)
1104 std::vector<types::global_dof_index> dof_indices;
1105 std::set<types::global_dof_index> global_dof_indices;
1107 const auto filtered_iterators_range =
1110 for (
const auto &cell : filtered_iterators_range)
1112 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1113 cell->get_mg_dof_indices(dof_indices);
1117 global_dof_indices.insert(dof_index);
1120 dof_set.
add_indices(global_dof_indices.begin(), global_dof_indices.end());
1129 template <
int dim,
int spacedim>
1134 const unsigned int level)
1141 template <
int dim,
int spacedim>
1156 std::vector<types::global_dof_index> dof_indices;
1157 std::vector<types::global_dof_index> dofs_on_ghosts;
1160 if (cell->is_ghost())
1162 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1163 cell->get_dof_indices(dof_indices);
1164 for (
const auto dof_index : dof_indices)
1166 dofs_on_ghosts.push_back(dof_index);
1170 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1172 std::unique(dofs_on_ghosts.begin(),
1173 dofs_on_ghosts.end()));
1181 template <
int dim,
int spacedim>
1191 template <
int dim,
int spacedim>
1195 const unsigned int level)
1208 std::vector<types::global_dof_index> dof_indices;
1209 std::vector<types::global_dof_index> dofs_on_ghosts;
1220 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1221 cell->get_mg_dof_indices(dof_indices);
1222 for (
const auto dof_index : dof_indices)
1224 dofs_on_ghosts.push_back(dof_index);
1228 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1230 std::unique(dofs_on_ghosts.begin(),
1231 dofs_on_ghosts.end()));
1240 template <
int dim,
int spacedim>
1244 const unsigned int level,
1252 template <
int dim,
int spacedim>
1256 std::vector<std::vector<bool>> & constant_modes)
1262 constant_modes = std::vector<std::vector<bool>>(0);
1270 std::vector<unsigned char> dofs_by_component(
1275 unsigned int n_selected_dofs = 0;
1276 for (
unsigned int i = 0; i < n_components; ++i)
1277 if (component_mask[i] ==
true)
1279 std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
1283 std::vector<unsigned int> component_numbering(
1285 for (
unsigned int i = 0, count = 0; i < locally_owned_dofs.
n_elements();
1287 if (component_mask[dofs_by_component[i]])
1288 component_numbering[i] = count++;
1295 const ::hp::FECollection<dim, spacedim> &fe_collection =
1297 std::vector<Table<2, bool>> element_constant_modes;
1298 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
1299 constant_mode_to_component_translation(n_components);
1301 unsigned int n_constant_modes = 0;
1302 int first_non_empty_constant_mode = -1;
1303 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
1305 std::pair<Table<2, bool>, std::vector<unsigned int>> data =
1306 fe_collection[f].get_constant_modes();
1310 if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
1312 first_non_empty_constant_mode = f;
1316 for (
unsigned int i = 0; i < data.second.size(); ++i)
1317 if (component_mask[data.second[i]])
1318 constant_mode_to_component_translation[data.second[i]]
1319 .emplace_back(n_constant_modes++, i);
1325 element_constant_modes.push_back(data.first);
1327 element_constant_modes.back().n_rows() == 0 ||
1328 element_constant_modes.back().n_rows() ==
1329 element_constant_modes[first_non_empty_constant_mode].n_rows(),
1336 constant_modes.clear();
1337 constant_modes.resize(n_constant_modes,
1338 std::vector<bool>(n_selected_dofs,
false));
1342 std::vector<types::global_dof_index> dof_indices;
1346 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1347 cell->get_dof_indices(dof_indices);
1349 for (
unsigned int i = 0; i < dof_indices.size(); ++i)
1350 if (locally_owned_dofs.
is_element(dof_indices[i]))
1352 const unsigned int loc_index =
1354 const unsigned int comp = dofs_by_component[loc_index];
1355 if (component_mask[comp])
1356 for (
auto &indices :
1357 constant_mode_to_component_translation[comp])
1358 constant_modes[indices
1359 .first][component_numbering[loc_index]] =
1360 element_constant_modes[cell->active_fe_index()](
1368 template <
int dim,
int spacedim>
1371 std::vector<unsigned int> & active_fe_indices)
1377 active_fe_indices[cell->active_cell_index()] = cell->active_fe_index();
1380 template <
int dim,
int spacedim>
1381 std::vector<IndexSet>
1385 ExcMessage(
"The given DoFHandler has no DoFs."));
1393 "For parallel::distributed::Triangulation objects and "
1394 "associated DoF handler objects, asking for any information "
1395 "related to a subdomain other than the locally owned one does "
1396 "not make sense."));
1400 std::vector<::types::subdomain_id> subdomain_association(
1403 subdomain_association);
1417 const unsigned int n_subdomains =
1421 unsigned int max_subdomain_id = 0;
1424 std::max(max_subdomain_id, cell->subdomain_id());
1425 return max_subdomain_id + 1;
1430 ->get_communicator()));
1431 Assert(n_subdomains > *std::max_element(subdomain_association.begin(),
1432 subdomain_association.end()),
1435 std::vector<::IndexSet> index_sets(
1444 index < subdomain_association.size();
1448 if (subdomain_association[index] != this_subdomain)
1450 index_sets[this_subdomain].add_range(i_min, index);
1452 this_subdomain = subdomain_association[index];
1457 if (i_min == subdomain_association.size() - 1)
1459 index_sets[this_subdomain].add_index(i_min);
1465 index_sets[this_subdomain].add_range(i_min,
1466 subdomain_association.size());
1469 for (
unsigned int i = 0; i < n_subdomains; ++i)
1475 template <
int dim,
int spacedim>
1476 std::vector<IndexSet>
1486 "For parallel::distributed::Triangulation objects and "
1487 "associated DoF handler objects, asking for any information "
1488 "related to a subdomain other than the locally owned one does "
1489 "not make sense."));
1497 std::vector<IndexSet> dof_set =
1518 std::vector<types::global_dof_index> local_dof_indices;
1519 std::set<types::global_dof_index> subdomain_halo_global_dof_indices;
1520 for (
typename std::vector<
1522 const_iterator it_cell = active_halo_layer.
begin();
1523 it_cell != active_halo_layer.end();
1531 "The subdomain ID of the halo cell should not match that of the vector entry."));
1534 cell->get_dof_indices(local_dof_indices);
1538 subdomain_halo_global_dof_indices.insert(local_dof_index);
1542 subdomain_halo_global_dof_indices.begin(),
1543 subdomain_halo_global_dof_indices.end());
1551 template <
int dim,
int spacedim>
1555 std::vector<types::subdomain_id> &subdomain_association)
1563 "For parallel::distributed::Triangulation objects and "
1564 "associated DoF handler objects, asking for any subdomain other "
1565 "than the locally owned one does not make sense."));
1567 Assert(subdomain_association.size() == dof_handler.
n_dofs(),
1581 std::vector<types::subdomain_id> cell_owners(
1587 cell_owners = tr->get_true_subdomain_ids_of_cells();
1588 Assert(tr->get_true_subdomain_ids_of_cells().size() ==
1589 tr->n_active_cells(),
1596 cell_owners[cell->active_cell_index()] = cell->subdomain_id();
1600 std::fill_n(subdomain_association.begin(),
1604 std::vector<types::global_dof_index> local_dof_indices;
1605 local_dof_indices.reserve(
1613 cell_owners[cell->active_cell_index()];
1614 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1615 local_dof_indices.resize(dofs_per_cell);
1616 cell->get_dof_indices(local_dof_indices);
1622 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1623 if (subdomain_association[local_dof_indices[i]] ==
1625 subdomain_association[local_dof_indices[i]] =
subdomain_id;
1626 else if (subdomain_association[local_dof_indices[i]] >
subdomain_id)
1628 subdomain_association[local_dof_indices[i]] =
subdomain_id;
1632 Assert(std::find(subdomain_association.begin(),
1633 subdomain_association.end(),
1635 subdomain_association.end(),
1641 template <
int dim,
int spacedim>
1647 std::vector<types::subdomain_id> subdomain_association(
1651 return std::count(subdomain_association.begin(),
1652 subdomain_association.end(),
1658 template <
int dim,
int spacedim>
1671 "For parallel::distributed::Triangulation objects and "
1672 "associated DoF handler objects, asking for any subdomain other "
1673 "than the locally owned one does not make sense."));
1677 std::vector<types::global_dof_index> local_dof_indices;
1678 local_dof_indices.reserve(
1686 std::vector<types::global_dof_index> subdomain_indices;
1689 if ((cell->is_artificial() ==
false) &&
1690 (cell->subdomain_id() == subdomain))
1692 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1693 local_dof_indices.resize(dofs_per_cell);
1694 cell->get_dof_indices(local_dof_indices);
1695 subdomain_indices.insert(subdomain_indices.end(),
1696 local_dof_indices.begin(),
1697 local_dof_indices.end());
1700 std::sort(subdomain_indices.begin(), subdomain_indices.end());
1701 subdomain_indices.erase(std::unique(subdomain_indices.begin(),
1702 subdomain_indices.end()),
1703 subdomain_indices.end());
1706 index_set.
add_indices(subdomain_indices.begin(), subdomain_indices.end());
1714 template <
int dim,
int spacedim>
1719 std::vector<unsigned int> & n_dofs_on_subdomain)
1724 std::fill(n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
1733 return cell.subdomain_id() == subdomain;
1735 ExcMessage(
"There are no cells for the given subdomain!"));
1737 std::vector<types::subdomain_id> subdomain_association(
1741 std::vector<unsigned char> component_association(dof_handler.
n_dofs());
1743 std::vector<bool>(),
1744 component_association);
1746 for (
unsigned int c = 0; c < dof_handler.
get_fe(0).n_components(); ++c)
1749 if ((subdomain_association[i] == subdomain) &&
1750 (component_association[i] ==
static_cast<unsigned char>(c)))
1751 ++n_dofs_on_subdomain[c];
1763 template <
int dim,
int spacedim>
1766 const std::vector<unsigned char> & dofs_by_component,
1767 const std::vector<unsigned int> & target_component,
1768 const bool only_once,
1769 std::vector<types::global_dof_index> &dofs_per_component,
1770 unsigned int & component)
1789 for (
unsigned int dd = 0; dd <
d; ++dd, ++component)
1790 dofs_per_component[target_component[component]] +=
1791 std::count(dofs_by_component.begin(),
1792 dofs_by_component.end(),
1799 for (
unsigned int dd = 1; dd <
d; ++dd)
1800 dofs_per_component[target_component[component -
d + dd]] =
1801 dofs_per_component[target_component[component -
d]];
1808 template <
int dim,
int spacedim>
1811 const std::vector<unsigned char> & dofs_by_component,
1812 const std::vector<unsigned int> & target_component,
1813 const bool only_once,
1814 std::vector<types::global_dof_index> &dofs_per_component,
1815 unsigned int & component)
1820 for (
unsigned int fe = 1; fe < fe_collection.
size(); ++fe)
1822 Assert(fe_collection[fe].n_components() ==
1823 fe_collection[0].n_components(),
1825 Assert(fe_collection[fe].n_base_elements() ==
1826 fe_collection[0].n_base_elements(),
1828 for (
unsigned int b = 0;
b < fe_collection[0].n_base_elements(); ++
b)
1830 Assert(fe_collection[fe].base_element(
b).n_components() ==
1831 fe_collection[0].base_element(
b).n_components(),
1833 Assert(fe_collection[fe].base_element(
b).n_base_elements() ==
1834 fe_collection[0].base_element(
b).n_base_elements(),
1857 template <
int dim,
int spacedim>
1869 template <
int dim,
int spacedim>
1871 all_elements_are_primitive(
1872 const ::hp::FECollection<dim, spacedim> &fe_collection)
1874 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
1875 if (fe_collection[i].is_primitive() ==
false)
1885 template <
int dim,
int spacedim>
1886 std::vector<types::global_dof_index>
1889 const bool only_once,
1890 const std::vector<unsigned int> &target_component_)
1896 std::vector<unsigned int> target_component = target_component_;
1897 if (target_component.size() == 0)
1899 target_component.resize(n_components);
1900 for (
unsigned int i = 0; i < n_components; ++i)
1901 target_component[i] = i;
1904 Assert(target_component.size() == n_components,
1908 const unsigned int max_component =
1909 *std::max_element(target_component.begin(), target_component.end());
1910 const unsigned int n_target_components = max_component + 1;
1912 std::vector<types::global_dof_index> dofs_per_component(
1917 if (n_components == 1)
1920 return dofs_per_component;
1926 std::vector<unsigned char> dofs_by_component(
1933 unsigned int component = 0;
1944 Assert((internal::all_elements_are_primitive(
1946 (std::accumulate(dofs_per_component.begin(),
1947 dofs_per_component.end(),
1953 #ifdef DEAL_II_WITH_MPI
1959 std::vector<types::global_dof_index> local_dof_count =
1962 const int ierr = MPI_Allreduce(local_dof_count.data(),
1963 dofs_per_component.data(),
1964 n_target_components,
1972 return dofs_per_component;
1977 template <
int dim,
int spacedim>
1978 std::vector<types::global_dof_index>
1980 const std::vector<unsigned int> &target_block_)
1982 const ::hp::FECollection<dim, spacedim> &fe_collection =
1991 const unsigned int n_blocks = fe_collection[0].n_blocks();
1993 std::vector<unsigned int> target_block = target_block_;
1994 if (target_block.size() == 0)
1996 target_block.resize(fe_collection[0].
n_blocks());
1997 for (
unsigned int i = 0; i <
n_blocks; ++i)
1998 target_block[i] = i;
2003 for (
unsigned int f = 1; f < fe_collection.size(); ++f)
2005 ExcMessage(
"This function can only work if all elements in a "
2006 "collection have the same number of blocks."));
2012 std::vector<types::global_dof_index> dofs_per_block(1);
2013 dofs_per_block[0] = dof_handler.
n_dofs();
2014 return dofs_per_block;
2018 const unsigned int max_block =
2019 *std::max_element(target_block.begin(), target_block.end());
2020 const unsigned int n_target_blocks = max_block + 1;
2022 std::vector<types::global_dof_index> dofs_per_block(n_target_blocks);
2026 for (
unsigned int this_fe = fe_collection.size() - 1;
2027 this_fe < fe_collection.size();
2032 std::vector<unsigned char> dofs_by_block(
2037 for (
unsigned int block = 0; block < fe.
n_blocks(); ++block)
2038 dofs_per_block[target_block[block]] +=
2039 std::count(dofs_by_block.begin(), dofs_by_block.end(), block);
2041 #ifdef DEAL_II_WITH_MPI
2048 std::vector<types::global_dof_index> local_dof_count =
2050 const int ierr = MPI_Allreduce(local_dof_count.data(),
2051 dofs_per_block.data(),
2061 return dofs_per_block;
2066 template <
int dim,
int spacedim>
2069 std::vector<types::global_dof_index> &mapping)
2072 mapping.insert(mapping.end(),
2076 std::vector<types::global_dof_index> dofs_on_face;
2087 for (
const unsigned int f : cell->face_indices())
2088 if (cell->at_boundary(f))
2090 const unsigned int dofs_per_face =
2091 cell->get_fe().n_dofs_per_face(f);
2092 dofs_on_face.resize(dofs_per_face);
2093 cell->face(f)->get_dof_indices(dofs_on_face,
2094 cell->active_fe_index());
2095 for (
unsigned int i = 0; i < dofs_per_face; ++i)
2097 mapping[dofs_on_face[i]] = next_boundary_index++;
2105 template <
int dim,
int spacedim>
2108 const std::set<types::boundary_id> &boundary_ids,
2109 std::vector<types::global_dof_index> &mapping)
2116 mapping.insert(mapping.end(),
2121 if (boundary_ids.size() == 0)
2124 std::vector<types::global_dof_index> dofs_on_face;
2129 for (
const unsigned int f : cell->face_indices())
2130 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
2133 const unsigned int dofs_per_face =
2134 cell->get_fe().n_dofs_per_face(f);
2135 dofs_on_face.resize(dofs_per_face);
2136 cell->face(f)->get_dof_indices(dofs_on_face,
2137 cell->active_fe_index());
2138 for (
unsigned int i = 0; i < dofs_per_face; ++i)
2140 mapping[dofs_on_face[i]] = next_boundary_index++;
2151 template <
int dim,
int spacedim>
2163 for (
unsigned int fe_index = 0; fe_index < fe_collection.
size();
2167 Assert((fe_collection[fe_index].n_dofs_per_cell() == 0) ||
2168 (fe_collection[fe_index].has_support_points()),
2171 fe_collection[fe_index].get_unit_support_points()));
2176 (in_mask.
size() == 0 ?
2192 std::vector<types::global_dof_index> local_dof_indices;
2195 if (cell->is_artificial() ==
false)
2197 hp_fe_values.reinit(cell);
2201 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2202 cell->get_dof_indices(local_dof_indices);
2204 const std::vector<Point<spacedim>> &points =
2206 for (
unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell();
2209 const unsigned int dof_comp =
2210 cell->get_fe().system_to_component_index(i).first;
2214 support_points[local_dof_indices[i]] = points[i];
2220 template <
int dim,
int spacedim>
2229 std::map<types::global_dof_index, Point<spacedim>> x_support_points;
2239 Assert(x_support_points.find(i) != x_support_points.end(),
2242 support_points[i] = x_support_points[i];
2248 template <
int dim,
int spacedim>
2260 "This function can not be used with distributed triangulations. "
2261 "See the documentation for more information."));
2274 template <
int dim,
int spacedim>
2287 "This function can not be used with distributed triangulations. "
2288 "See the documentation for more information."));
2299 template <
int dim,
int spacedim>
2307 support_points.clear();
2320 template <
int dim,
int spacedim>
2328 support_points.clear();
2338 template <
int spacedim>
2346 std::map<Point<spacedim>,
2347 std::vector<types::global_dof_index>,
2352 for (
const auto &it : support_points)
2354 std::vector<types::global_dof_index> &v = point_map[it.second];
2355 v.push_back(it.first);
2359 for (
const auto &it : point_map)
2361 out << it.first <<
" \"";
2362 const std::vector<types::global_dof_index> &v = it.second;
2363 for (
unsigned int i = 0; i < v.size(); ++i)
2377 template <
int dim,
int spacedim>
2386 const unsigned int nb = fe.
n_blocks();
2388 tables_by_block.resize(1);
2389 tables_by_block[0].reinit(nb, nb);
2390 tables_by_block[0].fill(
none);
2398 tables_by_block[0](ib, jb) |= table(i, j);
2406 tables_by_block.resize(fe_collection.
size());
2408 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
2412 const unsigned int nb = fe.
n_blocks();
2413 tables_by_block[f].reinit(nb, nb);
2414 tables_by_block[f].fill(
none);
2421 tables_by_block[f](ib, jb) |= table(i, j);
2430 template <
int dim,
int spacedim>
2434 const unsigned int level,
2435 const std::vector<bool> & selected_dofs,
2438 std::vector<types::global_dof_index> indices;
2443 if (cell->is_locally_owned_on_level())
2450 if (cell->is_locally_owned_on_level())
2452 indices.resize(cell->get_fe().n_dofs_per_cell());
2453 cell->get_mg_dof_indices(indices);
2455 if (selected_dofs.size() != 0)
2460 if (selected_dofs.size() == 0)
2461 block_list.
add(i, indices[j] - offset);
2464 if (selected_dofs[j])
2465 block_list.
add(i, indices[j] - offset);
2473 template <
int dim,
int spacedim>
2477 const unsigned int level,
2478 const bool interior_only)
2483 std::vector<types::global_dof_index> indices;
2484 std::vector<bool> exclude;
2488 indices.resize(cell->get_fe().n_dofs_per_cell());
2489 cell->get_mg_dof_indices(indices);
2495 std::fill(exclude.begin(), exclude.end(),
false);
2497 for (
const unsigned int face : cell->face_indices())
2498 if (cell->at_boundary(face) ||
2499 cell->neighbor(face)->level() != cell->level())
2504 block_list.
add(0, indices[j]);
2508 for (
const auto index : indices)
2509 block_list.
add(0, index);
2515 template <
int dim,
int spacedim>
2519 const unsigned int level,
2520 const bool interior_dofs_only,
2521 const bool boundary_dofs)
2526 std::vector<types::global_dof_index> indices;
2527 std::vector<bool> exclude;
2529 unsigned int block = 0;
2532 if (pcell->is_active())
2535 for (
unsigned int child = 0; child < pcell->n_children(); ++child)
2537 const auto cell = pcell->child(child);
2542 indices.resize(n_dofs);
2543 exclude.resize(n_dofs);
2544 std::fill(exclude.begin(), exclude.end(),
false);
2545 cell->get_mg_dof_indices(indices);
2547 if (interior_dofs_only)
2551 for (
unsigned int d = 0;
d < dim; ++
d)
2553 const unsigned int face =
2562 for (
const unsigned int face :
2564 if (cell->at_boundary(face))
2570 for (
unsigned int i = 0; i < n_dofs; ++i)
2572 block_list.
add(block, indices[i]);
2578 template <
int dim,
int spacedim>
2579 std::vector<unsigned int>
2582 const unsigned int level,
2583 const bool interior_only,
2584 const bool boundary_patches,
2585 const bool level_boundary_patches,
2586 const bool single_cell_patches,
2587 const bool invert_vertex_mapping)
2594 exclude_boundary_dofs,
2596 level_boundary_patches,
2597 single_cell_patches,
2598 invert_vertex_mapping);
2601 template <
int dim,
int spacedim>
2602 std::vector<unsigned int>
2605 const unsigned int level,
2606 const BlockMask & exclude_boundary_dofs,
2607 const bool boundary_patches,
2608 const bool level_boundary_patches,
2609 const bool single_cell_patches,
2610 const bool invert_vertex_mapping)
2614 std::vector<unsigned int> vertex_cell_count(
2618 std::vector<bool> vertex_boundary(
2621 std::vector<unsigned int> vertex_mapping(
2626 std::vector<unsigned int> vertex_dof_count(
2632 for (
const unsigned int v : cell->vertex_indices())
2634 const unsigned int vg = cell->vertex_index(v);
2635 vertex_dof_count[vg] += cell->get_fe().n_dofs_per_cell();
2636 ++vertex_cell_count[vg];
2637 for (
unsigned int d = 0;
d < dim; ++
d)
2640 if (cell->at_boundary(face))
2641 vertex_boundary[vg] =
true;
2642 else if ((!level_boundary_patches) &&
2643 (cell->neighbor(face)->level() !=
2644 static_cast<int>(
level)))
2645 vertex_boundary[vg] =
true;
2651 for (
unsigned int vg = 0; vg < vertex_dof_count.size(); ++vg)
2652 if ((!single_cell_patches && vertex_cell_count[vg] < 2) ||
2653 (!boundary_patches && vertex_boundary[vg]))
2654 vertex_dof_count[vg] = 0;
2657 unsigned int n_vertex_count = 0;
2658 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2659 if (vertex_dof_count[vg] != 0)
2660 vertex_mapping[vg] = n_vertex_count++;
2663 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2664 if (vertex_dof_count[vg] != 0)
2665 vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
2669 vertex_dof_count.resize(n_vertex_count);
2673 block_list.
reinit(vertex_dof_count.size(),
2677 std::vector<types::global_dof_index> indices;
2678 std::vector<bool> exclude;
2684 cell->get_mg_dof_indices(indices);
2686 for (
const unsigned int v : cell->vertex_indices())
2688 const unsigned int vg = cell->vertex_index(v);
2689 const unsigned int block = vertex_mapping[vg];
2695 if (exclude_boundary_dofs.
size() == 0 ||
2701 std::fill(exclude.begin(), exclude.end(),
false);
2703 for (
unsigned int d = 0;
d < dim; ++
d)
2705 const unsigned int a_face =
2707 const unsigned int face =
2720 for (
unsigned int j = 0; j < indices.size(); ++j)
2722 block_list.
add(block, indices[j]);
2726 for (
const auto index : indices)
2727 block_list.
add(block, index);
2732 if (invert_vertex_mapping)
2735 unsigned int n_vertex_count = 0;
2736 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2738 vertex_mapping[n_vertex_count++] = vg;
2741 vertex_mapping.resize(n_vertex_count);
2744 return vertex_mapping;
2748 template <
int dim,
int spacedim>
2754 std::set<types::global_dof_index> dofs_on_patch;
2755 std::vector<types::global_dof_index> local_dof_indices;
2760 for (
unsigned int i = 0; i < patch.size(); ++i)
2764 Assert(cell->is_artificial() ==
false,
2765 ExcMessage(
"This function can not be called with cells that are "
2766 "not either locally owned or ghost cells."));
2768 cell->get_dof_indices(local_dof_indices);
2769 dofs_on_patch.insert(local_dof_indices.begin(),
2770 local_dof_indices.end());
2774 return dofs_on_patch.size();
2779 template <
int dim,
int spacedim>
2780 std::vector<types::global_dof_index>
2785 std::set<types::global_dof_index> dofs_on_patch;
2786 std::vector<types::global_dof_index> local_dof_indices;
2791 for (
unsigned int i = 0; i < patch.size(); ++i)
2795 Assert(cell->is_artificial() ==
false,
2796 ExcMessage(
"This function can not be called with cells that are "
2797 "not either locally owned or ghost cells."));
2799 cell->get_dof_indices(local_dof_indices);
2800 dofs_on_patch.insert(local_dof_indices.begin(),
2801 local_dof_indices.end());
2804 Assert((dofs_on_patch.size() == count_dofs_on_patch<dim, spacedim>(patch)),
2811 return std::vector<types::global_dof_index>(dofs_on_patch.begin(),
2812 dofs_on_patch.end());
2822 #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
types::global_dof_index n_boundary_dofs() const
const Triangulation< dim, spacedim > & get_triangulation() const
active_cell_iterator begin_active(const unsigned int level=0) const
bool has_hp_capabilities() const
types::global_dof_index n_dofs() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
cell_iterator begin(const unsigned int level=0) const
const IndexSet & locally_owned_dofs() const
const hp::FECollection< dim, spacedim > & get_fe_collection() 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, unsigned int > system_to_component_index(const unsigned int index) const
ComponentMask component_mask(const FEValuesExtractors::Scalar &scalar) const
std::pair< unsigned int, types::global_dof_index > system_to_block_index(const unsigned int component) const
const ComponentMask & get_nonzero_components(const unsigned int i) const
bool is_primitive() const
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
unsigned int component_to_block_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index, const unsigned int face_no=0) const
unsigned int n_base_elements() 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)
void add(const size_type i, const size_type j)
virtual void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths) override
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
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)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
@ update_quadrature_points
Transformed quadrature points.
IteratorRange< FilteredIterator< BaseIterator > > filter_iterators(IteratorRange< BaseIterator > i, const Predicate &p)
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcInvalidBoundaryIndicator()
#define AssertThrow(cond, exc)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if< IsBlockVector< VectorType >::value, unsigned int >::type n_blocks(const VectorType &vector)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
std::string compress(const std::string &input)
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
unsigned int global_dof_index
unsigned int subdomain_id
const ::Triangulation< dim, spacedim > & tria
#define DEAL_II_DOF_INDEX_MPI_TYPE