54 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
86 [](
bool p,
bool q) { return !p || q; }),
88 "marked_vertices should be a subset of used vertices in the triangulation "
89 "but marked_vertices contains one or more vertices that are not used vertices!"));
97 const std::vector<bool> &
used =
102 std::vector<bool>::const_iterator
first =
103 std::find(
used.begin(),
used.end(),
true);
117 double dist = (
p -
vertices[j]).norm_square();
130 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
134 const MeshType<dim, spacedim> &
mesh,
139 if (
mapping.preserves_vertex_locations() ==
true)
167 [](
bool p,
bool q) { return !p || q; }),
169 "marked_vertices should be a subset of used vertices in the triangulation "
170 "but marked_vertices contains one or more vertices that are not used vertices!"));
191 template <
typename MeshType>
194 std::vector<typename MeshType::active_cell_iterator>
197 typename ::internal::ActiveCellIterator<MeshType::dimension,
198 MeshType::space_dimension,
202 const unsigned int vertex)
204 const int dim = MeshType::dimension;
205 const int spacedim = MeshType::space_dimension;
211 Assert(
mesh.get_triangulation().get_used_vertices()[vertex],
217 std::set<typename ::internal::
218 ActiveCellIterator<dim, spacedim, MeshType>::type>
221 typename ::internal::ActiveCellIterator<dim, spacedim, MeshType>::type
222 cell =
mesh.begin_active(),
258 for (; cell != endc; ++cell)
260 for (
const unsigned int v : cell->vertex_indices())
261 if (cell->vertex_index(v) == vertex)
273 for (
const auto face :
274 cell->reference_cell().faces_for_given_vertex(v))
275 if (!cell->at_boundary(face) &&
276 cell->neighbor(face)->is_active())
298 for (
unsigned int e = 0; e < cell->n_lines(); ++e)
299 if (cell->line(e)->has_children())
303 if (cell->line(e)->child(0)->vertex_index(1) == vertex)
326 return std::vector<typename ::internal::
327 ActiveCellIterator<dim, spacedim, MeshType>::type>(
335 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
338 void find_active_cell_around_point_internal(
339 const MeshType<dim, spacedim> &
mesh,
341 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>
343 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>
347 typename ::internal::
348 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
352 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
357 using cell_iterator =
358 typename MeshType<dim, spacedim>::active_cell_iterator;
360 using cell_iterator = typename ::internal::
361 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
369 std::set<cell_iterator> adjacent_cells_new;
373 std::vector<cell_iterator> active_neighbors;
374 get_active_neighbors<MeshType<dim, spacedim>>(cell, active_neighbors);
375 for (
unsigned int i = 0; i < active_neighbors.size(); ++i)
376 if (searched_cells.find(active_neighbors[i]) ==
377 searched_cells.end())
378 adjacent_cells_new.insert(active_neighbors[i]);
382 adjacent_cells_new.end());
391 cell_iterator
it =
mesh.begin_active();
393 if (searched_cells.find(
it) == searched_cells.end())
404 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
408 typename MeshType<dim, spacedim>::active_cell_iterator
410 typename ::internal::
411 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
416 const double tolerance)
418 return find_active_cell_around_point<dim, MeshType, spacedim>(
429 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
433 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
435 std::pair<typename ::internal::
436 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
440 const MeshType<dim, spacedim> &
mesh,
443 const double tolerance)
445 using active_cell_iterator = typename ::internal::
446 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
452 double best_distance = tolerance;
454 std::pair<active_cell_iterator, Point<dim>> best_cell;
457 best_cell.first =
mesh.end();
461 std::vector<active_cell_iterator> adjacent_cells_tmp =
470 std::set<active_cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
471 adjacent_cells_tmp.end());
472 std::set<active_cell_iterator> searched_cells;
480 const auto n_active_cells =
mesh.get_triangulation().n_active_cells();
482 unsigned int cells_searched = 0;
483 while (!found && cells_searched < n_active_cells)
487 if (cell->is_artificial() ==
false)
492 bool any_vertex_marked =
false;
497 any_vertex_marked =
true;
501 if (!any_vertex_marked)
508 mapping.transform_real_to_unit_cell(cell,
p);
513 cell->reference_cell().closest_point(p_cell).
distance(
520 if ((dist < best_distance) ||
521 ((dist == best_distance) &&
522 (cell->level() > best_level)))
525 best_distance = dist;
526 best_level = cell->level();
527 best_cell = std::make_pair(cell, p_cell);
558 if (!found && cells_searched < n_active_cells)
560 find_active_cell_around_point_internal<dim, MeshType, spacedim>(
570 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
574 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
577 std::vector<std::pair<
578 typename ::internal::
579 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
583 const MeshType<dim, spacedim> &
mesh,
585 const double tolerance,
591 if (cell_and_point.first ==
mesh.end())
594 return find_all_active_cells_around_point(
600 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
604 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
607 std::vector<std::pair<
608 typename ::internal::
609 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
612 find_all_active_cells_around_point(
614 const MeshType<dim, spacedim> &
mesh,
616 const double tolerance,
617 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
620 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
624 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
629 cells_and_points.push_back(first_cell);
631 const Point<dim> unit_point = cells_and_points.front().second;
632 const auto my_cell = cells_and_points.front().first;
634 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
637 if (my_cell->reference_cell().is_hyper_cube())
643 unsigned int n_dirs_at_threshold = 0;
645 for (
unsigned int d = 0;
d < dim; ++
d)
647 distance_to_center[
d] =
std::abs(unit_point[d] - 0.5);
648 if (distance_to_center[d] > 0.5 - tolerance)
650 ++n_dirs_at_threshold;
651 last_point_at_threshold =
d;
656 if (n_dirs_at_threshold == 1)
658 unsigned int neighbor_index =
659 2 * last_point_at_threshold +
660 (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
661 if (!my_cell->at_boundary(neighbor_index))
663 const auto neighbor_cell = my_cell->neighbor(neighbor_index);
665 if (neighbor_cell->is_active())
666 cells_to_add.push_back(neighbor_cell);
668 for (
const auto &child_cell :
669 neighbor_cell->child_iterators())
671 if (child_cell->is_active())
672 cells_to_add.push_back(child_cell);
677 else if (n_dirs_at_threshold == dim)
679 unsigned int local_vertex_index = 0;
680 for (
unsigned int d = 0;
d < dim; ++
d)
681 local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) <<
d;
683 const auto fu = [&](
const auto &tentative_cells) {
684 for (
const auto &cell : tentative_cells)
686 cells_to_add.push_back(cell);
689 const auto vertex_index = my_cell->vertex_index(local_vertex_index);
691 if (vertex_to_cells !=
nullptr)
692 fu((*vertex_to_cells)[vertex_index]);
701 else if (n_dirs_at_threshold == 2)
704 unsigned int count_vertex_indices = 0;
706 for (
unsigned int d = 0;
d < dim; ++
d)
708 if (distance_to_center[d] > 0.5 - tolerance)
712 unit_point[
d] > 0.5 ? 1 : 0;
713 ++count_vertex_indices;
723 const unsigned int first_vertex =
726 for (
unsigned int d = 0;
d < 2; ++
d)
728 const auto fu = [&](
const auto &tentative_cells) {
729 for (
const auto &cell : tentative_cells)
731 bool cell_not_yet_present =
true;
732 for (
const auto &other_cell : cells_to_add)
733 if (cell == other_cell)
735 cell_not_yet_present =
false;
738 if (cell_not_yet_present)
739 cells_to_add.push_back(cell);
743 const auto vertex_index =
744 my_cell->vertex_index(first_vertex + (d << free_direction));
746 if (vertex_to_cells !=
nullptr)
747 fu((*vertex_to_cells)[vertex_index]);
763 const auto fu = [&](
const auto &tentative_cells) {
764 for (
const auto &cell : tentative_cells)
766 bool cell_not_yet_present =
true;
767 for (
const auto &other_cell : cells_to_add)
768 if (cell == other_cell)
770 cell_not_yet_present =
false;
773 if (cell_not_yet_present)
774 cells_to_add.push_back(cell);
778 const auto vertex_index = my_cell->vertex_index(v);
780 if (vertex_to_cells !=
nullptr)
781 fu((*vertex_to_cells)[vertex_index]);
787 for (
const auto &cell : cells_to_add)
793 mapping.transform_real_to_unit_cell(cell,
p);
794 if (cell->reference_cell().contains_point(p_unit, tolerance))
795 cells_and_points.emplace_back(cell, p_unit);
802 cells_and_points.begin(),
803 cells_and_points.end(),
804 [](
const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
806 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
807 Point<dim>> &b) { return a.first < b.first; });
809 return cells_and_points;
814 template <
typename MeshType>
818 const MeshType &
mesh,
819 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
822 std::vector<typename MeshType::active_cell_iterator> active_halo_layer;
823 std::vector<bool> locally_active_vertices_on_subdomain(
824 mesh.get_triangulation().n_vertices(),
false);
826 std::map<unsigned int, std::vector<unsigned int>> coinciding_vertex_groups;
827 std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
829 coinciding_vertex_groups,
830 vertex_to_coinciding_vertex_group);
835 for (
const auto &cell :
mesh.active_cell_iterators())
837 for (
const auto v : cell->vertex_indices())
839 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
true;
840 for (
const auto vv : coinciding_vertex_groups
841 [vertex_to_coinciding_vertex_group[cell->vertex_index(v)]])
842 locally_active_vertices_on_subdomain[vv] =
true;
848 for (
const auto &cell :
mesh.active_cell_iterators())
849 if (!predicate(cell))
850 for (
const auto v : cell->vertex_indices())
851 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
854 active_halo_layer.push_back(cell);
858 return active_halo_layer;
863 template <
typename MeshType>
867 const MeshType &
mesh,
868 const std::function<
bool(
const typename MeshType::cell_iterator &)>
870 const unsigned int level)
872 std::vector<typename MeshType::cell_iterator> level_halo_layer;
873 std::vector<bool> locally_active_vertices_on_level_subdomain(
874 mesh.get_triangulation().n_vertices(),
false);
879 for (
typename MeshType::cell_iterator cell =
mesh.begin(
level);
883 for (
const unsigned int v : cell->vertex_indices())
884 locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
890 for (
typename MeshType::cell_iterator cell =
mesh.begin(
level);
893 if (!predicate(cell))
894 for (
const unsigned int v : cell->vertex_indices())
895 if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
898 level_halo_layer.push_back(cell);
902 return level_halo_layer;
908 template <
typename MeshType>
910 bool contains_locally_owned_cells(
911 const std::vector<typename MeshType::active_cell_iterator> &cells)
913 for (
typename std::vector<
914 typename MeshType::active_cell_iterator>::const_iterator
it =
919 if ((*it)->is_locally_owned())
925 template <
typename MeshType>
927 bool contains_artificial_cells(
928 const std::vector<typename MeshType::active_cell_iterator> &cells)
930 for (
typename std::vector<
931 typename MeshType::active_cell_iterator>::const_iterator
it =
936 if ((*it)->is_artificial())
945 template <
typename MeshType>
951 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
954 const std::vector<typename MeshType::active_cell_iterator>
959 Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) ==
false,
960 ExcMessage(
"Halo layer contains locally owned cells"));
961 Assert(contains_artificial_cells<MeshType>(active_halo_layer) ==
false,
962 ExcMessage(
"Halo layer contains artificial cells"));
964 return active_halo_layer;
969 template <
typename MeshType>
973 const MeshType &
mesh,
974 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
976 const double layer_thickness)
978 std::vector<typename MeshType::active_cell_iterator>
979 subdomain_boundary_cells, active_cell_layer_within_distance;
980 std::vector<bool> vertices_outside_subdomain(
981 mesh.get_triangulation().n_vertices(),
false);
983 const unsigned int spacedim = MeshType::space_dimension;
985 unsigned int n_non_predicate_cells = 0;
993 for (
const auto &cell :
mesh.active_cell_iterators())
994 if (!predicate(cell))
996 for (
const unsigned int v : cell->vertex_indices())
997 vertices_outside_subdomain[cell->vertex_index(v)] =
true;
998 ++n_non_predicate_cells;
1005 if (n_non_predicate_cells == 0 ||
1006 n_non_predicate_cells ==
mesh.get_triangulation().n_active_cells())
1007 return std::vector<typename MeshType::active_cell_iterator>();
1011 for (
const auto &cell :
mesh.active_cell_iterators())
1012 if (predicate(cell))
1014 for (
const unsigned int v : cell->vertex_indices())
1015 if (vertices_outside_subdomain[cell->vertex_index(v)] ==
true)
1017 subdomain_boundary_cells.push_back(cell);
1027 const double DOUBLE_EPSILON = 100. * std::numeric_limits<double>::epsilon();
1030 for (
unsigned int d = 0; d < spacedim; ++d)
1032 bounding_box.first[d] -= (layer_thickness + DOUBLE_EPSILON);
1033 bounding_box.second[d] += (layer_thickness + DOUBLE_EPSILON);
1036 std::vector<Point<spacedim>>
1037 subdomain_boundary_cells_centers;
1040 subdomain_boundary_cells_radii;
1042 subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size());
1043 subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size());
1045 for (
typename std::vector<typename MeshType::active_cell_iterator>::
1046 const_iterator subdomain_boundary_cell_iterator =
1047 subdomain_boundary_cells.begin();
1048 subdomain_boundary_cell_iterator != subdomain_boundary_cells.end();
1049 ++subdomain_boundary_cell_iterator)
1051 const std::pair<Point<spacedim>,
double>
1052 &subdomain_boundary_cell_enclosing_ball =
1053 (*subdomain_boundary_cell_iterator)->enclosing_ball();
1055 subdomain_boundary_cells_centers.push_back(
1056 subdomain_boundary_cell_enclosing_ball.first);
1057 subdomain_boundary_cells_radii.push_back(
1058 subdomain_boundary_cell_enclosing_ball.second);
1060 AssertThrow(subdomain_boundary_cells_radii.size() ==
1061 subdomain_boundary_cells_centers.size(),
1070 for (
const auto &cell :
mesh.active_cell_iterators())
1073 if (predicate(cell))
1076 const std::pair<Point<spacedim>,
double> &cell_enclosing_ball =
1077 cell->enclosing_ball();
1080 cell_enclosing_ball.first;
1081 const double cell_enclosing_ball_radius = cell_enclosing_ball.second;
1083 bool cell_inside =
true;
1085 for (
unsigned int d = 0; d < spacedim; ++d)
1087 (cell_enclosing_ball_center[d] + cell_enclosing_ball_radius >
1088 bounding_box.first[d]) &&
1089 (cell_enclosing_ball_center[d] - cell_enclosing_ball_radius <
1090 bounding_box.second[d]);
1096 for (
unsigned int i = 0; i < subdomain_boundary_cells_radii.size();
1099 subdomain_boundary_cells_centers[i]) <
1100 Utilities::fixed_power<2>(cell_enclosing_ball_radius +
1101 subdomain_boundary_cells_radii[i] +
1102 layer_thickness + DOUBLE_EPSILON))
1104 active_cell_layer_within_distance.push_back(cell);
1109 return active_cell_layer_within_distance;
1114 template <
typename MeshType>
1124 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1125 predicate(locally_owned_cell_predicate);
1127 const std::vector<typename MeshType::active_cell_iterator>
1128 ghost_cell_layer_within_distance =
1136 contains_locally_owned_cells<MeshType>(
1137 ghost_cell_layer_within_distance) ==
false,
1139 "Ghost cells within layer_thickness contains locally owned cells."));
1141 contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) ==
1144 "Ghost cells within layer_thickness contains artificial cells. "
1145 "The function compute_ghost_cell_layer_within_distance "
1146 "is probably called while using parallel::distributed::Triangulation. "
1147 "In such case please refer to the description of this function."));
1149 return ghost_cell_layer_within_distance;
1154 template <
typename MeshType>
1160 const std::function<
bool(
1161 const typename MeshType::
1162 active_cell_iterator &)>
1165 std::vector<bool> locally_active_vertices_on_subdomain(
1166 mesh.get_triangulation().n_vertices(),
false);
1168 const unsigned int spacedim = MeshType::space_dimension;
1175 for (
const auto &cell :
mesh.active_cell_iterators())
1176 if (predicate(cell))
1178 minp = cell->center();
1179 maxp = cell->center();
1185 for (
const auto &cell :
mesh.active_cell_iterators())
1186 if (predicate(cell))
1187 for (
const unsigned int v : cell->vertex_indices())
1188 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
1191 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
1193 for (
unsigned int d = 0; d < spacedim; ++d)
1195 minp[d] =
std::min(minp[d], cell->vertex(v)[d]);
1196 maxp[d] =
std::max(maxp[d], cell->vertex(v)[d]);
1200 return std::make_pair(minp, maxp);
1205 template <
typename MeshType>
1207 std::list<std::pair<
1208 typename MeshType::cell_iterator,
1215 ExcMessage(
"The two meshes must be represent triangulations that "
1216 "have the same coarse meshes"));
1223 bool remove_ghost_cells =
false;
1224#ifdef DEAL_II_WITH_MPI
1226 constexpr int dim = MeshType::dimension;
1227 constexpr int spacedim = MeshType::space_dimension;
1229 *
>(&mesh_1.get_triangulation()) !=
nullptr ||
1231 *
>(&mesh_2.get_triangulation()) !=
nullptr)
1233 Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(),
1234 ExcMessage(
"This function can only be used with meshes "
1235 "corresponding to distributed Triangulations when "
1236 "both Triangulations are equal."));
1237 remove_ghost_cells =
true;
1249 using CellList = std::list<std::pair<
typename MeshType::cell_iterator,
1250 typename MeshType::cell_iterator>>;
1254 typename MeshType::cell_iterator cell_1 = mesh_1.begin(0),
1255 cell_2 = mesh_2.begin(0);
1256 for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
1257 cell_list.emplace_back(cell_1, cell_2);
1260 typename CellList::iterator cell_pair = cell_list.begin();
1261 while (cell_pair != cell_list.end())
1265 if (cell_pair->first->has_children() &&
1266 cell_pair->second->has_children())
1268 Assert(cell_pair->first->refinement_case() ==
1269 cell_pair->second->refinement_case(),
1271 for (
unsigned int c = 0; c < cell_pair->first->n_children(); ++c)
1272 cell_list.emplace_back(cell_pair->first->child(c),
1273 cell_pair->second->child(c));
1278 const auto previous_cell_pair = cell_pair;
1280 cell_list.erase(previous_cell_pair);
1285 if (remove_ghost_cells &&
1286 ((cell_pair->first->is_active() &&
1287 !cell_pair->first->is_locally_owned()) ||
1288 (cell_pair->second->is_active() &&
1289 !cell_pair->second->is_locally_owned())))
1292 const auto previous_cell_pair = cell_pair;
1294 cell_list.erase(previous_cell_pair);
1303 for (cell_pair = cell_list.begin(); cell_pair != cell_list.end();
1305 Assert(cell_pair->first->is_active() || cell_pair->second->is_active() ||
1306 (cell_pair->first->refinement_case() !=
1307 cell_pair->second->refinement_case()),
1315 template <
int dim,
int spacedim>
1332 endc = mesh_1.
end(0);
1333 for (; cell_1 != endc; ++cell_1, ++cell_2)
1335 if (cell_1->n_vertices() != cell_2->n_vertices())
1337 for (
const unsigned int v : cell_1->vertex_indices())
1338 if (cell_1->vertex(v) != cell_2->vertex(v))
1351 template <
typename MeshType>
1356 mesh_2.get_triangulation());
1361 template <
int dim,
int spacedim>
1362 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1368 const double tolerance)
1371 (
mapping.size() ==
mesh.get_fe_collection().size()),
1372 ExcMessage(
"Mapping collection needs to have either size 1 "
1373 "or size equal to the number of elements in "
1374 "the FECollection."));
1376 using cell_iterator =
1379 std::pair<cell_iterator, Point<dim>> best_cell;
1394 double best_distance = tolerance;
1395 int best_level = -1;
1402 std::vector<cell_iterator> adjacent_cells_tmp =
1410 std::set<cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
1411 adjacent_cells_tmp.end());
1412 std::set<cell_iterator> searched_cells;
1420 const auto n_cells =
mesh.get_triangulation().n_cells();
1422 unsigned int cells_searched = 0;
1423 while (!found && cells_searched < n_cells)
1430 mapping[cell->active_fe_index()]
1431 .transform_real_to_unit_cell(cell,
p);
1437 cell->reference_cell().closest_point(p_cell).
distance(
1444 if (dist < best_distance ||
1445 (dist == best_distance && cell->level() > best_level))
1448 best_distance = dist;
1449 best_level = cell->level();
1450 best_cell = std::make_pair(cell, p_cell);
1476 if (!found && cells_searched < n_cells)
1478 find_active_cell_around_point_internal<dim,
1490 template <
typename MeshType>
1493 const typename MeshType::active_cell_iterator &cell)
1495 Assert(cell->is_locally_owned(),
1496 ExcMessage(
"This function only makes sense if the cell for "
1497 "which you are asking for a patch, is locally "
1500 std::vector<typename MeshType::active_cell_iterator> patch;
1501 patch.push_back(cell);
1502 for (
const unsigned int face_number : cell->face_indices())
1503 if (cell->face(face_number)->at_boundary() ==
false)
1505 if (cell->neighbor(face_number)->has_children() ==
false)
1506 patch.push_back(cell->neighbor(face_number));
1511 if (MeshType::dimension > 1)
1513 for (
unsigned int subface = 0;
1514 subface < cell->face(face_number)->n_children();
1517 cell->neighbor_child_on_subface(face_number, subface));
1523 typename MeshType::cell_iterator neighbor =
1524 cell->neighbor(face_number);
1525 while (neighbor->has_children())
1526 neighbor = neighbor->child(1 - face_number);
1528 Assert(neighbor->neighbor(1 - face_number) == cell,
1530 patch.push_back(neighbor);
1538 template <
class Container>
1539 std::vector<typename Container::cell_iterator>
1541 const std::vector<typename Container::active_cell_iterator> &patch)
1545 "Vector containing patch cells should not be an empty vector!"));
1549 int min_level = patch[0]->level();
1551 for (
unsigned int i = 0; i < patch.size(); ++i)
1553 std::set<typename Container::cell_iterator> uniform_cells;
1554 typename std::vector<
1555 typename Container::active_cell_iterator>::const_iterator patch_cell;
1557 for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell)
1562 if ((*patch_cell)->level() == min_level)
1563 uniform_cells.insert(*patch_cell);
1570 typename Container::cell_iterator parent = *patch_cell;
1572 while (parent->level() > min_level)
1573 parent = parent->parent();
1574 uniform_cells.insert(parent);
1578 return std::vector<typename Container::cell_iterator>(uniform_cells.begin(),
1579 uniform_cells.end());
1584 template <
class Container>
1587 const std::vector<typename Container::active_cell_iterator> &patch,
1589 &local_triangulation,
1592 Container::space_dimension>::active_cell_iterator,
1593 typename Container::active_cell_iterator> &patch_to_global_tria_map)
1596 const std::vector<typename Container::cell_iterator> uniform_cells =
1597 get_cells_at_coarsest_common_level<Container>(patch);
1599 local_triangulation.
clear();
1600 std::vector<Point<Container::space_dimension>>
vertices;
1601 const unsigned int n_uniform_cells = uniform_cells.size();
1602 std::vector<CellData<Container::dimension>> cells(n_uniform_cells);
1605 typename std::vector<typename Container::cell_iterator>::const_iterator
1607 for (uniform_cell = uniform_cells.begin();
1608 uniform_cell != uniform_cells.end();
1611 for (
const unsigned int v : (*uniform_cell)->vertex_indices())
1614 (*uniform_cell)->vertex(v);
1615 bool repeat_vertex =
false;
1617 for (
unsigned int m = 0; m < i; ++m)
1621 repeat_vertex =
true;
1622 cells[k].vertices[v] = m;
1626 if (repeat_vertex ==
false)
1629 cells[k].vertices[v] = i;
1640 unsigned int index = 0;
1643 Container::space_dimension>::cell_iterator,
1644 typename Container::cell_iterator>
1645 patch_to_global_tria_map_tmp;
1647 Container::space_dimension>::cell_iterator
1648 coarse_cell = local_triangulation.
begin();
1649 coarse_cell != local_triangulation.
end();
1650 ++coarse_cell, ++index)
1652 patch_to_global_tria_map_tmp.insert(
1653 std::make_pair(coarse_cell, uniform_cells[index]));
1657 Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <=
1658 1e-15 * coarse_cell->diameter(),
1661 bool refinement_necessary;
1666 refinement_necessary =
false;
1667 for (
const auto &active_tria_cell :
1670 if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
1672 active_tria_cell->set_refine_flag();
1673 refinement_necessary =
true;
1676 for (
unsigned int i = 0; i < patch.size(); ++i)
1681 if (patch_to_global_tria_map_tmp[active_tria_cell] ==
1686 for (
const unsigned int v :
1687 active_tria_cell->vertex_indices())
1688 active_tria_cell->vertex(v) = patch[i]->vertex(v);
1690 Assert(active_tria_cell->center().distance(
1691 patch_to_global_tria_map_tmp[active_tria_cell]
1693 1e-15 * active_tria_cell->diameter(),
1696 active_tria_cell->set_user_flag();
1702 if (refinement_necessary)
1707 Container::dimension,
1708 Container::space_dimension>::cell_iterator cell =
1709 local_triangulation.
begin();
1710 cell != local_triangulation.
end();
1713 if (patch_to_global_tria_map_tmp.find(cell) !=
1714 patch_to_global_tria_map_tmp.end())
1716 if (cell->has_children())
1723 for (
unsigned int c = 0; c < cell->n_children(); ++c)
1725 if (patch_to_global_tria_map_tmp.find(cell->child(
1726 c)) == patch_to_global_tria_map_tmp.end())
1728 patch_to_global_tria_map_tmp.insert(
1731 patch_to_global_tria_map_tmp[cell]->child(
1753 patch_to_global_tria_map_tmp.erase(cell);
1759 while (refinement_necessary);
1765 Container::space_dimension>::cell_iterator
1766 cell = local_triangulation.
begin();
1767 cell != local_triangulation.
end();
1770 if (cell->user_flag_set())
1772 Assert(patch_to_global_tria_map_tmp.find(cell) !=
1773 patch_to_global_tria_map_tmp.end(),
1776 Assert(cell->center().distance(
1777 patch_to_global_tria_map_tmp[cell]->center()) <=
1778 1e-15 * cell->diameter(),
1786 Container::space_dimension>::cell_iterator,
1787 typename Container::cell_iterator>::iterator
1788 map_tmp_it = patch_to_global_tria_map_tmp.
begin(),
1789 map_tmp_end = patch_to_global_tria_map_tmp.end();
1793 for (; map_tmp_it != map_tmp_end; ++map_tmp_it)
1794 patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
1799 template <
int dim,
int spacedim>
1802 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1816 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1817 dof_to_set_of_cells_map;
1819 std::vector<types::global_dof_index> local_dof_indices;
1820 std::vector<types::global_dof_index> local_face_dof_indices;
1821 std::vector<types::global_dof_index> local_line_dof_indices;
1825 std::vector<bool> user_flags;
1830 std::map<typename DoFHandler<dim, spacedim>::active_line_iterator,
1832 lines_to_parent_lines_map;
1839 .clear_user_flags();
1844 endc = dof_handler.
end();
1845 for (; cell != endc; ++cell)
1851 if (cell->is_artificial() ==
false)
1853 for (
unsigned int l = 0; l < cell->n_lines(); ++l)
1854 if (cell->line(l)->has_children())
1855 for (
unsigned int c = 0; c < cell->line(l)->n_children();
1858 lines_to_parent_lines_map[cell->line(l)->child(c)] =
1862 cell->line(l)->child(c)->set_user_flag();
1877 endc = dof_handler.
end();
1878 for (; cell != endc; ++cell)
1883 if (cell->is_artificial() ==
false)
1885 const unsigned int n_dofs_per_cell =
1887 local_dof_indices.resize(n_dofs_per_cell);
1891 cell->get_dof_indices(local_dof_indices);
1892 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i)
1893 dof_to_set_of_cells_map[local_dof_indices[i]].insert(cell);
1903 for (
const unsigned int f : cell->face_indices())
1905 if (cell->face(f)->has_children())
1907 for (
unsigned int c = 0; c < cell->face(f)->n_children();
1922 Assert(cell->face(f)->child(c)->has_children() ==
false,
1925 const unsigned int n_dofs_per_face =
1927 local_face_dof_indices.resize(n_dofs_per_face);
1929 cell->face(f)->child(c)->get_dof_indices(
1930 local_face_dof_indices);
1931 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1932 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1936 else if ((cell->face(f)->at_boundary() ==
false) &&
1937 (cell->neighbor_is_coarser(f)))
1954 std::pair<unsigned int, unsigned int>
1955 neighbor_face_no_subface_no =
1956 cell->neighbor_of_coarser_neighbor(f);
1957 unsigned int face_no = neighbor_face_no_subface_no.first;
1958 unsigned int subface = neighbor_face_no_subface_no.second;
1960 const unsigned int n_dofs_per_face =
1962 local_face_dof_indices.resize(n_dofs_per_face);
1964 cell->neighbor(f)->face(face_no)->get_dof_indices(
1965 local_face_dof_indices);
1966 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1967 dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(
1972 for (
unsigned int c = 0;
1973 c < cell->neighbor(f)->face(face_no)->n_children();
1979 const unsigned int n_dofs_per_face =
1981 local_face_dof_indices.resize(n_dofs_per_face);
1986 ->has_children() ==
false,
1991 ->get_dof_indices(local_face_dof_indices);
1992 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1993 dof_to_set_of_cells_map[local_face_dof_indices[i]]
2010 for (
unsigned int l = 0; l < cell->n_lines(); ++l)
2012 if (cell->line(l)->has_children())
2014 for (
unsigned int c = 0;
2015 c < cell->line(l)->n_children();
2018 Assert(cell->line(l)->child(c)->has_children() ==
2024 const unsigned int n_dofs_per_line =
2027 local_line_dof_indices.resize(n_dofs_per_line);
2029 cell->line(l)->child(c)->get_dof_indices(
2030 local_line_dof_indices);
2031 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2032 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2040 else if (cell->line(l)->user_flag_set() ==
true)
2044 lines_to_parent_lines_map[cell->line(l)];
2049 const unsigned int n_dofs_per_line =
2052 local_line_dof_indices.resize(n_dofs_per_line);
2054 parent_line->get_dof_indices(local_line_dof_indices);
2055 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2056 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2059 for (
unsigned int c = 0; c < parent_line->n_children();
2062 Assert(parent_line->child(c)->has_children() ==
2066 const unsigned int n_dofs_per_line =
2069 local_line_dof_indices.resize(n_dofs_per_line);
2071 parent_line->child(c)->get_dof_indices(
2072 local_line_dof_indices);
2073 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2074 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2091 .load_user_flags(user_flags);
2098 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2099 dof_to_cell_patches;
2103 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>::
2104 iterator
it = dof_to_set_of_cells_map.begin(),
2105 it_end = dof_to_set_of_cells_map.end();
2106 for (;
it != it_end; ++
it)
2107 dof_to_cell_patches[
it->first].assign(
it->second.begin(),
2110 return dof_to_cell_patches;
2116 template <
typename CellIterator>
2119 std::set<std::pair<CellIterator, unsigned int>> &pairs1,
2122 const unsigned int direction,
2124 const ::Tensor<1, CellIterator::AccessorType::space_dimension>
2128 static const int space_dim = CellIterator::AccessorType::space_dimension;
2134 constexpr int dim = CellIterator::AccessorType::dimension;
2135 constexpr int spacedim = CellIterator::AccessorType::space_dimension;
2140 if (!(((pairs1.size() > 0) &&
2141 (
dynamic_cast<const parallel::fullydistributed::
2142 Triangulation<dim, spacedim> *
>(
2143 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2144 ((pairs2.size() > 0) &&
2147 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2148 Assert(pairs1.size() == pairs2.size(),
2149 ExcMessage(
"Unmatched faces on periodic boundaries"));
2153 unsigned int n_matches = 0;
2156 using PairIterator =
2157 typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator;
2158 for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
2160 for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
2162 const CellIterator cell1 = it1->first;
2163 const CellIterator cell2 = it2->first;
2164 const unsigned int face_idx1 = it1->second;
2165 const unsigned int face_idx2 = it2->second;
2166 if (
const std::optional<unsigned char> orientation =
2168 cell2->face(face_idx2),
2178 {face_idx1, face_idx2},
2179 orientation.value(),
2181 matched_pairs.push_back(matched_face);
2195 constexpr int dim = CellIterator::AccessorType::dimension;
2196 constexpr int spacedim = CellIterator::AccessorType::space_dimension;
2197 if (!(((pairs1.size() > 0) &&
2198 (
dynamic_cast<const parallel::fullydistributed::
2199 Triangulation<dim, spacedim> *
>(
2200 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2201 ((pairs2.size() > 0) &&
2204 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2205 AssertThrow(n_matches == pairs1.size() && pairs2.empty(),
2206 ExcMessage(
"Unmatched faces on periodic boundaries"));
2212 template <
typename MeshType>
2215 const MeshType &
mesh,
2216 const
types::boundary_id b_id,
2217 const
unsigned int direction,
2220 const
Tensor<1, MeshType::space_dimension> &offset,
2223 static const int dim = MeshType::dimension;
2224 static const int space_dim = MeshType::space_dimension;
2234 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2235 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2237 for (
typename MeshType::cell_iterator cell =
mesh.begin(0);
2238 cell !=
mesh.end(0);
2241 const typename MeshType::face_iterator face_1 =
2242 cell->face(2 * direction);
2243 const typename MeshType::face_iterator face_2 =
2244 cell->face(2 * direction + 1);
2246 if (face_1->at_boundary() && face_1->boundary_id() == b_id)
2248 const std::pair<typename MeshType::cell_iterator, unsigned int>
2249 pair1 = std::make_pair(cell, 2 * direction);
2250 pairs1.insert(pair1);
2253 if (face_2->at_boundary() && face_2->boundary_id() == b_id)
2255 const std::pair<typename MeshType::cell_iterator, unsigned int>
2256 pair2 = std::make_pair(cell, 2 * direction + 1);
2257 pairs2.insert(pair2);
2261 Assert(pairs1.size() == pairs2.size(),
2262 ExcMessage(
"Unmatched faces on periodic boundaries"));
2264 Assert(pairs1.size() > 0,
2265 ExcMessage(
"No new periodic face pairs have been found. "
2266 "Are you sure that you've selected the correct boundary "
2267 "id's and that the coarsest level mesh is colorized?"));
2270 const unsigned int size_old = matched_pairs.size();
2275 pairs1, pairs2, direction, matched_pairs, offset, matrix);
2279 const unsigned int size_new = matched_pairs.size();
2280 for (
unsigned int i = size_old; i < size_new; ++i)
2282 Assert(matched_pairs[i].orientation ==
2285 "Found a face match with non standard orientation. "
2286 "This function is only suitable for meshes with cells "
2287 "in default orientation"));
2294 template <
typename MeshType>
2297 const MeshType &
mesh,
2298 const
types::boundary_id b_id1,
2299 const
types::boundary_id b_id2,
2300 const
unsigned int direction,
2303 const
Tensor<1, MeshType::space_dimension> &offset,
2306 static const int dim = MeshType::dimension;
2307 static const int space_dim = MeshType::space_dimension;
2315 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2316 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2318 for (
typename MeshType::cell_iterator cell =
mesh.begin(0);
2319 cell !=
mesh.end(0);
2322 for (
const unsigned int i : cell->face_indices())
2324 const typename MeshType::face_iterator face = cell->face(i);
2325 if (face->at_boundary() && face->boundary_id() == b_id1)
2327 const std::pair<typename MeshType::cell_iterator, unsigned int>
2328 pair1 = std::make_pair(cell, i);
2329 pairs1.insert(pair1);
2332 if (face->at_boundary() && face->boundary_id() == b_id2)
2334 const std::pair<typename MeshType::cell_iterator, unsigned int>
2335 pair2 = std::make_pair(cell, i);
2336 pairs2.insert(pair2);
2347 if (!(((pairs1.size() > 0) &&
2350 *
>(&pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2351 ((pairs2.size() > 0) &&
2354 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2355 Assert(pairs1.size() == pairs2.size(),
2356 ExcMessage(
"Unmatched faces on periodic boundaries"));
2359 (pairs1.size() > 0 ||
2362 &
mesh.begin()->get_triangulation()) !=
nullptr)),
2363 ExcMessage(
"No new periodic face pairs have been found. "
2364 "Are you sure that you've selected the correct boundary "
2365 "id's and that the coarsest level mesh is colorized?"));
2369 pairs1, pairs2, direction, matched_pairs, offset, matrix);
2383 template <
int spacedim>
2387 const unsigned int direction,
2397 if (matrix.m() == spacedim)
2398 for (
unsigned int i = 0; i < spacedim; ++i)
2399 for (
unsigned int j = 0; j < spacedim; ++j)
2400 distance[i] += matrix(i, j) * point1[j];
2404 distance += offset - point2;
2406 for (
unsigned int i = 0; i < spacedim; ++i)
2412 if (
std::abs(distance[i]) > 1.e-10)
2421 template <
typename FaceIterator>
2422 std::optional<unsigned char>
2424 const FaceIterator &face1,
2425 const FaceIterator &face2,
2426 const unsigned int direction,
2430 Assert(matrix.m() == matrix.n(),
2431 ExcMessage(
"The supplied matrix must be a square matrix"));
2433 static const int dim = FaceIterator::AccessorType::dimension;
2437 std::array<unsigned int, GeometryInfo<dim>::vertices_per_face>
2438 face1_vertices, face2_vertices;
2442 std::set<unsigned int> face2_vertices_set;
2443 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2444 face2_vertices_set.insert(i);
2446 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2448 for (
auto it = face2_vertices_set.begin();
2449 it != face2_vertices_set.end();
2458 face1_vertices[i] = *
it;
2459 face2_vertices[i] = i;
2460 face2_vertices_set.erase(
it);
2466 if (face2_vertices_set.empty())
2468 const auto reference_cell = face1->reference_cell();
2471 return std::make_optional(reference_cell.get_combined_orientation(
2473 face2_vertices.cbegin() + face2->n_vertices()),
2475 face1_vertices.cbegin() + face1->n_vertices())));
2478 return std::nullopt;
2483#include "grid_tools_dof_handlers.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
cell_iterator end() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const Triangulation< dim, spacedim > & get_triangulation() const
active_cell_iterator begin_active(const unsigned int level=0) 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
Abstract base class for mapping classes.
constexpr numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
static constexpr unsigned char default_combined_face_orientation()
cell_iterator begin(const unsigned int level=0) const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
unsigned int n_active_cells() const
void save_user_flags(std::ostream &out) const
const std::vector< Point< spacedim > > & get_vertices() const
cell_iterator end() const
virtual void execute_coarsening_and_refinement()
unsigned int n_cells() const
const std::vector< bool > & get_used_vertices() const
Triangulation< dim, spacedim > & get_triangulation()
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
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::active_cell_iterator active_cell_iterator
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static const unsigned int invalid_unsigned_int
typename type_identity< T >::type type_identity_t
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index