54 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
58 const
Point<spacedim> &p,
59 const std::vector<
bool> &marked_vertices)
70 marked_vertices.empty(),
72 marked_vertices.size()));
82 marked_vertices.empty() ||
83 std::equal(marked_vertices.begin(),
84 marked_vertices.end(),
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);
109 unsigned int best_vertex = std::distance(used.begin(),
first);
110 double best_dist = (p -
vertices[best_vertex]).norm_square();
114 for (
unsigned int j = best_vertex + 1; j <
vertices.size(); ++j)
117 double dist = (p -
vertices[j]).norm_square();
118 if (dist < best_dist)
130 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
134 const MeshType<dim, spacedim> &mesh,
135 const
Point<spacedim> &p,
136 const std::vector<
bool> &marked_vertices)
139 if (mapping.preserves_vertex_locations() ==
true)
151 marked_vertices.empty(),
153 marked_vertices.size()));
163 marked_vertices.empty() ||
164 std::equal(marked_vertices.begin(),
165 marked_vertices.end(),
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!"));
173 if (marked_vertices.size())
176 if (marked_vertices[it->first] ==
false)
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],
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>
351 typename ::internal::
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();
392 for (; it != mesh.end(); ++it)
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
415 const std::vector<bool> &marked_vertices,
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,
442 const std::vector<bool> &marked_vertices,
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;
487 if (cell->is_artificial() ==
false)
490 if (marked_vertices.size() > 0)
492 bool any_vertex_marked =
false;
493 for (
const auto &v : cell->vertex_indices())
495 if (marked_vertices[cell->vertex_index(v)])
497 any_vertex_marked =
true;
501 if (!any_vertex_marked)
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);
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,
586 const std::vector<bool> &marked_vertices)
589 mapping, mesh, p, marked_vertices, tolerance);
591 if (cell_and_point.first == mesh.end())
595 mapping, mesh, p, tolerance, cell_and_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,
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]);
761 for (
const auto v : my_cell->vertex_indices())
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)
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,
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);
829 for (
const auto &cell : mesh.active_cell_iterators())
831 for (
const auto v : cell->vertex_indices())
832 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
true;
837 for (
const auto &cell : mesh.active_cell_iterators())
838 if (!predicate(cell))
839 for (
const auto v : cell->vertex_indices())
840 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
843 active_halo_layer.push_back(cell);
847 return active_halo_layer;
852 template <
typename MeshType>
856 const MeshType &mesh,
857 const std::function<
bool(
const typename MeshType::cell_iterator &)>
859 const unsigned int level)
861 std::vector<typename MeshType::cell_iterator> level_halo_layer;
862 std::vector<bool> locally_active_vertices_on_level_subdomain(
863 mesh.get_triangulation().n_vertices(),
false);
868 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
869 cell != mesh.end(
level);
872 for (
const unsigned int v : cell->vertex_indices())
873 locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
879 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
880 cell != mesh.end(
level);
882 if (!predicate(cell))
883 for (
const unsigned int v : cell->vertex_indices())
884 if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
887 level_halo_layer.push_back(cell);
891 return level_halo_layer;
897 template <
typename MeshType>
899 bool contains_locally_owned_cells(
900 const std::vector<typename MeshType::active_cell_iterator> &cells)
902 for (
typename std::vector<
903 typename MeshType::active_cell_iterator>::const_iterator it =
908 if ((*it)->is_locally_owned())
914 template <
typename MeshType>
916 bool contains_artificial_cells(
917 const std::vector<typename MeshType::active_cell_iterator> &cells)
919 for (
typename std::vector<
920 typename MeshType::active_cell_iterator>::const_iterator it =
925 if ((*it)->is_artificial())
934 template <
typename MeshType>
940 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
943 const std::vector<typename MeshType::active_cell_iterator>
948 Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) ==
false,
949 ExcMessage(
"Halo layer contains locally owned cells"));
950 Assert(contains_artificial_cells<MeshType>(active_halo_layer) ==
false,
951 ExcMessage(
"Halo layer contains artificial cells"));
953 return active_halo_layer;
958 template <
typename MeshType>
962 const MeshType &mesh,
963 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
965 const double layer_thickness)
967 std::vector<typename MeshType::active_cell_iterator>
968 subdomain_boundary_cells, active_cell_layer_within_distance;
969 std::vector<bool> vertices_outside_subdomain(
970 mesh.get_triangulation().n_vertices(),
false);
972 const unsigned int spacedim = MeshType::space_dimension;
974 unsigned int n_non_predicate_cells = 0;
982 for (
const auto &cell : mesh.active_cell_iterators())
983 if (!predicate(cell))
985 for (
const unsigned int v : cell->vertex_indices())
986 vertices_outside_subdomain[cell->vertex_index(v)] =
true;
987 n_non_predicate_cells++;
994 if (n_non_predicate_cells == 0 ||
995 n_non_predicate_cells == mesh.get_triangulation().n_active_cells())
996 return std::vector<typename MeshType::active_cell_iterator>();
1000 for (
const auto &cell : mesh.active_cell_iterators())
1001 if (predicate(cell))
1003 for (
const unsigned int v : cell->vertex_indices())
1004 if (vertices_outside_subdomain[cell->vertex_index(v)] ==
true)
1006 subdomain_boundary_cells.push_back(cell);
1019 for (
unsigned int d = 0;
d < spacedim; ++
d)
1021 bounding_box.first[
d] -= (layer_thickness + DOUBLE_EPSILON);
1022 bounding_box.second[
d] += (layer_thickness + DOUBLE_EPSILON);
1025 std::vector<Point<spacedim>>
1026 subdomain_boundary_cells_centers;
1029 subdomain_boundary_cells_radii;
1031 subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size());
1032 subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size());
1034 for (
typename std::vector<typename MeshType::active_cell_iterator>::
1035 const_iterator subdomain_boundary_cell_iterator =
1036 subdomain_boundary_cells.begin();
1037 subdomain_boundary_cell_iterator != subdomain_boundary_cells.end();
1038 ++subdomain_boundary_cell_iterator)
1040 const std::pair<Point<spacedim>,
double>
1041 &subdomain_boundary_cell_enclosing_ball =
1042 (*subdomain_boundary_cell_iterator)->enclosing_ball();
1044 subdomain_boundary_cells_centers.push_back(
1045 subdomain_boundary_cell_enclosing_ball.first);
1046 subdomain_boundary_cells_radii.push_back(
1047 subdomain_boundary_cell_enclosing_ball.second);
1049 AssertThrow(subdomain_boundary_cells_radii.size() ==
1050 subdomain_boundary_cells_centers.size(),
1059 for (
const auto &cell : mesh.active_cell_iterators())
1062 if (predicate(cell))
1065 const std::pair<Point<spacedim>,
double> &cell_enclosing_ball =
1066 cell->enclosing_ball();
1069 cell_enclosing_ball.first;
1070 const double cell_enclosing_ball_radius = cell_enclosing_ball.second;
1072 bool cell_inside =
true;
1074 for (
unsigned int d = 0;
d < spacedim; ++
d)
1076 (cell_enclosing_ball_center[
d] + cell_enclosing_ball_radius >
1077 bounding_box.first[
d]) &&
1078 (cell_enclosing_ball_center[
d] - cell_enclosing_ball_radius <
1079 bounding_box.second[
d]);
1085 for (
unsigned int i = 0; i < subdomain_boundary_cells_radii.size();
1088 subdomain_boundary_cells_centers[i]) <
1089 Utilities::fixed_power<2>(cell_enclosing_ball_radius +
1090 subdomain_boundary_cells_radii[i] +
1091 layer_thickness + DOUBLE_EPSILON))
1093 active_cell_layer_within_distance.push_back(cell);
1098 return active_cell_layer_within_distance;
1103 template <
typename MeshType>
1113 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1114 predicate(locally_owned_cell_predicate);
1116 const std::vector<typename MeshType::active_cell_iterator>
1117 ghost_cell_layer_within_distance =
1125 contains_locally_owned_cells<MeshType>(
1126 ghost_cell_layer_within_distance) ==
false,
1128 "Ghost cells within layer_thickness contains locally owned cells."));
1130 contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) ==
1133 "Ghost cells within layer_thickness contains artificial cells. "
1134 "The function compute_ghost_cell_layer_within_distance "
1135 "is probably called while using parallel::distributed::Triangulation. "
1136 "In such case please refer to the description of this function."));
1138 return ghost_cell_layer_within_distance;
1143 template <
typename MeshType>
1149 const std::function<
bool(
1150 const typename MeshType::
1151 active_cell_iterator &)>
1154 std::vector<bool> locally_active_vertices_on_subdomain(
1155 mesh.get_triangulation().n_vertices(),
false);
1157 const unsigned int spacedim = MeshType::space_dimension;
1164 for (
const auto &cell : mesh.active_cell_iterators())
1165 if (predicate(cell))
1167 minp = cell->center();
1168 maxp = cell->center();
1174 for (
const auto &cell : mesh.active_cell_iterators())
1175 if (predicate(cell))
1176 for (
const unsigned int v : cell->vertex_indices())
1177 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
1180 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
1182 for (
unsigned int d = 0;
d < spacedim; ++
d)
1189 return std::make_pair(minp, maxp);
1194 template <
typename MeshType>
1196 std::list<std::pair<
1197 typename MeshType::cell_iterator,
1204 ExcMessage(
"The two meshes must be represent triangulations that "
1205 "have the same coarse meshes"));
1212 bool remove_ghost_cells =
false;
1213 #ifdef DEAL_II_WITH_MPI
1215 constexpr
int dim = MeshType::dimension;
1216 constexpr
int spacedim = MeshType::space_dimension;
1218 *
>(&mesh_1.get_triangulation()) !=
nullptr ||
1220 *
>(&mesh_2.get_triangulation()) !=
nullptr)
1222 Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(),
1223 ExcMessage(
"This function can only be used with meshes "
1224 "corresponding to distributed Triangulations when "
1225 "both Triangulations are equal."));
1226 remove_ghost_cells =
true;
1238 using CellList = std::list<std::pair<
typename MeshType::cell_iterator,
1239 typename MeshType::cell_iterator>>;
1243 typename MeshType::cell_iterator cell_1 = mesh_1.begin(0),
1244 cell_2 = mesh_2.begin(0);
1245 for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
1246 cell_list.emplace_back(cell_1, cell_2);
1249 typename CellList::iterator cell_pair = cell_list.begin();
1250 while (cell_pair != cell_list.end())
1254 if (cell_pair->first->has_children() &&
1255 cell_pair->second->has_children())
1257 Assert(cell_pair->first->refinement_case() ==
1258 cell_pair->second->refinement_case(),
1260 for (
unsigned int c = 0; c < cell_pair->first->n_children(); ++c)
1261 cell_list.emplace_back(cell_pair->first->child(c),
1262 cell_pair->second->child(c));
1267 const auto previous_cell_pair = cell_pair;
1269 cell_list.erase(previous_cell_pair);
1274 if (remove_ghost_cells &&
1275 ((cell_pair->first->is_active() &&
1276 !cell_pair->first->is_locally_owned()) ||
1277 (cell_pair->second->is_active() &&
1278 !cell_pair->second->is_locally_owned())))
1281 const auto previous_cell_pair = cell_pair;
1283 cell_list.erase(previous_cell_pair);
1292 for (cell_pair = cell_list.begin(); cell_pair != cell_list.end();
1294 Assert(cell_pair->first->is_active() || cell_pair->second->is_active() ||
1295 (cell_pair->first->refinement_case() !=
1296 cell_pair->second->refinement_case()),
1304 template <
int dim,
int spacedim>
1321 endc = mesh_1.
end(0);
1322 for (; cell_1 != endc; ++cell_1, ++cell_2)
1324 if (cell_1->
n_vertices() != cell_2->n_vertices())
1326 for (
const unsigned int v : cell_1->vertex_indices())
1327 if (cell_1->vertex(v) != cell_2->vertex(v))
1340 template <
typename MeshType>
1345 mesh_2.get_triangulation());
1350 template <
int dim,
int spacedim>
1351 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1357 const double tolerance)
1361 ExcMessage(
"Mapping collection needs to have either size 1 "
1362 "or size equal to the number of elements in "
1363 "the FECollection."));
1365 using cell_iterator =
1368 std::pair<cell_iterator, Point<dim>> best_cell;
1372 if (mapping.
size() == 1)
1374 const std::vector<bool> marked_vertices = {};
1376 mapping[0], mesh, p, marked_vertices, tolerance);
1383 double best_distance = tolerance;
1384 int best_level = -1;
1391 std::vector<cell_iterator> adjacent_cells_tmp =
1399 std::set<cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
1400 adjacent_cells_tmp.end());
1401 std::set<cell_iterator> searched_cells;
1411 unsigned int cells_searched = 0;
1412 while (!found && cells_searched <
n_cells)
1419 mapping[cell->active_fe_index()]
1420 .transform_real_to_unit_cell(cell, p);
1426 cell->reference_cell().closest_point(p_cell).distance(
1433 if (dist < best_distance ||
1434 (dist == best_distance && cell->level() > best_level))
1437 best_distance = dist;
1438 best_level = cell->level();
1439 best_cell = std::make_pair(cell, p_cell);
1465 if (!found && cells_searched <
n_cells)
1467 find_active_cell_around_point_internal<dim,
1479 template <
typename MeshType>
1482 const typename MeshType::active_cell_iterator &cell)
1484 Assert(cell->is_locally_owned(),
1485 ExcMessage(
"This function only makes sense if the cell for "
1486 "which you are asking for a patch, is locally "
1489 std::vector<typename MeshType::active_cell_iterator> patch;
1490 patch.push_back(cell);
1491 for (
const unsigned int face_number : cell->face_indices())
1492 if (cell->face(face_number)->at_boundary() ==
false)
1494 if (cell->neighbor(face_number)->has_children() ==
false)
1495 patch.push_back(cell->neighbor(face_number));
1500 if (MeshType::dimension > 1)
1502 for (
unsigned int subface = 0;
1503 subface < cell->face(face_number)->n_children();
1506 cell->neighbor_child_on_subface(face_number, subface));
1512 typename MeshType::cell_iterator neighbor =
1513 cell->neighbor(face_number);
1514 while (neighbor->has_children())
1515 neighbor = neighbor->child(1 - face_number);
1517 Assert(neighbor->neighbor(1 - face_number) == cell,
1519 patch.push_back(neighbor);
1527 template <
class Container>
1528 std::vector<typename Container::cell_iterator>
1530 const std::vector<typename Container::active_cell_iterator> &patch)
1534 "Vector containing patch cells should not be an empty vector!"));
1538 int min_level = patch[0]->level();
1540 for (
unsigned int i = 0; i < patch.size(); ++i)
1542 std::set<typename Container::cell_iterator> uniform_cells;
1543 typename std::vector<
1544 typename Container::active_cell_iterator>::const_iterator patch_cell;
1546 for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell)
1551 if ((*patch_cell)->level() == min_level)
1552 uniform_cells.insert(*patch_cell);
1559 typename Container::cell_iterator parent = *patch_cell;
1561 while (parent->level() > min_level)
1562 parent = parent->parent();
1563 uniform_cells.insert(parent);
1567 return std::vector<typename Container::cell_iterator>(uniform_cells.begin(),
1568 uniform_cells.end());
1573 template <
class Container>
1576 const std::vector<typename Container::active_cell_iterator> &patch,
1578 &local_triangulation,
1581 Container::space_dimension>::active_cell_iterator,
1582 typename Container::active_cell_iterator> &patch_to_global_tria_map)
1585 const std::vector<typename Container::cell_iterator> uniform_cells =
1586 get_cells_at_coarsest_common_level<Container>(patch);
1588 local_triangulation.
clear();
1589 std::vector<Point<Container::space_dimension>>
vertices;
1590 const unsigned int n_uniform_cells = uniform_cells.size();
1591 std::vector<CellData<Container::dimension>> cells(n_uniform_cells);
1594 typename std::vector<typename Container::cell_iterator>::const_iterator
1596 for (uniform_cell = uniform_cells.begin();
1597 uniform_cell != uniform_cells.end();
1600 for (
const unsigned int v : (*uniform_cell)->vertex_indices())
1603 (*uniform_cell)->vertex(v);
1604 bool repeat_vertex =
false;
1606 for (
unsigned int m = 0; m < i; ++m)
1610 repeat_vertex =
true;
1611 cells[k].vertices[v] = m;
1615 if (repeat_vertex ==
false)
1618 cells[k].vertices[v] = i;
1629 unsigned int index = 0;
1632 Container::space_dimension>::cell_iterator,
1633 typename Container::cell_iterator>
1634 patch_to_global_tria_map_tmp;
1636 Container::space_dimension>::cell_iterator
1637 coarse_cell = local_triangulation.
begin();
1638 coarse_cell != local_triangulation.
end();
1639 ++coarse_cell, ++index)
1641 patch_to_global_tria_map_tmp.insert(
1642 std::make_pair(coarse_cell, uniform_cells[index]));
1646 Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <=
1647 1
e-15 * coarse_cell->diameter(),
1650 bool refinement_necessary;
1655 refinement_necessary =
false;
1656 for (
const auto &active_tria_cell :
1659 if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
1661 active_tria_cell->set_refine_flag();
1662 refinement_necessary =
true;
1665 for (
unsigned int i = 0; i < patch.size(); ++i)
1670 if (patch_to_global_tria_map_tmp[active_tria_cell] ==
1675 for (
const unsigned int v :
1676 active_tria_cell->vertex_indices())
1677 active_tria_cell->vertex(v) = patch[i]->vertex(v);
1679 Assert(active_tria_cell->center().distance(
1680 patch_to_global_tria_map_tmp[active_tria_cell]
1682 1
e-15 * active_tria_cell->diameter(),
1685 active_tria_cell->set_user_flag();
1691 if (refinement_necessary)
1696 Container::dimension,
1697 Container::space_dimension>::cell_iterator cell =
1698 local_triangulation.
begin();
1699 cell != local_triangulation.
end();
1702 if (patch_to_global_tria_map_tmp.find(cell) !=
1703 patch_to_global_tria_map_tmp.end())
1705 if (cell->has_children())
1712 for (
unsigned int c = 0; c < cell->n_children(); ++c)
1714 if (patch_to_global_tria_map_tmp.find(cell->child(
1715 c)) == patch_to_global_tria_map_tmp.end())
1717 patch_to_global_tria_map_tmp.insert(
1720 patch_to_global_tria_map_tmp[cell]->child(
1742 patch_to_global_tria_map_tmp.erase(cell);
1748 while (refinement_necessary);
1754 Container::space_dimension>::cell_iterator
1755 cell = local_triangulation.
begin();
1756 cell != local_triangulation.
end();
1759 if (cell->user_flag_set())
1761 Assert(patch_to_global_tria_map_tmp.find(cell) !=
1762 patch_to_global_tria_map_tmp.end(),
1765 Assert(cell->center().distance(
1766 patch_to_global_tria_map_tmp[cell]->center()) <=
1767 1
e-15 * cell->diameter(),
1775 Container::space_dimension>::cell_iterator,
1776 typename Container::cell_iterator>::iterator
1777 map_tmp_it = patch_to_global_tria_map_tmp.
begin(),
1778 map_tmp_end = patch_to_global_tria_map_tmp.end();
1782 for (; map_tmp_it != map_tmp_end; ++map_tmp_it)
1783 patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
1788 template <
int dim,
int spacedim>
1791 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1805 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1806 dof_to_set_of_cells_map;
1808 std::vector<types::global_dof_index> local_dof_indices;
1809 std::vector<types::global_dof_index> local_face_dof_indices;
1810 std::vector<types::global_dof_index> local_line_dof_indices;
1814 std::vector<bool> user_flags;
1819 std::map<typename DoFHandler<dim, spacedim>::active_line_iterator,
1821 lines_to_parent_lines_map;
1828 .clear_user_flags();
1833 endc = dof_handler.
end();
1834 for (; cell != endc; ++cell)
1840 if (cell->is_artificial() ==
false)
1842 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
1843 if (cell->line(
l)->has_children())
1844 for (
unsigned int c = 0; c < cell->line(
l)->n_children();
1847 lines_to_parent_lines_map[cell->line(
l)->child(c)] =
1851 cell->line(
l)->child(c)->set_user_flag();
1866 endc = dof_handler.
end();
1867 for (; cell != endc; ++cell)
1872 if (cell->is_artificial() ==
false)
1874 const unsigned int n_dofs_per_cell =
1876 local_dof_indices.resize(n_dofs_per_cell);
1880 cell->get_dof_indices(local_dof_indices);
1881 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i)
1882 dof_to_set_of_cells_map[local_dof_indices[i]].
insert(cell);
1892 for (
const unsigned int f : cell->face_indices())
1894 if (cell->face(f)->has_children())
1896 for (
unsigned int c = 0; c < cell->face(f)->n_children();
1911 Assert(cell->face(f)->child(c)->has_children() ==
false,
1914 const unsigned int n_dofs_per_face =
1916 local_face_dof_indices.resize(n_dofs_per_face);
1918 cell->face(f)->child(c)->get_dof_indices(
1919 local_face_dof_indices);
1920 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1921 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1925 else if ((cell->face(f)->at_boundary() ==
false) &&
1926 (cell->neighbor_is_coarser(f)))
1943 std::pair<unsigned int, unsigned int>
1944 neighbor_face_no_subface_no =
1945 cell->neighbor_of_coarser_neighbor(f);
1946 unsigned int face_no = neighbor_face_no_subface_no.first;
1947 unsigned int subface = neighbor_face_no_subface_no.second;
1949 const unsigned int n_dofs_per_face =
1951 local_face_dof_indices.resize(n_dofs_per_face);
1953 cell->neighbor(f)->face(face_no)->get_dof_indices(
1954 local_face_dof_indices);
1955 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1956 dof_to_set_of_cells_map[local_face_dof_indices[i]].
insert(
1961 for (
unsigned int c = 0;
1962 c < cell->neighbor(f)->face(face_no)->n_children();
1968 const unsigned int n_dofs_per_face =
1970 local_face_dof_indices.resize(n_dofs_per_face);
1975 ->has_children() ==
false,
1980 ->get_dof_indices(local_face_dof_indices);
1981 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1982 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1999 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
2001 if (cell->line(
l)->has_children())
2003 for (
unsigned int c = 0;
2004 c < cell->line(
l)->n_children();
2007 Assert(cell->line(
l)->child(c)->has_children() ==
2013 const unsigned int n_dofs_per_line =
2016 local_line_dof_indices.resize(n_dofs_per_line);
2018 cell->line(
l)->child(c)->get_dof_indices(
2019 local_line_dof_indices);
2020 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2021 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2029 else if (cell->line(
l)->user_flag_set() ==
true)
2033 lines_to_parent_lines_map[cell->line(
l)];
2038 const unsigned int n_dofs_per_line =
2041 local_line_dof_indices.resize(n_dofs_per_line);
2043 parent_line->get_dof_indices(local_line_dof_indices);
2044 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2045 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2048 for (
unsigned int c = 0; c < parent_line->n_children();
2051 Assert(parent_line->child(c)->has_children() ==
2055 const unsigned int n_dofs_per_line =
2058 local_line_dof_indices.resize(n_dofs_per_line);
2060 parent_line->child(c)->get_dof_indices(
2061 local_line_dof_indices);
2062 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2063 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2080 .load_user_flags(user_flags);
2087 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2088 dof_to_cell_patches;
2092 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>::
2093 iterator it = dof_to_set_of_cells_map.begin(),
2094 it_end = dof_to_set_of_cells_map.end();
2095 for (; it != it_end; ++it)
2096 dof_to_cell_patches[it->first].assign(it->second.begin(),
2099 return dof_to_cell_patches;
2105 template <
typename CellIterator>
2108 std::set<std::pair<CellIterator, unsigned int>> &pairs1,
2111 const unsigned int direction,
2113 const ::Tensor<1, CellIterator::AccessorType::space_dimension>
2117 static const int space_dim = CellIterator::AccessorType::space_dimension;
2123 constexpr
int dim = CellIterator::AccessorType::dimension;
2124 constexpr
int spacedim = CellIterator::AccessorType::space_dimension;
2129 if (!(((pairs1.size() > 0) &&
2130 (
dynamic_cast<const parallel::fullydistributed::
2131 Triangulation<dim, spacedim> *
>(
2132 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2133 ((pairs2.size() > 0) &&
2136 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2137 Assert(pairs1.size() == pairs2.size(),
2138 ExcMessage(
"Unmatched faces on periodic boundaries"));
2142 unsigned int n_matches = 0;
2145 std::bitset<3> orientation;
2146 using PairIterator =
2147 typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator;
2148 for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
2150 for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
2152 const CellIterator cell1 = it1->first;
2153 const CellIterator cell2 = it2->first;
2154 const unsigned int face_idx1 = it1->second;
2155 const unsigned int face_idx2 = it2->second;
2157 cell1->face(face_idx1),
2158 cell2->face(face_idx2),
2167 {cell1, cell2}, {face_idx1, face_idx2}, orientation,
matrix};
2168 matched_pairs.push_back(matched_face);
2182 constexpr
int dim = CellIterator::AccessorType::dimension;
2183 constexpr
int spacedim = CellIterator::AccessorType::space_dimension;
2184 if (!(((pairs1.size() > 0) &&
2185 (
dynamic_cast<const parallel::fullydistributed::
2186 Triangulation<dim, spacedim> *
>(
2187 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2188 ((pairs2.size() > 0) &&
2191 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2192 AssertThrow(n_matches == pairs1.size() && pairs2.empty(),
2193 ExcMessage(
"Unmatched faces on periodic boundaries"));
2199 template <
typename MeshType>
2202 const MeshType &mesh,
2204 const
unsigned int direction,
2207 const
Tensor<1, MeshType::space_dimension> &offset,
2210 static const int dim = MeshType::dimension;
2211 static const int space_dim = MeshType::space_dimension;
2221 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2222 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2224 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2225 cell != mesh.end(0);
2228 const typename MeshType::face_iterator face_1 =
2229 cell->face(2 * direction);
2230 const typename MeshType::face_iterator face_2 =
2231 cell->face(2 * direction + 1);
2233 if (face_1->at_boundary() && face_1->boundary_id() == b_id)
2235 const std::pair<typename MeshType::cell_iterator, unsigned int>
2236 pair1 = std::make_pair(cell, 2 * direction);
2237 pairs1.insert(pair1);
2240 if (face_2->at_boundary() && face_2->boundary_id() == b_id)
2242 const std::pair<typename MeshType::cell_iterator, unsigned int>
2243 pair2 = std::make_pair(cell, 2 * direction + 1);
2244 pairs2.insert(pair2);
2248 Assert(pairs1.size() == pairs2.size(),
2249 ExcMessage(
"Unmatched faces on periodic boundaries"));
2251 Assert(pairs1.size() > 0,
2252 ExcMessage(
"No new periodic face pairs have been found. "
2253 "Are you sure that you've selected the correct boundary "
2254 "id's and that the coarsest level mesh is colorized?"));
2257 const unsigned int size_old = matched_pairs.size();
2262 pairs1, pairs2, direction, matched_pairs, offset,
matrix);
2266 const unsigned int size_new = matched_pairs.size();
2267 for (
unsigned int i = size_old; i < size_new; ++i)
2269 Assert(matched_pairs[i].orientation == 1,
2271 "Found a face match with non standard orientation. "
2272 "This function is only suitable for meshes with cells "
2273 "in default orientation"));
2280 template <
typename MeshType>
2283 const MeshType &mesh,
2286 const
unsigned int direction,
2289 const
Tensor<1, MeshType::space_dimension> &offset,
2292 static const int dim = MeshType::dimension;
2293 static const int space_dim = MeshType::space_dimension;
2301 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2302 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2304 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2305 cell != mesh.end(0);
2308 for (
const unsigned int i : cell->face_indices())
2310 const typename MeshType::face_iterator face = cell->face(i);
2311 if (face->at_boundary() && face->boundary_id() == b_id1)
2313 const std::pair<typename MeshType::cell_iterator, unsigned int>
2314 pair1 = std::make_pair(cell, i);
2315 pairs1.insert(pair1);
2318 if (face->at_boundary() && face->boundary_id() == b_id2)
2320 const std::pair<typename MeshType::cell_iterator, unsigned int>
2321 pair2 = std::make_pair(cell, i);
2322 pairs2.insert(pair2);
2333 if (!(((pairs1.size() > 0) &&
2336 *
>(&pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2337 ((pairs2.size() > 0) &&
2340 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2341 Assert(pairs1.size() == pairs2.size(),
2342 ExcMessage(
"Unmatched faces on periodic boundaries"));
2345 (pairs1.size() > 0 ||
2348 &mesh.begin()->get_triangulation()) !=
nullptr)),
2349 ExcMessage(
"No new periodic face pairs have been found. "
2350 "Are you sure that you've selected the correct boundary "
2351 "id's and that the coarsest level mesh is colorized?"));
2355 pairs1, pairs2, direction, matched_pairs, offset,
matrix);
2369 template <
int spacedim>
2373 const unsigned int direction,
2383 if (
matrix.m() == spacedim)
2384 for (
unsigned int i = 0; i < spacedim; ++i)
2385 for (
unsigned int j = 0; j < spacedim; ++j)
2386 distance(i) +=
matrix(i, j) * point1(j);
2390 distance += offset - point2;
2392 for (
unsigned int i = 0; i < spacedim; ++i)
2398 if (std::abs(distance(i)) > 1.e-10)
2424 std::array<unsigned int, GeometryInfo<1>::vertices_per_face>;
2425 static inline std::bitset<3>
2437 std::array<unsigned int, GeometryInfo<2>::vertices_per_face>;
2438 static inline std::bitset<3>
2446 static const MATCH_T m_tff = {{0, 1}};
2447 if (matching == m_tff)
2449 static const MATCH_T m_ttf = {{1, 0}};
2450 if (matching == m_ttf)
2463 std::array<unsigned int, GeometryInfo<3>::vertices_per_face>;
2465 static inline std::bitset<3>
2473 static const MATCH_T m_tff = {{0, 1, 2, 3}};
2474 if (matching == m_tff)
2476 static const MATCH_T m_tft = {{1, 3, 0, 2}};
2477 if (matching == m_tft)
2479 static const MATCH_T m_ttf = {{3, 2, 1, 0}};
2480 if (matching == m_ttf)
2482 static const MATCH_T m_ttt = {{2, 0, 3, 1}};
2483 if (matching == m_ttt)
2485 static const MATCH_T m_fff = {{0, 2, 1, 3}};
2486 if (matching == m_fff)
2488 static const MATCH_T m_fft = {{2, 3, 0, 1}};
2489 if (matching == m_fft)
2491 static const MATCH_T m_ftf = {{3, 1, 2, 0}};
2492 if (matching == m_ftf)
2494 static const MATCH_T m_ftt = {{1, 0, 3, 2}};
2495 if (matching == m_ftt)
2506 template <
typename FaceIterator>
2509 std::bitset<3> &orientation,
2510 const FaceIterator &face1,
2511 const FaceIterator &face2,
2512 const unsigned int direction,
2517 ExcMessage(
"The supplied matrix must be a square matrix"));
2519 static const int dim = FaceIterator::AccessorType::dimension;
2523 std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
2527 std::set<unsigned int> face2_vertices;
2528 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2529 face2_vertices.insert(i);
2531 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2533 for (std::set<unsigned int>::iterator it = face2_vertices.begin();
2534 it != face2_vertices.end();
2544 face2_vertices.erase(it);
2551 if (face2_vertices.empty())
2554 return face2_vertices.empty();
2559 template <
typename FaceIterator>
2562 const FaceIterator &face1,
2563 const FaceIterator &face2,
2564 const unsigned int direction,
2569 std::bitset<3> dummy;
2578 #include "grid_tools_dof_handlers.inst"
cell_iterator end() const
const Triangulation< dim, spacedim > & get_triangulation() const
active_cell_iterator begin_active(const unsigned int level=0) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() 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.
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
const std::vector< bool > & get_used_vertices() const
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
cell_iterator begin(const unsigned int level=0) const
unsigned int n_active_cells() const
void save_user_flags(std::ostream &out) const
cell_iterator end() const
virtual void execute_coarsening_and_refinement()
unsigned int n_cells() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int n_vertices() const
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
concept is_triangulation_or_dof_handler
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
static const unsigned int invalid_unsigned_int
typename type_identity< T >::type type_identity_t
unsigned int global_dof_index
const ::Triangulation< dim, spacedim > & tria