72 template <
int dim,
typename Number =
double>
84 double downstream_size = 0;
86 for (
unsigned int d = 0;
d < dim; ++
d)
88 downstream_size += (rhs[
d] - lhs[
d]) * weight;
91 if (downstream_size < 0)
93 else if (downstream_size > 0)
97 for (
unsigned int d = 0;
d < dim; ++
d)
101 return lhs[
d] < rhs[
d];
121 template <
int dim,
int spacedim>
122 std::vector<unsigned char>
126 std::vector<unsigned char> local_component_association(
135 local_component_association[i] =
146 const unsigned int first_comp =
151 local_component_association[i] = first_comp;
158 for (
unsigned int c = first_comp; c < fe.
n_components(); ++c)
159 if (component_mask[c] ==
true)
161 local_component_association[i] = c;
166 Assert(std::find(local_component_association.begin(),
167 local_component_association.end(),
168 static_cast<unsigned char>(-1)) ==
169 local_component_association.end(),
172 return local_component_association;
192 template <
int dim,
int spacedim>
196 std::vector<unsigned char> &dofs_by_component)
198 const ::hp::FECollection<dim, spacedim> &fe_collection =
212 std::vector<std::vector<unsigned char>> local_component_association(
213 fe_collection.size());
214 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
217 local_component_association[f] =
222 std::vector<types::global_dof_index> indices;
227 const unsigned int dofs_per_cell = c->get_fe().n_dofs_per_cell();
228 indices.resize(dofs_per_cell);
229 c->get_dof_indices(indices);
230 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
233 indices[i])] = local_component_association[
fe_index][i];
245 template <
int dim,
int spacedim>
248 std::vector<unsigned char> &dofs_by_block)
250 const ::hp::FECollection<dim, spacedim> &fe_collection =
264 std::vector<std::vector<unsigned char>> local_block_association(
265 fe_collection.size());
266 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
270 static_cast<unsigned char>(-1));
274 Assert(std::find(local_block_association[f].
begin(),
275 local_block_association[f].
end(),
276 static_cast<unsigned char>(-1)) ==
277 local_block_association[f].
end(),
282 std::vector<types::global_dof_index> indices;
284 if (cell->is_locally_owned())
287 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
288 indices.resize(dofs_per_cell);
289 cell->get_dof_indices(indices);
290 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
293 indices[i])] = local_block_association[
fe_index][i];
300 template <
int dim,
int spacedim,
typename Number>
305 const unsigned int component)
314 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
316 Assert(fe_collection[i].is_primitive() ==
true,
323 const bool consider_components =
327 if (consider_components ==
false)
331 std::vector<unsigned char> component_dofs(
338 for (
unsigned int i = 0; i < dof_data.
size(); ++i)
339 if (component_dofs[i] ==
static_cast<unsigned char>(component))
344 std::vector<unsigned char> touch_count(dof_handler.
n_dofs(), 0);
346 std::vector<types::global_dof_index> dof_indices;
347 dof_indices.reserve(fe_collection.max_dofs_per_cell());
351 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
352 dof_indices.resize(dofs_per_cell);
353 cell->get_dof_indices(dof_indices);
355 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
358 if (!consider_components ||
359 (cell->get_fe().system_to_component_index(i).first == component))
362 dof_data(dof_indices[i]) += cell_data(cell->active_cell_index());
365 ++touch_count[dof_indices[i]];
375 Assert(consider_components || (touch_count[i] != 0),
377 if (touch_count[i] != 0)
378 dof_data(i) /= touch_count[i];
384 template <
int dim,
int spacedim>
392 "The given component mask is not sized correctly to represent the "
393 "components of the given finite element."));
411 std::vector<types::global_dof_index> selected_dofs;
414 if (component_mask[dofs_by_component[i]] ==
true)
419 result.
add_indices(selected_dofs.begin(), selected_dofs.end());
425 template <
int dim,
int spacedim>
431 return extract_dofs<dim, spacedim>(
437 template <
int dim,
int spacedim>
438 std::vector<IndexSet>
445 "The given component mask is not sized correctly to represent the "
446 "components of the given finite element."));
455 std::vector<IndexSet> index_per_comp(n_comps,
IndexSet(dof.
n_dofs()));
459 const auto &comp_i = dofs_by_component[i];
460 if (component_mask[comp_i])
461 index_per_comp[comp_i].add_index(
462 locally_owned_dofs.nth_index_in_set(i));
464 for (
auto &c : index_per_comp)
466 return index_per_comp;
471 template <
int dim,
int spacedim>
476 std::vector<bool> &selected_dofs)
483 "The given component mask is not sized correctly to represent the "
484 "components of the given finite element."));
493 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
false);
500 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
true);
505 std::fill_n(selected_dofs.begin(), dof.
n_dofs(
level),
false);
509 std::vector<unsigned char> local_component_association =
513 local_selected_dofs[i] = component_mask[local_component_association[i]];
519 cell->get_mg_dof_indices(indices);
521 selected_dofs[indices[i]] = local_selected_dofs[i];
527 template <
int dim,
int spacedim>
532 std::vector<bool> &selected_dofs)
543 template <
int dim,
int spacedim>
547 std::vector<bool> &selected_dofs,
548 const std::set<types::boundary_id> &boundary_ids)
554 "This function can not be used with distributed triangulations. "
555 "See the documentation for more information."));
561 selected_dofs.clear();
562 selected_dofs.resize(dof_handler.
n_dofs(),
false);
570 template <
int dim,
int spacedim>
575 const std::set<types::boundary_id> &boundary_ids)
584 template <
int dim,
int spacedim>
588 const std::set<types::boundary_id> &boundary_ids)
592 ExcMessage(
"Component mask has invalid size."));
601 const bool check_boundary_id = (boundary_ids.size() != 0);
605 const bool check_vector_component =
611 std::vector<types::global_dof_index> face_dof_indices;
612 face_dof_indices.reserve(
624 if (cell->is_artificial() ==
false)
625 for (
const unsigned int face : cell->face_indices())
626 if (cell->at_boundary(face))
627 if (!check_boundary_id ||
628 (boundary_ids.find(cell->face(face)->boundary_id()) !=
635 const unsigned int n_vertices_per_cell =
638 const unsigned int n_vertices_per_face =
640 const unsigned int n_lines_per_face =
644 face_dof_indices.resize(dofs_per_face);
645 cell->face(face)->get_dof_indices(face_dof_indices,
646 cell->active_fe_index());
649 if (!check_vector_component)
650 selected_dofs.
add_index(face_dof_indices[i]);
665 (dim == 3 ? (i < n_vertices_per_face *
668 (i < n_vertices_per_face *
672 (i - n_vertices_per_face *
674 n_vertices_per_cell *
677 n_vertices_per_face *
681 n_vertices_per_cell *
691 selected_dofs.
add_index(face_dof_indices[i]);
695 const unsigned int first_nonzero_comp =
701 if (component_mask[first_nonzero_comp] ==
true)
702 selected_dofs.
add_index(face_dof_indices[i]);
707 return selected_dofs;
712 template <
int dim,
int spacedim>
717 std::vector<bool> &selected_dofs,
718 const std::set<types::boundary_id> &boundary_ids)
722 ExcMessage(
"This component mask has the wrong size."));
729 const bool check_boundary_id = (boundary_ids.size() != 0);
733 const bool check_vector_component =
737 selected_dofs.clear();
738 selected_dofs.resize(dof_handler.
n_dofs(),
false);
739 std::vector<types::global_dof_index> cell_dof_indices;
740 cell_dof_indices.reserve(
750 for (
const unsigned int face : cell->face_indices())
751 if (cell->at_boundary(face))
752 if (!check_boundary_id ||
753 (boundary_ids.find(cell->face(face)->boundary_id()) !=
759 cell_dof_indices.resize(dofs_per_cell);
760 cell->get_dof_indices(cell_dof_indices);
765 if (!check_vector_component)
766 selected_dofs[cell_dof_indices[i]] =
true;
773 selected_dofs[cell_dof_indices[i]] =
778 const unsigned int first_nonzero_comp =
784 selected_dofs[cell_dof_indices[i]] =
785 (component_mask[first_nonzero_comp] ==
true);
794 template <
int dim,
int spacedim,
typename number>
803 const std::function<
bool(
808 ->
bool {
return cell->is_locally_owned() && predicate(cell); };
810 std::vector<types::global_dof_index> local_dof_indices;
811 local_dof_indices.reserve(
815 std::set<types::global_dof_index> predicate_dofs;
818 if (!cell->is_artificial() && predicate(cell))
820 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
821 cell->get_dof_indices(local_dof_indices);
822 predicate_dofs.insert(local_dof_indices.begin(),
823 local_dof_indices.end());
827 std::set<types::global_dof_index> dofs_with_support_on_halo_cells;
829 const std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>
833 active_cell_iterator>::const_iterator it =
835 it != halo_cells.end();
847 const unsigned int dofs_per_cell = (*it)->get_fe().n_dofs_per_cell();
848 local_dof_indices.resize(dofs_per_cell);
849 (*it)->get_dof_indices(local_dof_indices);
850 dofs_with_support_on_halo_cells.insert(local_dof_indices.begin(),
851 local_dof_indices.end());
863 std::set<types::global_dof_index> extra_halo;
864 for (
const auto dof : dofs_with_support_on_halo_cells)
870 const unsigned int line_size = line_ptr->size();
871 for (
unsigned int j = 0; j < line_size; ++j)
872 extra_halo.insert((*line_ptr)[j].first);
875 dofs_with_support_on_halo_cells.insert(extra_halo.begin(),
882 support_set.
add_indices(predicate_dofs.begin(), predicate_dofs.end());
886 halo_set.
add_indices(dofs_with_support_on_halo_cells.begin(),
887 dofs_with_support_on_halo_cells.end());
902 template <
int spacedim>
911 template <
int spacedim>
915 const unsigned int dim = 2;
924 if (!cell->is_artificial())
926 for (
const unsigned int face : cell->face_indices())
927 if (cell->face(face)->has_children())
930 line = cell->face(face);
934 selected_dofs.add_index(
935 line->child(0)->vertex_dof_index(1, dof));
937 for (
unsigned int child = 0; child < 2; ++child)
939 if (cell->neighbor_child_on_subface(face, child)
944 selected_dofs.add_index(
945 line->child(child)->dof_index(dof));
950 selected_dofs.compress();
951 return selected_dofs;
955 template <
int spacedim>
959 const unsigned int dim = 3;
967 if (!cell->is_artificial())
968 for (
auto f : cell->face_indices())
972 if (cell->face(f)->has_children())
974 for (
unsigned int child = 0; child < 4; ++child)
975 if (!cell->neighbor_child_on_subface(f, child)
979 std::vector<types::global_dof_index> ldi(
981 face->child(child)->get_dof_indices(ldi);
982 selected_dofs.add_indices(ldi.begin(), ldi.end());
987 for (
unsigned int vertex = 0; vertex < 4; ++vertex)
990 unconstrained_dofs.add_index(
991 face->vertex_dof_index(vertex, dof));
994 selected_dofs.subtract_set(unconstrained_dofs);
995 return selected_dofs;
1002 template <
int dim,
int spacedim>
1011 template <
int dim,
int spacedim>
1015 std::vector<bool> &selected_dofs)
1021 std::fill_n(selected_dofs.begin(), dof_handler.
n_dofs(),
false);
1023 std::vector<types::global_dof_index> local_dof_indices;
1024 local_dof_indices.reserve(
1032 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1033 local_dof_indices.resize(dofs_per_cell);
1034 cell->get_dof_indices(local_dof_indices);
1035 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1036 selected_dofs[local_dof_indices[i]] =
true;
1042 template <
int dim,
int spacedim>
1052 std::vector<types::global_dof_index> dof_indices;
1053 std::set<types::global_dof_index> global_dof_indices;
1058 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1059 cell->get_dof_indices(dof_indices);
1063 global_dof_indices.insert(dof_index);
1066 dof_set.
add_indices(global_dof_indices.begin(), global_dof_indices.end());
1075 template <
int dim,
int spacedim>
1085 template <
int dim,
int spacedim>
1089 const unsigned int level)
1097 std::vector<types::global_dof_index> dof_indices;
1098 std::set<types::global_dof_index> global_dof_indices;
1100 const auto filtered_iterators_range =
1103 for (
const auto &cell : filtered_iterators_range)
1105 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1106 cell->get_mg_dof_indices(dof_indices);
1110 global_dof_indices.insert(dof_index);
1113 dof_set.
add_indices(global_dof_indices.begin(), global_dof_indices.end());
1122 template <
int dim,
int spacedim>
1127 const unsigned int level)
1134 template <
int dim,
int spacedim>
1149 std::vector<types::global_dof_index> dof_indices;
1150 std::vector<types::global_dof_index> dofs_on_ghosts;
1153 if (cell->is_ghost())
1155 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1156 cell->get_dof_indices(dof_indices);
1157 for (
const auto dof_index : dof_indices)
1159 dofs_on_ghosts.push_back(dof_index);
1163 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1164 dof_set.
add_indices(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1172 template <
int dim,
int spacedim>
1182 template <
int dim,
int spacedim>
1186 const unsigned int level)
1199 std::vector<types::global_dof_index> dof_indices;
1200 std::vector<types::global_dof_index> dofs_on_ghosts;
1211 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1212 cell->get_mg_dof_indices(dof_indices);
1213 for (
const auto dof_index : dof_indices)
1215 dofs_on_ghosts.push_back(dof_index);
1219 std::sort(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1220 dof_set.
add_indices(dofs_on_ghosts.begin(), dofs_on_ghosts.end());
1228 template <
int dim,
int spacedim>
1232 const unsigned int level,
1240 template <
int dim,
int spacedim>
1244 std::vector<std::vector<bool>> &constant_modes)
1250 constant_modes = std::vector<std::vector<bool>>(0);
1258 std::vector<unsigned char> dofs_by_component(
1263 unsigned int n_selected_dofs = 0;
1264 for (
unsigned int i = 0; i < n_components; ++i)
1265 if (component_mask[i] ==
true)
1267 std::count(dofs_by_component.begin(), dofs_by_component.end(), i);
1271 std::vector<unsigned int> component_numbering(
1273 for (
unsigned int i = 0, count = 0; i < locally_owned_dofs.
n_elements();
1275 if (component_mask[dofs_by_component[i]])
1276 component_numbering[i] = count++;
1283 const ::hp::FECollection<dim, spacedim> &fe_collection =
1285 std::vector<Table<2, bool>> element_constant_modes;
1286 std::vector<std::vector<std::pair<unsigned int, unsigned int>>>
1287 constant_mode_to_component_translation(n_components);
1289 unsigned int n_constant_modes = 0;
1290 int first_non_empty_constant_mode = -1;
1291 for (
unsigned int f = 0; f < fe_collection.size(); ++f)
1293 std::pair<Table<2, bool>, std::vector<unsigned int>> data =
1294 fe_collection[f].get_constant_modes();
1298 if (first_non_empty_constant_mode < 0 && data.first.n_rows() > 0)
1300 first_non_empty_constant_mode = f;
1304 for (
unsigned int i = 0; i < data.second.size(); ++i)
1305 if (component_mask[data.second[i]])
1306 constant_mode_to_component_translation[data.second[i]]
1307 .emplace_back(n_constant_modes++, i);
1313 element_constant_modes.push_back(data.first);
1315 element_constant_modes.back().n_rows() == 0 ||
1316 element_constant_modes.back().n_rows() ==
1317 element_constant_modes[first_non_empty_constant_mode].n_rows(),
1324 constant_modes.clear();
1325 constant_modes.resize(n_constant_modes,
1326 std::vector<bool>(n_selected_dofs,
false));
1330 std::vector<types::global_dof_index> dof_indices;
1334 dof_indices.resize(cell->get_fe().n_dofs_per_cell());
1335 cell->get_dof_indices(dof_indices);
1337 for (
unsigned int i = 0; i < dof_indices.size(); ++i)
1338 if (locally_owned_dofs.
is_element(dof_indices[i]))
1340 const unsigned int loc_index =
1342 const unsigned int comp = dofs_by_component[loc_index];
1343 if (component_mask[comp])
1344 for (
auto &indices :
1345 constant_mode_to_component_translation[comp])
1346 constant_modes[indices
1347 .first][component_numbering[loc_index]] =
1348 element_constant_modes[cell->active_fe_index()](
1356 template <
int dim,
int spacedim>
1357 std::map<
typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
1358 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1369 std::map<
typename DoFHandler<dim - 1, spacedim>::active_cell_iterator,
1370 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1372 c1_to_c0_cell_and_face;
1375 if (c1_to_c0_face.empty())
1376 return c1_to_c0_cell_and_face;
1379 std::map<typename Triangulation<dim, spacedim>::face_iterator,
1380 typename DoFHandler<dim - 1, spacedim>::active_cell_iterator>
1384 for (
const auto &[c1_cell, c0_cell] : c1_to_c0_face)
1385 if (!c1_cell->has_children())
1386 c0_to_c1[c0_cell] = c1_cell->as_dof_handler_iterator(c1_dh);
1390 for (
const auto &cell :
1392 for (
const auto f : cell->face_indices())
1393 if (cell->face(f)->at_boundary())
1395 const auto &it = c0_to_c1.find(cell->face(f));
1396 if (it != c0_to_c1.end())
1398 const auto &c1_cell = it->second;
1399 c1_to_c0_cell_and_face[c1_cell] = {cell, f};
1405 return c1_to_c0_cell_and_face;
1410 template <
int dim,
int spacedim>
1413 std::vector<unsigned int> &active_fe_indices)
1420 active_fe_indices.assign(indices.begin(), indices.end());
1423 template <
int dim,
int spacedim>
1424 std::vector<IndexSet>
1428 ExcMessage(
"The given DoFHandler has no DoFs."));
1436 "For parallel::distributed::Triangulation objects and "
1437 "associated DoF handler objects, asking for any information "
1438 "related to a subdomain other than the locally owned one does "
1439 "not make sense."));
1443 std::vector<::types::subdomain_id> subdomain_association(
1446 subdomain_association);
1460 const unsigned int n_subdomains =
1464 unsigned int max_subdomain_id = 0;
1467 std::max(max_subdomain_id, cell->subdomain_id());
1468 return max_subdomain_id + 1;
1473 ->get_communicator()));
1474 Assert(n_subdomains > *std::max_element(subdomain_association.begin(),
1475 subdomain_association.end()),
1478 std::vector<::IndexSet> index_sets(
1487 index < subdomain_association.size();
1491 if (subdomain_association[index] != this_subdomain)
1493 index_sets[this_subdomain].add_range(i_min, index);
1495 this_subdomain = subdomain_association[index];
1500 if (i_min == subdomain_association.size() - 1)
1502 index_sets[this_subdomain].add_index(i_min);
1508 index_sets[this_subdomain].add_range(i_min,
1509 subdomain_association.size());
1512 for (
unsigned int i = 0; i < n_subdomains; ++i)
1518 template <
int dim,
int spacedim>
1519 std::vector<IndexSet>
1529 "For parallel::distributed::Triangulation objects and "
1530 "associated DoF handler objects, asking for any information "
1531 "related to a subdomain other than the locally owned one does "
1532 "not make sense."));
1540 std::vector<IndexSet> dof_set =
1561 std::vector<types::global_dof_index> local_dof_indices;
1562 std::set<types::global_dof_index> subdomain_halo_global_dof_indices;
1563 for (
typename std::vector<
1565 const_iterator it_cell = active_halo_layer.
begin();
1566 it_cell != active_halo_layer.end();
1574 "The subdomain ID of the halo cell should not match that of the vector entry."));
1577 cell->get_dof_indices(local_dof_indices);
1581 subdomain_halo_global_dof_indices.insert(local_dof_index);
1585 subdomain_halo_global_dof_indices.begin(),
1586 subdomain_halo_global_dof_indices.end());
1594 template <
int dim,
int spacedim>
1598 std::vector<types::subdomain_id> &subdomain_association)
1606 "For parallel::distributed::Triangulation objects and "
1607 "associated DoF handler objects, asking for any subdomain other "
1608 "than the locally owned one does not make sense."));
1610 Assert(subdomain_association.size() == dof_handler.
n_dofs(),
1624 std::vector<types::subdomain_id> cell_owners(
1630 cell_owners = tr->get_true_subdomain_ids_of_cells();
1631 Assert(tr->get_true_subdomain_ids_of_cells().size() ==
1632 tr->n_active_cells(),
1639 cell_owners[cell->active_cell_index()] = cell->subdomain_id();
1643 std::fill_n(subdomain_association.begin(),
1647 std::vector<types::global_dof_index> local_dof_indices;
1648 local_dof_indices.reserve(
1656 cell_owners[cell->active_cell_index()];
1657 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1658 local_dof_indices.resize(dofs_per_cell);
1659 cell->get_dof_indices(local_dof_indices);
1665 for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1666 if (subdomain_association[local_dof_indices[i]] ==
1668 subdomain_association[local_dof_indices[i]] =
subdomain_id;
1669 else if (subdomain_association[local_dof_indices[i]] >
subdomain_id)
1671 subdomain_association[local_dof_indices[i]] =
subdomain_id;
1675 Assert(std::find(subdomain_association.begin(),
1676 subdomain_association.end(),
1678 subdomain_association.end(),
1684 template <
int dim,
int spacedim>
1690 std::vector<types::subdomain_id> subdomain_association(
1694 return std::count(subdomain_association.begin(),
1695 subdomain_association.end(),
1701 template <
int dim,
int spacedim>
1714 "For parallel::distributed::Triangulation objects and "
1715 "associated DoF handler objects, asking for any subdomain other "
1716 "than the locally owned one does not make sense."));
1720 std::vector<types::global_dof_index> local_dof_indices;
1721 local_dof_indices.reserve(
1729 std::vector<types::global_dof_index> subdomain_indices;
1732 if ((cell->is_artificial() ==
false) &&
1733 (cell->subdomain_id() == subdomain))
1735 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1736 local_dof_indices.resize(dofs_per_cell);
1737 cell->get_dof_indices(local_dof_indices);
1738 subdomain_indices.insert(subdomain_indices.end(),
1739 local_dof_indices.begin(),
1740 local_dof_indices.end());
1743 std::sort(subdomain_indices.begin(), subdomain_indices.end());
1744 index_set.
add_indices(subdomain_indices.begin(), subdomain_indices.end());
1752 template <
int dim,
int spacedim>
1757 std::vector<unsigned int> &n_dofs_on_subdomain)
1762 std::fill(n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0);
1771 return cell.subdomain_id() == subdomain;
1773 ExcMessage(
"There are no cells for the given subdomain!"));
1775 std::vector<types::subdomain_id> subdomain_association(
1779 std::vector<unsigned char> component_association(dof_handler.
n_dofs());
1782 component_association);
1784 for (
unsigned int c = 0; c < dof_handler.
get_fe(0).n_components(); ++c)
1787 if ((subdomain_association[i] == subdomain) &&
1788 (component_association[i] ==
static_cast<unsigned char>(c)))
1789 ++n_dofs_on_subdomain[c];
1801 template <
int dim,
int spacedim>
1804 const std::vector<unsigned char> &dofs_by_component,
1805 const std::vector<unsigned int> &target_component,
1806 const bool only_once,
1807 std::vector<types::global_dof_index> &dofs_per_component,
1808 unsigned int &component)
1827 for (
unsigned int dd = 0; dd <
d; ++dd, ++component)
1828 dofs_per_component[target_component[component]] +=
1829 std::count(dofs_by_component.begin(),
1830 dofs_by_component.end(),
1837 for (
unsigned int dd = 1; dd <
d; ++dd)
1838 dofs_per_component[target_component[component -
d + dd]] =
1839 dofs_per_component[target_component[component -
d]];
1846 template <
int dim,
int spacedim>
1849 const std::vector<unsigned char> &dofs_by_component,
1850 const std::vector<unsigned int> &target_component,
1851 const bool only_once,
1852 std::vector<types::global_dof_index> &dofs_per_component,
1853 unsigned int &component)
1858 for (
unsigned int fe = 1; fe < fe_collection.
size(); ++fe)
1860 Assert(fe_collection[fe].n_components() ==
1861 fe_collection[0].n_components(),
1863 Assert(fe_collection[fe].n_base_elements() ==
1864 fe_collection[0].n_base_elements(),
1866 for (
unsigned int b = 0;
b < fe_collection[0].n_base_elements(); ++
b)
1868 Assert(fe_collection[fe].base_element(
b).n_components() ==
1869 fe_collection[0].base_element(
b).n_components(),
1871 Assert(fe_collection[fe].base_element(
b).n_base_elements() ==
1872 fe_collection[0].base_element(
b).n_base_elements(),
1895 template <
int dim,
int spacedim>
1907 template <
int dim,
int spacedim>
1909 all_elements_are_primitive(
1910 const ::hp::FECollection<dim, spacedim> &fe_collection)
1912 for (
unsigned int i = 0; i < fe_collection.size(); ++i)
1913 if (fe_collection[i].is_primitive() ==
false)
1923 template <
int dim,
int spacedim>
1924 std::vector<types::global_dof_index>
1927 const bool only_once,
1928 const std::vector<unsigned int> &target_component_)
1934 std::vector<unsigned int> target_component = target_component_;
1935 if (target_component.empty())
1937 target_component.resize(n_components);
1938 for (
unsigned int i = 0; i < n_components; ++i)
1939 target_component[i] = i;
1942 Assert(target_component.size() == n_components,
1946 const unsigned int max_component =
1947 *std::max_element(target_component.begin(), target_component.end());
1948 const unsigned int n_target_components = max_component + 1;
1950 std::vector<types::global_dof_index> dofs_per_component(
1955 if (n_components == 1)
1958 return dofs_per_component;
1964 std::vector<unsigned char> dofs_by_component(
1971 unsigned int component = 0;
1982 Assert((internal::all_elements_are_primitive(
1984 (std::accumulate(dofs_per_component.begin(),
1985 dofs_per_component.end(),
1991 #ifdef DEAL_II_WITH_MPI
1997 std::vector<types::global_dof_index> local_dof_count =
2000 const int ierr = MPI_Allreduce(local_dof_count.data(),
2001 dofs_per_component.data(),
2002 n_target_components,
2010 return dofs_per_component;
2015 template <
int dim,
int spacedim>
2016 std::vector<types::global_dof_index>
2018 const std::vector<unsigned int> &target_block_)
2020 const ::hp::FECollection<dim, spacedim> &fe_collection =
2029 const unsigned int n_blocks = fe_collection[0].n_blocks();
2031 std::vector<unsigned int> target_block = target_block_;
2032 if (target_block.empty())
2034 target_block.resize(fe_collection[0].
n_blocks());
2035 for (
unsigned int i = 0; i <
n_blocks; ++i)
2036 target_block[i] = i;
2041 for (
unsigned int f = 1; f < fe_collection.size(); ++f)
2043 ExcMessage(
"This function can only work if all elements in a "
2044 "collection have the same number of blocks."));
2050 std::vector<types::global_dof_index> dofs_per_block(1);
2051 dofs_per_block[0] = dof_handler.
n_dofs();
2052 return dofs_per_block;
2056 const unsigned int max_block =
2057 *std::max_element(target_block.begin(), target_block.end());
2058 const unsigned int n_target_blocks = max_block + 1;
2060 std::vector<types::global_dof_index> dofs_per_block(n_target_blocks);
2064 for (
unsigned int this_fe = fe_collection.size() - 1;
2065 this_fe < fe_collection.size();
2070 std::vector<unsigned char> dofs_by_block(
2075 for (
unsigned int block = 0; block < fe.
n_blocks(); ++block)
2076 dofs_per_block[target_block[block]] +=
2077 std::count(dofs_by_block.begin(), dofs_by_block.end(), block);
2079 #ifdef DEAL_II_WITH_MPI
2086 std::vector<types::global_dof_index> local_dof_count =
2088 const int ierr = MPI_Allreduce(local_dof_count.data(),
2089 dofs_per_block.data(),
2099 return dofs_per_block;
2104 template <
int dim,
int spacedim>
2107 std::vector<types::global_dof_index> &mapping)
2110 mapping.insert(mapping.end(),
2114 std::vector<types::global_dof_index> dofs_on_face;
2125 for (
const unsigned int f : cell->face_indices())
2126 if (cell->at_boundary(f))
2128 const unsigned int dofs_per_face =
2129 cell->get_fe().n_dofs_per_face(f);
2130 dofs_on_face.resize(dofs_per_face);
2131 cell->face(f)->get_dof_indices(dofs_on_face,
2132 cell->active_fe_index());
2133 for (
unsigned int i = 0; i < dofs_per_face; ++i)
2135 mapping[dofs_on_face[i]] = next_boundary_index++;
2143 template <
int dim,
int spacedim>
2146 const std::set<types::boundary_id> &boundary_ids,
2147 std::vector<types::global_dof_index> &mapping)
2154 mapping.insert(mapping.end(),
2159 if (boundary_ids.empty())
2162 std::vector<types::global_dof_index> dofs_on_face;
2167 for (
const unsigned int f : cell->face_indices())
2168 if (boundary_ids.find(cell->face(f)->boundary_id()) !=
2171 const unsigned int dofs_per_face =
2172 cell->get_fe().n_dofs_per_face(f);
2173 dofs_on_face.resize(dofs_per_face);
2174 cell->face(f)->get_dof_indices(dofs_on_face,
2175 cell->active_fe_index());
2176 for (
unsigned int i = 0; i < dofs_per_face; ++i)
2178 mapping[dofs_on_face[i]] = next_boundary_index++;
2189 template <
int dim,
int spacedim>
2190 std::map<types::global_dof_index, Point<spacedim>>
2196 std::map<types::global_dof_index, Point<spacedim>> support_points;
2207 (fe_collection[
fe_index].has_support_points()),
2210 fe_collection[
fe_index].get_unit_support_points()));
2215 (in_mask.
size() == 0 ?
2231 std::vector<types::global_dof_index> local_dof_indices;
2234 if (cell->is_artificial() ==
false)
2236 hp_fe_values.reinit(cell);
2240 local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
2241 cell->get_dof_indices(local_dof_indices);
2243 const std::vector<Point<spacedim>> &points =
2245 for (
unsigned int i = 0; i < cell->get_fe().n_dofs_per_cell();
2248 const unsigned int dof_comp =
2249 cell->get_fe().system_to_component_index(i).first;
2253 support_points[local_dof_indices[i]] = points[i];
2257 return support_points;
2261 template <
int dim,
int spacedim>
2262 std::vector<Point<spacedim>>
2263 map_dofs_to_support_points_vector(
2268 std::vector<Point<spacedim>> support_points(dof_handler.
n_dofs());
2271 const std::map<types::global_dof_index, Point<spacedim>>
2279 Assert(x_support_points.find(i) != x_support_points.end(),
2282 support_points[i] = x_support_points.find(i)->second;
2285 return support_points;
2291 template <
int dim,
int spacedim>
2303 "This function can not be used with distributed triangulations. "
2304 "See the documentation for more information."));
2311 internal::map_dofs_to_support_points_vector(mapping_collection,
2317 template <
int dim,
int spacedim>
2330 "This function can not be used with distributed triangulations. "
2331 "See the documentation for more information."));
2336 internal::map_dofs_to_support_points_vector(mapping, dof_handler, mask);
2341 template <
int dim,
int spacedim>
2349 support_points.clear();
2362 template <
int dim,
int spacedim>
2370 support_points.clear();
2379 template <
int dim,
int spacedim>
2380 std::map<types::global_dof_index, Point<spacedim>>
2395 template <
int dim,
int spacedim>
2396 std::map<types::global_dof_index, Point<spacedim>>
2406 template <
int spacedim>
2414 std::map<Point<spacedim>,
2415 std::vector<types::global_dof_index>,
2420 for (
const auto &it : support_points)
2422 std::vector<types::global_dof_index> &v = point_map[it.second];
2423 v.push_back(it.first);
2427 for (
const auto &it : point_map)
2429 out << it.first <<
" \"";
2430 const std::vector<types::global_dof_index> &v = it.second;
2431 for (
unsigned int i = 0; i < v.size(); ++i)
2445 template <
int dim,
int spacedim>
2454 const unsigned int nb = fe.
n_blocks();
2456 tables_by_block.resize(1);
2457 tables_by_block[0].reinit(nb, nb);
2458 tables_by_block[0].fill(
none);
2466 tables_by_block[0](ib, jb) |= table(i, j);
2474 tables_by_block.resize(fe_collection.
size());
2476 for (
unsigned int f = 0; f < fe_collection.
size(); ++f)
2480 const unsigned int nb = fe.
n_blocks();
2481 tables_by_block[f].reinit(nb, nb);
2482 tables_by_block[f].fill(
none);
2489 tables_by_block[f](ib, jb) |= table(i, j);
2498 template <
int dim,
int spacedim>
2502 const unsigned int level,
2503 const std::vector<bool> &selected_dofs,
2506 std::vector<types::global_dof_index> indices;
2511 if (cell->is_locally_owned_on_level())
2518 if (cell->is_locally_owned_on_level())
2520 indices.resize(cell->get_fe().n_dofs_per_cell());
2521 cell->get_mg_dof_indices(indices);
2523 if (selected_dofs.size() != 0)
2528 if (selected_dofs.empty())
2529 block_list.
add(i, indices[j] - offset);
2532 if (selected_dofs[j])
2533 block_list.
add(i, indices[j] - offset);
2541 template <
int dim,
int spacedim>
2545 const unsigned int level,
2546 const bool interior_only)
2551 std::vector<types::global_dof_index> indices;
2552 std::vector<bool> exclude;
2556 indices.resize(cell->get_fe().n_dofs_per_cell());
2557 cell->get_mg_dof_indices(indices);
2563 std::fill(exclude.begin(), exclude.end(),
false);
2565 for (
const unsigned int face : cell->face_indices())
2566 if (cell->at_boundary(face) ||
2567 cell->neighbor(face)->level() != cell->level())
2572 block_list.
add(0, indices[j]);
2576 for (
const auto index : indices)
2577 block_list.
add(0, index);
2583 template <
int dim,
int spacedim>
2587 const unsigned int level,
2588 const bool interior_dofs_only,
2589 const bool boundary_dofs)
2594 std::vector<types::global_dof_index> indices;
2595 std::vector<bool> exclude;
2597 unsigned int block = 0;
2600 if (pcell->is_active())
2603 for (
unsigned int child = 0; child < pcell->n_children(); ++child)
2605 const auto cell = pcell->child(child);
2610 indices.resize(n_dofs);
2611 exclude.resize(n_dofs);
2612 std::fill(exclude.begin(), exclude.end(),
false);
2613 cell->get_mg_dof_indices(indices);
2615 if (interior_dofs_only)
2619 for (
unsigned int d = 0;
d < dim; ++
d)
2621 const unsigned int face =
2630 for (
const unsigned int face :
2632 if (cell->at_boundary(face))
2638 for (
unsigned int i = 0; i < n_dofs; ++i)
2640 block_list.
add(block, indices[i]);
2646 template <
int dim,
int spacedim>
2647 std::vector<unsigned int>
2650 const unsigned int level,
2651 const bool interior_only,
2652 const bool boundary_patches,
2653 const bool level_boundary_patches,
2654 const bool single_cell_patches,
2655 const bool invert_vertex_mapping)
2662 exclude_boundary_dofs,
2664 level_boundary_patches,
2665 single_cell_patches,
2666 invert_vertex_mapping);
2669 template <
int dim,
int spacedim>
2670 std::vector<unsigned int>
2673 const unsigned int level,
2675 const bool boundary_patches,
2676 const bool level_boundary_patches,
2677 const bool single_cell_patches,
2678 const bool invert_vertex_mapping)
2682 std::vector<unsigned int> vertex_cell_count(
2686 std::vector<bool> vertex_boundary(
2689 std::vector<unsigned int> vertex_mapping(
2694 std::vector<unsigned int> vertex_dof_count(
2700 for (
const unsigned int v : cell->vertex_indices())
2702 const unsigned int vg = cell->vertex_index(v);
2703 vertex_dof_count[vg] += cell->get_fe().n_dofs_per_cell();
2704 ++vertex_cell_count[vg];
2705 for (
unsigned int d = 0;
d < dim; ++
d)
2708 if (cell->at_boundary(face))
2709 vertex_boundary[vg] =
true;
2710 else if ((!level_boundary_patches) &&
2711 (cell->neighbor(face)->level() !=
2712 static_cast<int>(
level)))
2713 vertex_boundary[vg] =
true;
2719 for (
unsigned int vg = 0; vg < vertex_dof_count.size(); ++vg)
2720 if ((!single_cell_patches && vertex_cell_count[vg] < 2) ||
2721 (!boundary_patches && vertex_boundary[vg]))
2722 vertex_dof_count[vg] = 0;
2725 unsigned int n_vertex_count = 0;
2726 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2727 if (vertex_dof_count[vg] != 0)
2728 vertex_mapping[vg] = n_vertex_count++;
2731 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2732 if (vertex_dof_count[vg] != 0)
2733 vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
2737 vertex_dof_count.resize(n_vertex_count);
2741 block_list.
reinit(vertex_dof_count.size(),
2745 std::vector<types::global_dof_index> indices;
2746 std::vector<bool> exclude;
2752 cell->get_mg_dof_indices(indices);
2754 for (
const unsigned int v : cell->vertex_indices())
2756 const unsigned int vg = cell->vertex_index(v);
2757 const unsigned int block = vertex_mapping[vg];
2763 if (exclude_boundary_dofs.
size() == 0 ||
2769 std::fill(exclude.begin(), exclude.end(),
false);
2771 for (
unsigned int d = 0;
d < dim; ++
d)
2773 const unsigned int a_face =
2775 const unsigned int face =
2788 for (
unsigned int j = 0; j < indices.size(); ++j)
2790 block_list.
add(block, indices[j]);
2794 for (
const auto index : indices)
2795 block_list.
add(block, index);
2800 if (invert_vertex_mapping)
2803 unsigned int n_vertex_count = 0;
2804 for (
unsigned int vg = 0; vg < vertex_mapping.size(); ++vg)
2806 vertex_mapping[n_vertex_count++] = vg;
2809 vertex_mapping.resize(n_vertex_count);
2812 return vertex_mapping;
2816 template <
int dim,
int spacedim>
2822 std::set<types::global_dof_index> dofs_on_patch;
2823 std::vector<types::global_dof_index> local_dof_indices;
2828 for (
unsigned int i = 0; i < patch.size(); ++i)
2832 Assert(cell->is_artificial() ==
false,
2833 ExcMessage(
"This function can not be called with cells that are "
2834 "not either locally owned or ghost cells."));
2836 cell->get_dof_indices(local_dof_indices);
2837 dofs_on_patch.insert(local_dof_indices.begin(),
2838 local_dof_indices.end());
2842 return dofs_on_patch.size();
2847 template <
int dim,
int spacedim>
2848 std::vector<types::global_dof_index>
2853 std::set<types::global_dof_index> dofs_on_patch;
2854 std::vector<types::global_dof_index> local_dof_indices;
2859 for (
unsigned int i = 0; i < patch.size(); ++i)
2863 Assert(cell->is_artificial() ==
false,
2864 ExcMessage(
"This function can not be called with cells that are "
2865 "not either locally owned or ghost cells."));
2867 cell->get_dof_indices(local_dof_indices);
2868 dofs_on_patch.insert(local_dof_indices.begin(),
2869 local_dof_indices.end());
2872 Assert((dofs_on_patch.size() == count_dofs_on_patch<dim, spacedim>(patch)),
2879 return std::vector<types::global_dof_index>(dofs_on_patch.begin(),
2880 dofs_on_patch.end());
2890 #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
std::vector< types::fe_index > get_active_fe_indices() 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 IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) 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 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)
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
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)
@ update_quadrature_points
Transformed quadrature points.
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > 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
unsigned short int fe_index
const ::Triangulation< dim, spacedim > & tria
#define DEAL_II_DOF_INDEX_MPI_TYPE