55 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
59 const std::vector<bool> & marked_vertices)
70 marked_vertices.size() == 0,
72 marked_vertices.size()));
82 marked_vertices.size() == 0 ||
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 = (marked_vertices.size() == 0) ?
103 std::vector<bool>::const_iterator
first =
104 std::find(used.begin(), used.end(),
true);
110 unsigned int best_vertex = std::distance(used.begin(),
first);
111 double best_dist = (p -
vertices[best_vertex]).norm_square();
115 for (
unsigned int j = best_vertex + 1; j <
vertices.size(); ++j)
118 double dist = (p -
vertices[j]).norm_square();
119 if (dist < best_dist)
131 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
134 const MeshType<dim, spacedim> &mesh,
136 const std::vector<bool> & marked_vertices)
151 marked_vertices.size() == 0,
153 marked_vertices.size()));
163 marked_vertices.size() == 0 ||
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 <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
193 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
196 typename ::internal::
197 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
200 const unsigned int vertex)
206 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
212 std::set<typename ::internal::
213 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
216 typename ::internal::
217 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
218 cell = mesh.begin_active(),
254 for (; cell != endc; ++cell)
256 for (
const unsigned int v : cell->vertex_indices())
257 if (cell->vertex_index(v) == vertex)
269 for (
const auto face :
270 cell->reference_cell().faces_for_given_vertex(v))
271 if (!cell->at_boundary(face) &&
272 cell->neighbor(face)->is_active())
294 for (
unsigned int e = 0;
e < cell->n_lines(); ++
e)
295 if (cell->line(e)->has_children())
299 if (cell->line(e)->child(0)->vertex_index(1) == vertex)
323 typename ::internal::
324 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
332 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
334 find_active_cell_around_point_internal(
335 const MeshType<dim, spacedim> &mesh,
337 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>
339 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>
343 typename ::internal::
344 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
347 typename ::internal::
348 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
353 using cell_iterator =
354 typename MeshType<dim, spacedim>::active_cell_iterator;
356 using cell_iterator = typename ::internal::
357 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
365 std::set<cell_iterator> adjacent_cells_new;
367 typename std::set<cell_iterator>::const_iterator cell =
371 for (; cell != endc; ++cell)
373 std::vector<cell_iterator> active_neighbors;
374 get_active_neighbors<MeshType<dim, spacedim>>(*cell,
376 for (
unsigned int i = 0; i < active_neighbors.size(); ++i)
377 if (searched_cells.find(active_neighbors[i]) ==
378 searched_cells.end())
379 adjacent_cells_new.insert(active_neighbors[i]);
383 adjacent_cells_new.end());
392 cell_iterator it = mesh.begin_active();
393 for (; it != mesh.end(); ++it)
394 if (searched_cells.find(it) == searched_cells.end())
405 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
407 typename MeshType<dim, spacedim>::active_cell_iterator
409 typename ::internal::
410 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
414 const std::vector<bool> & marked_vertices,
415 const double tolerance)
417 return find_active_cell_around_point<dim, MeshType, spacedim>(
428 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
430 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
432 std::pair<typename ::internal::
433 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
437 const MeshType<dim, spacedim> &mesh,
439 const std::vector<bool> & marked_vertices,
440 const double tolerance)
442 using active_cell_iterator = typename ::internal::
443 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
449 double best_distance = tolerance;
451 std::pair<active_cell_iterator, Point<dim>> best_cell;
454 best_cell.first = mesh.end();
458 std::vector<active_cell_iterator> adjacent_cells_tmp =
467 std::set<active_cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
468 adjacent_cells_tmp.end());
469 std::set<active_cell_iterator> searched_cells;
477 const unsigned int n_active_cells =
478 mesh.get_triangulation().n_active_cells();
480 unsigned int cells_searched = 0;
481 while (!found && cells_searched < n_active_cells)
483 typename std::set<active_cell_iterator>::const_iterator
486 for (; cell != endc; ++cell)
488 if ((*cell)->is_artificial() ==
false)
504 if ((dist < best_distance) ||
505 ((dist == best_distance) &&
506 ((*cell)->level() > best_level)))
509 best_distance = dist;
510 best_level = (*cell)->level();
511 best_cell = std::make_pair(*cell, p_cell);
538 if (marked_vertices.size() > 0)
539 cells_searched = n_active_cells;
548 if (!found && cells_searched < n_active_cells)
550 find_active_cell_around_point_internal<dim, MeshType, spacedim>(
560 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
562 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
565 std::vector<std::pair<
566 typename ::internal::
567 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
571 const MeshType<dim, spacedim> &mesh,
573 const double tolerance,
574 const std::vector<bool> &marked_vertices)
577 mapping, mesh, p, marked_vertices, tolerance);
579 if (cell_and_point.first == mesh.end())
583 mapping, mesh, p, tolerance, cell_and_point);
588 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
590 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
593 std::vector<std::pair<
594 typename ::internal::
595 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
600 const MeshType<dim, spacedim> &mesh,
602 const double tolerance,
603 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
607 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
612 cells_and_points.push_back(first_cell);
616 const Point<dim> unit_point = cells_and_points.front().second;
617 const auto my_cell = cells_and_points.front().first;
619 unsigned int n_dirs_at_threshold = 0;
621 for (
unsigned int d = 0; d < dim; ++d)
623 distance_to_center[d] =
std::abs(unit_point[d] - 0.5);
624 if (distance_to_center[d] > 0.5 - tolerance)
626 ++n_dirs_at_threshold;
627 last_point_at_threshold = d;
631 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
634 if (n_dirs_at_threshold == 1)
636 unsigned int neighbor_index =
637 2 * last_point_at_threshold +
638 (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
639 if (!my_cell->at_boundary(neighbor_index))
640 cells_to_add.push_back(my_cell->neighbor(neighbor_index));
643 else if (n_dirs_at_threshold == dim)
645 unsigned int local_vertex_index = 0;
646 for (
unsigned int d = 0; d < dim; ++d)
647 local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d;
648 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
650 mesh, my_cell->vertex_index(local_vertex_index));
651 for (
const auto &cell : cells)
653 cells_to_add.push_back(cell);
660 else if (n_dirs_at_threshold == 2)
663 unsigned int count_vertex_indices = 0;
665 for (
unsigned int d = 0; d < dim; ++d)
667 if (distance_to_center[d] > 0.5 - tolerance)
671 unit_point[d] > 0.5 ? 1 : 0;
672 count_vertex_indices++;
682 const unsigned int first_vertex =
685 for (
unsigned int d = 0; d < 2; ++d)
689 my_cell->vertex_index(first_vertex + (d << free_direction)));
690 for (
const auto &cell : tentative_cells)
692 bool cell_not_yet_present =
true;
693 for (
const auto &other_cell : cells_to_add)
694 if (cell == other_cell)
696 cell_not_yet_present =
false;
699 if (cell_not_yet_present)
700 cells_to_add.push_back(cell);
705 const double original_distance_to_unit_cell =
707 for (
const auto &cell : cells_to_add)
715 original_distance_to_unit_cell + tolerance)
716 cells_and_points.emplace_back(cell, p_unit);
723 cells_and_points.begin(),
724 cells_and_points.end(),
725 [](
const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
727 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
728 Point<dim>> &b) { return a.first < b.first; });
730 return cells_and_points;
735 template <
class MeshType>
736 std::vector<typename MeshType::active_cell_iterator>
738 const MeshType &mesh,
739 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
742 std::vector<typename MeshType::active_cell_iterator> active_halo_layer;
743 std::vector<bool> locally_active_vertices_on_subdomain(
744 mesh.get_triangulation().n_vertices(),
false);
749 for (
const auto &cell : mesh.active_cell_iterators())
751 for (
const auto v : cell->vertex_indices())
752 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
true;
757 for (
const auto &cell : mesh.active_cell_iterators())
758 if (!predicate(cell))
759 for (
const auto v : cell->vertex_indices())
760 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
763 active_halo_layer.push_back(cell);
767 return active_halo_layer;
772 template <
class MeshType>
773 std::vector<typename MeshType::cell_iterator>
775 const MeshType &mesh,
776 const std::function<
bool(
const typename MeshType::cell_iterator &)>
778 const unsigned int level)
780 std::vector<typename MeshType::cell_iterator> level_halo_layer;
781 std::vector<bool> locally_active_vertices_on_level_subdomain(
782 mesh.get_triangulation().n_vertices(),
false);
787 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
788 cell != mesh.end(
level);
791 for (
const unsigned int v : cell->vertex_indices())
792 locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
798 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
799 cell != mesh.end(
level);
801 if (!predicate(cell))
802 for (
const unsigned int v : cell->vertex_indices())
803 if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
806 level_halo_layer.push_back(cell);
810 return level_halo_layer;
816 template <
class MeshType>
818 contains_locally_owned_cells(
819 const std::vector<typename MeshType::active_cell_iterator> &cells)
821 for (
typename std::vector<
822 typename MeshType::active_cell_iterator>::const_iterator it =
827 if ((*it)->is_locally_owned())
833 template <
class MeshType>
835 contains_artificial_cells(
836 const std::vector<typename MeshType::active_cell_iterator> &cells)
838 for (
typename std::vector<
839 typename MeshType::active_cell_iterator>::const_iterator it =
844 if ((*it)->is_artificial())
853 template <
class MeshType>
854 std::vector<typename MeshType::active_cell_iterator>
857 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
860 const std::vector<typename MeshType::active_cell_iterator>
865 Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) ==
false,
866 ExcMessage(
"Halo layer contains locally owned cells"));
867 Assert(contains_artificial_cells<MeshType>(active_halo_layer) ==
false,
868 ExcMessage(
"Halo layer contains artificial cells"));
870 return active_halo_layer;
875 template <
class MeshType>
876 std::vector<typename MeshType::active_cell_iterator>
878 const MeshType &mesh,
879 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
881 const double layer_thickness)
883 std::vector<typename MeshType::active_cell_iterator>
884 subdomain_boundary_cells, active_cell_layer_within_distance;
885 std::vector<bool> vertices_outside_subdomain(
886 mesh.get_triangulation().n_vertices(),
false);
888 const unsigned int spacedim = MeshType::space_dimension;
890 unsigned int n_non_predicate_cells = 0;
898 for (
const auto &cell : mesh.active_cell_iterators())
899 if (!predicate(cell))
901 for (
const unsigned int v : cell->vertex_indices())
902 vertices_outside_subdomain[cell->vertex_index(v)] =
true;
903 n_non_predicate_cells++;
910 if (n_non_predicate_cells == 0 ||
911 n_non_predicate_cells == mesh.get_triangulation().n_active_cells())
912 return std::vector<typename MeshType::active_cell_iterator>();
916 for (
const auto &cell : mesh.active_cell_iterators())
919 for (
const unsigned int v : cell->vertex_indices())
920 if (vertices_outside_subdomain[cell->vertex_index(v)] ==
true)
922 subdomain_boundary_cells.push_back(cell);
932 const double DOUBLE_EPSILON = 100. * std::numeric_limits<double>::epsilon();
935 for (
unsigned int d = 0; d < spacedim; ++d)
937 bounding_box.first[d] -= (layer_thickness + DOUBLE_EPSILON);
938 bounding_box.second[d] += (layer_thickness + DOUBLE_EPSILON);
941 std::vector<Point<spacedim>>
942 subdomain_boundary_cells_centers;
945 subdomain_boundary_cells_radii;
947 subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size());
948 subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size());
950 for (
typename std::vector<typename MeshType::active_cell_iterator>::
951 const_iterator subdomain_boundary_cell_iterator =
952 subdomain_boundary_cells.begin();
953 subdomain_boundary_cell_iterator != subdomain_boundary_cells.end();
954 ++subdomain_boundary_cell_iterator)
956 const std::pair<Point<spacedim>,
double>
957 &subdomain_boundary_cell_enclosing_ball =
958 (*subdomain_boundary_cell_iterator)->enclosing_ball();
960 subdomain_boundary_cells_centers.push_back(
961 subdomain_boundary_cell_enclosing_ball.first);
962 subdomain_boundary_cells_radii.push_back(
963 subdomain_boundary_cell_enclosing_ball.second);
965 AssertThrow(subdomain_boundary_cells_radii.size() ==
966 subdomain_boundary_cells_centers.size(),
975 for (
const auto &cell : mesh.active_cell_iterators())
981 const std::pair<Point<spacedim>,
double> &cell_enclosing_ball =
982 cell->enclosing_ball();
985 cell_enclosing_ball.first;
986 const double cell_enclosing_ball_radius = cell_enclosing_ball.second;
988 bool cell_inside =
true;
990 for (
unsigned int d = 0; d < spacedim; ++d)
992 (cell_enclosing_ball_center[d] + cell_enclosing_ball_radius >
993 bounding_box.first[d]) &&
994 (cell_enclosing_ball_center[d] - cell_enclosing_ball_radius <
995 bounding_box.second[d]);
1001 for (
unsigned int i = 0; i < subdomain_boundary_cells_radii.size();
1004 subdomain_boundary_cells_centers[i]) <
1005 Utilities::fixed_power<2>(cell_enclosing_ball_radius +
1006 subdomain_boundary_cells_radii[i] +
1007 layer_thickness + DOUBLE_EPSILON))
1009 active_cell_layer_within_distance.push_back(cell);
1014 return active_cell_layer_within_distance;
1019 template <
class MeshType>
1020 std::vector<typename MeshType::active_cell_iterator>
1022 const double layer_thickness)
1025 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1026 predicate(locally_owned_cell_predicate);
1028 const std::vector<typename MeshType::active_cell_iterator>
1029 ghost_cell_layer_within_distance =
1037 contains_locally_owned_cells<MeshType>(
1038 ghost_cell_layer_within_distance) ==
false,
1040 "Ghost cells within layer_thickness contains locally owned cells."));
1042 contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) ==
1045 "Ghost cells within layer_thickness contains artificial cells. "
1046 "The function compute_ghost_cell_layer_within_distance "
1047 "is probably called while using parallel::distributed::Triangulation. "
1048 "In such case please refer to the description of this function."));
1050 return ghost_cell_layer_within_distance;
1055 template <
class MeshType>
1058 const MeshType &mesh,
1059 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1062 std::vector<bool> locally_active_vertices_on_subdomain(
1063 mesh.get_triangulation().n_vertices(),
false);
1065 const unsigned int spacedim = MeshType::space_dimension;
1072 for (
const auto &cell : mesh.active_cell_iterators())
1073 if (predicate(cell))
1075 minp = cell->center();
1076 maxp = cell->center();
1082 for (
const auto &cell : mesh.active_cell_iterators())
1083 if (predicate(cell))
1084 for (
const unsigned int v : cell->vertex_indices())
1085 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
1088 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
1090 for (
unsigned int d = 0; d < spacedim; ++d)
1092 minp[d] =
std::min(minp[d], cell->vertex(v)[d]);
1093 maxp[d] =
std::max(maxp[d], cell->vertex(v)[d]);
1097 return std::make_pair(minp, maxp);
1102 template <
typename MeshType>
1103 std::list<std::pair<
typename MeshType::cell_iterator,
1104 typename MeshType::cell_iterator>>
1108 ExcMessage(
"The two meshes must be represent triangulations that "
1109 "have the same coarse meshes"));
1116 bool remove_ghost_cells =
false;
1117#ifdef DEAL_II_WITH_MPI
1119 constexpr int dim = MeshType::dimension;
1120 constexpr int spacedim = MeshType::space_dimension;
1122 *
>(&mesh_1.get_triangulation()) !=
nullptr ||
1124 *
>(&mesh_2.get_triangulation()) !=
nullptr)
1126 Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(),
1127 ExcMessage(
"This function can only be used with meshes "
1128 "corresponding to distributed Triangulations when "
1129 "both Triangulations are equal."));
1130 remove_ghost_cells =
true;
1142 using CellList = std::list<std::pair<
typename MeshType::cell_iterator,
1143 typename MeshType::cell_iterator>>;
1147 typename MeshType::cell_iterator cell_1 = mesh_1.begin(0),
1148 cell_2 = mesh_2.begin(0);
1149 for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
1150 cell_list.emplace_back(cell_1, cell_2);
1153 typename CellList::iterator cell_pair = cell_list.begin();
1154 while (cell_pair != cell_list.end())
1158 if (cell_pair->first->has_children() &&
1159 cell_pair->second->has_children())
1161 Assert(cell_pair->first->refinement_case() ==
1162 cell_pair->second->refinement_case(),
1164 for (
unsigned int c = 0; c < cell_pair->first->n_children(); ++c)
1165 cell_list.emplace_back(cell_pair->first->child(c),
1166 cell_pair->second->child(c));
1171 const auto previous_cell_pair = cell_pair;
1173 cell_list.erase(previous_cell_pair);
1178 if (remove_ghost_cells &&
1179 ((cell_pair->first->is_active() &&
1180 !cell_pair->first->is_locally_owned()) ||
1181 (cell_pair->second->is_active() &&
1182 !cell_pair->second->is_locally_owned())))
1185 const auto previous_cell_pair = cell_pair;
1187 cell_list.erase(previous_cell_pair);
1196 for (cell_pair = cell_list.begin(); cell_pair != cell_list.end();
1198 Assert(cell_pair->first->is_active() || cell_pair->second->is_active() ||
1199 (cell_pair->first->refinement_case() !=
1200 cell_pair->second->refinement_case()),
1208 template <
int dim,
int spacedim>
1225 endc = mesh_1.
end(0);
1226 for (; cell_1 != endc; ++cell_1, ++cell_2)
1228 if (cell_1->n_vertices() != cell_2->n_vertices())
1230 for (
const unsigned int v : cell_1->vertex_indices())
1231 if (cell_1->vertex(v) != cell_2->vertex(v))
1244 template <
typename MeshType>
1249 mesh_2.get_triangulation());
1254 template <
int dim,
int spacedim>
1255 std::pair<typename ::DoFHandler<dim, spacedim>::active_cell_iterator,
1261 const double tolerance)
1265 ExcMessage(
"Mapping collection needs to have either size 1 "
1266 "or size equal to the number of elements in "
1267 "the FECollection."));
1269 using cell_iterator =
1270 typename ::DoFHandler<dim, spacedim>::active_cell_iterator;
1272 std::pair<cell_iterator, Point<dim>> best_cell;
1276 if (mapping.
size() == 1)
1278 const std::vector<bool> marked_vertices = {};
1280 mapping[0], mesh, p, marked_vertices, tolerance);
1287 double best_distance = tolerance;
1288 int best_level = -1;
1295 std::vector<cell_iterator> adjacent_cells_tmp =
1303 std::set<cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
1304 adjacent_cells_tmp.end());
1305 std::set<cell_iterator> searched_cells;
1315 unsigned int cells_searched = 0;
1316 while (!found && cells_searched < n_cells)
1318 typename std::set<cell_iterator>::const_iterator
1321 for (; cell != endc; ++cell)
1326 mapping[(*cell)->active_fe_index()]
1327 .transform_real_to_unit_cell(*cell, p);
1339 if (dist < best_distance || (dist == best_distance &&
1340 (*cell)->level() > best_level))
1343 best_distance = dist;
1344 best_level = (*cell)->level();
1345 best_cell = std::make_pair(*cell, p_cell);
1371 if (!found && cells_searched < n_cells)
1373 find_active_cell_around_point_internal<dim,
1385 template <
class MeshType>
1386 std::vector<typename MeshType::active_cell_iterator>
1389 Assert(cell->is_locally_owned(),
1390 ExcMessage(
"This function only makes sense if the cell for "
1391 "which you are asking for a patch, is locally "
1394 std::vector<typename MeshType::active_cell_iterator> patch;
1395 patch.push_back(cell);
1396 for (
const unsigned int face_number : cell->face_indices())
1397 if (cell->face(face_number)->at_boundary() ==
false)
1399 if (cell->neighbor(face_number)->has_children() ==
false)
1400 patch.push_back(cell->neighbor(face_number));
1405 if (MeshType::dimension > 1)
1407 for (
unsigned int subface = 0;
1408 subface < cell->face(face_number)->n_children();
1411 cell->neighbor_child_on_subface(face_number, subface));
1417 typename MeshType::cell_iterator neighbor =
1418 cell->neighbor(face_number);
1419 while (neighbor->has_children())
1420 neighbor = neighbor->child(1 - face_number);
1422 Assert(neighbor->neighbor(1 - face_number) == cell,
1424 patch.push_back(neighbor);
1432 template <
class Container>
1433 std::vector<typename Container::cell_iterator>
1435 const std::vector<typename Container::active_cell_iterator> &patch)
1439 "Vector containing patch cells should not be an empty vector!"));
1443 int min_level = patch[0]->level();
1445 for (
unsigned int i = 0; i < patch.size(); ++i)
1447 std::set<typename Container::cell_iterator> uniform_cells;
1448 typename std::vector<
1449 typename Container::active_cell_iterator>::const_iterator patch_cell;
1451 for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell)
1456 if ((*patch_cell)->level() == min_level)
1457 uniform_cells.insert(*patch_cell);
1464 typename Container::cell_iterator parent = *patch_cell;
1466 while (parent->level() > min_level)
1467 parent = parent->parent();
1468 uniform_cells.insert(parent);
1472 return std::vector<typename Container::cell_iterator>(uniform_cells.begin(),
1473 uniform_cells.end());
1478 template <
class Container>
1481 const std::vector<typename Container::active_cell_iterator> &patch,
1483 &local_triangulation,
1486 Container::space_dimension>::active_cell_iterator,
1487 typename Container::active_cell_iterator> &patch_to_global_tria_map)
1490 const std::vector<typename Container::cell_iterator> uniform_cells =
1491 get_cells_at_coarsest_common_level<Container>(patch);
1493 local_triangulation.
clear();
1494 std::vector<Point<Container::space_dimension>>
vertices;
1495 const unsigned int n_uniform_cells = uniform_cells.size();
1496 std::vector<CellData<Container::dimension>> cells(n_uniform_cells);
1499 typename std::vector<typename Container::cell_iterator>::const_iterator
1501 for (uniform_cell = uniform_cells.begin();
1502 uniform_cell != uniform_cells.end();
1505 for (
const unsigned int v : (*uniform_cell)->vertex_indices())
1508 (*uniform_cell)->vertex(v);
1509 bool repeat_vertex =
false;
1511 for (
unsigned int m = 0; m < i; ++m)
1515 repeat_vertex =
true;
1516 cells[k].vertices[v] = m;
1520 if (repeat_vertex ==
false)
1523 cells[k].vertices[v] = i;
1534 unsigned int index = 0;
1537 Container::space_dimension>::cell_iterator,
1538 typename Container::cell_iterator>
1539 patch_to_global_tria_map_tmp;
1541 Container::space_dimension>::cell_iterator
1542 coarse_cell = local_triangulation.
begin();
1543 coarse_cell != local_triangulation.
end();
1544 ++coarse_cell, ++index)
1546 patch_to_global_tria_map_tmp.insert(
1547 std::make_pair(coarse_cell, uniform_cells[index]));
1551 Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <=
1552 1e-15 * coarse_cell->diameter(),
1555 bool refinement_necessary;
1560 refinement_necessary =
false;
1561 for (
const auto &active_tria_cell :
1564 if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
1566 active_tria_cell->set_refine_flag();
1567 refinement_necessary =
true;
1570 for (
unsigned int i = 0; i < patch.size(); ++i)
1575 if (patch_to_global_tria_map_tmp[active_tria_cell] ==
1580 for (
const unsigned int v :
1581 active_tria_cell->vertex_indices())
1582 active_tria_cell->vertex(v) = patch[i]->vertex(v);
1584 Assert(active_tria_cell->center().distance(
1585 patch_to_global_tria_map_tmp[active_tria_cell]
1587 1e-15 * active_tria_cell->diameter(),
1590 active_tria_cell->set_user_flag();
1596 if (refinement_necessary)
1601 Container::dimension,
1602 Container::space_dimension>::cell_iterator cell =
1603 local_triangulation.
begin();
1604 cell != local_triangulation.
end();
1607 if (patch_to_global_tria_map_tmp.find(cell) !=
1608 patch_to_global_tria_map_tmp.end())
1610 if (cell->has_children())
1617 for (
unsigned int c = 0; c < cell->n_children(); ++c)
1619 if (patch_to_global_tria_map_tmp.find(cell->child(
1620 c)) == patch_to_global_tria_map_tmp.end())
1622 patch_to_global_tria_map_tmp.insert(
1625 patch_to_global_tria_map_tmp[cell]->child(
1647 patch_to_global_tria_map_tmp.erase(cell);
1653 while (refinement_necessary);
1659 Container::space_dimension>::cell_iterator
1660 cell = local_triangulation.
begin();
1661 cell != local_triangulation.
end();
1664 if (cell->user_flag_set())
1666 Assert(patch_to_global_tria_map_tmp.find(cell) !=
1667 patch_to_global_tria_map_tmp.end(),
1670 Assert(cell->center().distance(
1671 patch_to_global_tria_map_tmp[cell]->center()) <=
1672 1e-15 * cell->diameter(),
1680 Container::space_dimension>::cell_iterator,
1681 typename Container::cell_iterator>::iterator
1682 map_tmp_it = patch_to_global_tria_map_tmp.
begin(),
1683 map_tmp_end = patch_to_global_tria_map_tmp.end();
1687 for (; map_tmp_it != map_tmp_end; ++map_tmp_it)
1688 patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
1693 template <
int dim,
int spacedim>
1696 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1710 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1711 dof_to_set_of_cells_map;
1713 std::vector<types::global_dof_index> local_dof_indices;
1714 std::vector<types::global_dof_index> local_face_dof_indices;
1715 std::vector<types::global_dof_index> local_line_dof_indices;
1719 std::vector<bool> user_flags;
1724 std::map<typename DoFHandler<dim, spacedim>::active_line_iterator,
1726 lines_to_parent_lines_map;
1733 .clear_user_flags();
1738 endc = dof_handler.
end();
1739 for (; cell != endc; ++cell)
1745 if (cell->is_artificial() ==
false)
1747 for (
unsigned int l = 0; l < cell->n_lines(); ++l)
1748 if (cell->line(l)->has_children())
1749 for (
unsigned int c = 0; c < cell->line(l)->n_children();
1752 lines_to_parent_lines_map[cell->line(l)->child(c)] =
1756 cell->line(l)->child(c)->set_user_flag();
1771 endc = dof_handler.
end();
1772 for (; cell != endc; ++cell)
1777 if (cell->is_artificial() ==
false)
1779 const unsigned int n_dofs_per_cell =
1781 local_dof_indices.resize(n_dofs_per_cell);
1785 cell->get_dof_indices(local_dof_indices);
1786 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i)
1787 dof_to_set_of_cells_map[local_dof_indices[i]].insert(cell);
1797 for (
const unsigned int f : cell->face_indices())
1799 if (cell->face(f)->has_children())
1801 for (
unsigned int c = 0; c < cell->face(f)->n_children();
1816 Assert(cell->face(f)->child(c)->has_children() ==
false,
1819 const unsigned int n_dofs_per_face =
1821 local_face_dof_indices.resize(n_dofs_per_face);
1823 cell->face(f)->child(c)->get_dof_indices(
1824 local_face_dof_indices);
1825 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1826 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1830 else if ((cell->face(f)->at_boundary() ==
false) &&
1831 (cell->neighbor_is_coarser(f)))
1848 std::pair<unsigned int, unsigned int>
1849 neighbor_face_no_subface_no =
1850 cell->neighbor_of_coarser_neighbor(f);
1851 unsigned int face_no = neighbor_face_no_subface_no.first;
1852 unsigned int subface = neighbor_face_no_subface_no.second;
1854 const unsigned int n_dofs_per_face =
1856 local_face_dof_indices.resize(n_dofs_per_face);
1858 cell->neighbor(f)->face(face_no)->get_dof_indices(
1859 local_face_dof_indices);
1860 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1861 dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(
1866 for (
unsigned int c = 0;
1867 c < cell->neighbor(f)->face(face_no)->n_children();
1873 const unsigned int n_dofs_per_face =
1875 local_face_dof_indices.resize(n_dofs_per_face);
1880 ->has_children() ==
false,
1885 ->get_dof_indices(local_face_dof_indices);
1886 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1887 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1904 for (
unsigned int l = 0; l < cell->n_lines(); ++l)
1906 if (cell->line(l)->has_children())
1908 for (
unsigned int c = 0;
1909 c < cell->line(l)->n_children();
1912 Assert(cell->line(l)->child(c)->has_children() ==
1918 const unsigned int n_dofs_per_line =
1921 local_line_dof_indices.resize(n_dofs_per_line);
1923 cell->line(l)->child(c)->get_dof_indices(
1924 local_line_dof_indices);
1925 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1926 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1934 else if (cell->line(l)->user_flag_set() ==
true)
1938 lines_to_parent_lines_map[cell->line(l)];
1943 const unsigned int n_dofs_per_line =
1946 local_line_dof_indices.resize(n_dofs_per_line);
1948 parent_line->get_dof_indices(local_line_dof_indices);
1949 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1950 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1953 for (
unsigned int c = 0; c < parent_line->n_children();
1956 Assert(parent_line->child(c)->has_children() ==
1960 const unsigned int n_dofs_per_line =
1963 local_line_dof_indices.resize(n_dofs_per_line);
1965 parent_line->child(c)->get_dof_indices(
1966 local_line_dof_indices);
1967 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1968 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1985 .load_user_flags(user_flags);
1992 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1993 dof_to_cell_patches;
1997 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>::
1998 iterator it = dof_to_set_of_cells_map.begin(),
1999 it_end = dof_to_set_of_cells_map.end();
2000 for (; it != it_end; ++it)
2001 dof_to_cell_patches[it->first].assign(it->second.begin(),
2004 return dof_to_cell_patches;
2010 template <
typename CellIterator>
2013 std::set<std::pair<CellIterator, unsigned int>> &pairs1,
2016 const unsigned int direction,
2018 const ::Tensor<1, CellIterator::AccessorType::space_dimension>
2022 static const int space_dim = CellIterator::AccessorType::space_dimension;
2028 constexpr int dim = CellIterator::AccessorType::dimension;
2029 constexpr int spacedim = CellIterator::AccessorType::space_dimension;
2034 if (!(((pairs1.size() > 0) &&
2035 (
dynamic_cast<const parallel::fullydistributed::
2036 Triangulation<dim, spacedim> *
>(
2037 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2038 ((pairs2.size() > 0) &&
2041 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2042 Assert(pairs1.size() == pairs2.size(),
2043 ExcMessage(
"Unmatched faces on periodic boundaries"));
2047 unsigned int n_matches = 0;
2050 std::bitset<3> orientation;
2051 using PairIterator =
2052 typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator;
2053 for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
2055 for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
2057 const CellIterator cell1 = it1->first;
2058 const CellIterator cell2 = it2->first;
2059 const unsigned int face_idx1 = it1->second;
2060 const unsigned int face_idx2 = it2->second;
2062 cell1->face(face_idx1),
2063 cell2->face(face_idx2),
2072 {cell1, cell2}, {face_idx1, face_idx2}, orientation, matrix};
2073 matched_pairs.push_back(matched_face);
2087 constexpr int dim = CellIterator::AccessorType::dimension;
2088 constexpr int spacedim = CellIterator::AccessorType::space_dimension;
2089 if (!(((pairs1.size() > 0) &&
2090 (
dynamic_cast<const parallel::fullydistributed::
2091 Triangulation<dim, spacedim> *
>(
2092 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2093 ((pairs2.size() > 0) &&
2096 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2097 AssertThrow(n_matches == pairs1.size() && pairs2.size() == 0,
2098 ExcMessage(
"Unmatched faces on periodic boundaries"));
2104 template <
typename MeshType>
2107 const MeshType & mesh,
2109 const unsigned int direction,
2115 static const int dim = MeshType::dimension;
2116 static const int space_dim = MeshType::space_dimension;
2126 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2127 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2129 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2130 cell != mesh.end(0);
2133 const typename MeshType::face_iterator face_1 =
2134 cell->face(2 * direction);
2135 const typename MeshType::face_iterator face_2 =
2136 cell->face(2 * direction + 1);
2138 if (face_1->at_boundary() && face_1->boundary_id() == b_id)
2140 const std::pair<typename MeshType::cell_iterator, unsigned int>
2141 pair1 = std::make_pair(cell, 2 * direction);
2142 pairs1.insert(pair1);
2145 if (face_2->at_boundary() && face_2->boundary_id() == b_id)
2147 const std::pair<typename MeshType::cell_iterator, unsigned int>
2148 pair2 = std::make_pair(cell, 2 * direction + 1);
2149 pairs2.insert(pair2);
2153 Assert(pairs1.size() == pairs2.size(),
2154 ExcMessage(
"Unmatched faces on periodic boundaries"));
2156 Assert(pairs1.size() > 0,
2157 ExcMessage(
"No new periodic face pairs have been found. "
2158 "Are you sure that you've selected the correct boundary "
2159 "id's and that the coarsest level mesh is colorized?"));
2162 const unsigned int size_old = matched_pairs.size();
2167 pairs1, pairs2, direction, matched_pairs, offset, matrix);
2171 const unsigned int size_new = matched_pairs.size();
2172 for (
unsigned int i = size_old; i < size_new; ++i)
2174 Assert(matched_pairs[i].orientation == 1,
2176 "Found a face match with non standard orientation. "
2177 "This function is only suitable for meshes with cells "
2178 "in default orientation"));
2185 template <
typename MeshType>
2188 const MeshType & mesh,
2191 const unsigned int direction,
2197 static const int dim = MeshType::dimension;
2198 static const int space_dim = MeshType::space_dimension;
2206 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2207 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2209 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2210 cell != mesh.end(0);
2213 for (
unsigned int i : cell->face_indices())
2215 const typename MeshType::face_iterator face = cell->face(i);
2216 if (face->at_boundary() && face->boundary_id() == b_id1)
2218 const std::pair<typename MeshType::cell_iterator, unsigned int>
2219 pair1 = std::make_pair(cell, i);
2220 pairs1.insert(pair1);
2223 if (face->at_boundary() && face->boundary_id() == b_id2)
2225 const std::pair<typename MeshType::cell_iterator, unsigned int>
2226 pair2 = std::make_pair(cell, i);
2227 pairs2.insert(pair2);
2238 if (!(((pairs1.size() > 0) &&
2241 *
>(&pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2242 ((pairs2.size() > 0) &&
2245 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2246 Assert(pairs1.size() == pairs2.size(),
2247 ExcMessage(
"Unmatched faces on periodic boundaries"));
2250 (pairs1.size() > 0 ||
2253 &mesh.begin()->get_triangulation()) !=
nullptr)),
2254 ExcMessage(
"No new periodic face pairs have been found. "
2255 "Are you sure that you've selected the correct boundary "
2256 "id's and that the coarsest level mesh is colorized?"));
2260 pairs1, pairs2, direction, matched_pairs, offset, matrix);
2274 template <
int spacedim>
2278 const unsigned int direction,
2288 if (matrix.m() == spacedim)
2289 for (
unsigned int i = 0; i < spacedim; ++i)
2290 for (
unsigned int j = 0; j < spacedim; ++j)
2291 distance(i) += matrix(i, j) * point1(j);
2295 distance += offset - point2;
2297 for (
unsigned int i = 0; i < spacedim; ++i)
2303 if (
std::abs(distance(i)) > 1.e-10)
2329 std::array<unsigned int, GeometryInfo<1>::vertices_per_face>;
2330 static inline std::bitset<3>
2342 std::array<unsigned int, GeometryInfo<2>::vertices_per_face>;
2343 static inline std::bitset<3>
2351 static const MATCH_T m_tff = {{0, 1}};
2352 if (matching == m_tff)
2354 static const MATCH_T m_ttf = {{1, 0}};
2355 if (matching == m_ttf)
2368 std::array<unsigned int, GeometryInfo<3>::vertices_per_face>;
2370 static inline std::bitset<3>
2378 static const MATCH_T m_tff = {{0, 1, 2, 3}};
2379 if (matching == m_tff)
2381 static const MATCH_T m_tft = {{1, 3, 0, 2}};
2382 if (matching == m_tft)
2384 static const MATCH_T m_ttf = {{3, 2, 1, 0}};
2385 if (matching == m_ttf)
2387 static const MATCH_T m_ttt = {{2, 0, 3, 1}};
2388 if (matching == m_ttt)
2390 static const MATCH_T m_fff = {{0, 2, 1, 3}};
2391 if (matching == m_fff)
2393 static const MATCH_T m_fft = {{2, 3, 0, 1}};
2394 if (matching == m_fft)
2396 static const MATCH_T m_ftf = {{3, 1, 2, 0}};
2397 if (matching == m_ftf)
2399 static const MATCH_T m_ftt = {{1, 0, 3, 2}};
2400 if (matching == m_ftt)
2411 template <
typename FaceIterator>
2414 std::bitset<3> & orientation,
2415 const FaceIterator & face1,
2416 const FaceIterator & face2,
2417 const unsigned int direction,
2421 Assert(matrix.m() == matrix.n(),
2422 ExcMessage(
"The supplied matrix must be a square matrix"));
2424 static const int dim = FaceIterator::AccessorType::dimension;
2428 std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
2432 std::set<unsigned int> face2_vertices;
2433 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2434 face2_vertices.insert(i);
2436 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2438 for (std::set<unsigned int>::iterator it = face2_vertices.begin();
2439 it != face2_vertices.end();
2449 face2_vertices.erase(it);
2456 if (face2_vertices.empty())
2459 return face2_vertices.empty();
2464 template <
typename FaceIterator>
2467 const FaceIterator & face1,
2468 const FaceIterator & face2,
2469 const unsigned int direction,
2474 std::bitset<3> dummy;
2483#include "grid_tools_dof_handlers.inst"
cell_iterator end() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() 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.
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
virtual bool preserves_vertex_locations() const =0
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
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()
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#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 > e(const Tensor< 2, dim, Number > &F)
static const unsigned int invalid_unsigned_int
::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
static double distance_to_unit_cell(const Point< dim > &p)
const ::Triangulation< dim, spacedim > & tria