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.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>
135 const MeshType<dim, spacedim> &mesh,
136 const
Point<spacedim> & p,
137 const std::vector<
bool> &marked_vertices)
140 if (mapping.preserves_vertex_locations() ==
true)
152 marked_vertices.size() == 0,
154 marked_vertices.size()));
164 marked_vertices.size() == 0 ||
165 std::equal(marked_vertices.begin(),
166 marked_vertices.end(),
168 [](
bool p,
bool q) { return !p || q; }),
170 "marked_vertices should be a subset of used vertices in the triangulation "
171 "but marked_vertices contains one or more vertices that are not used vertices!"));
174 if (marked_vertices.size())
177 if (marked_vertices[it->first] ==
false)
192 template <
class MeshType>
195 std::vector<typename MeshType::active_cell_iterator>
198 typename ::internal::ActiveCellIterator<MeshType::dimension,
199 MeshType::space_dimension,
203 const unsigned int vertex)
205 const int dim = MeshType::dimension;
206 const int spacedim = MeshType::space_dimension;
212 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
219 ActiveCellIterator<dim, spacedim, MeshType>::type>
222 typename ::internal::ActiveCellIterator<dim, spacedim, MeshType>::type
223 cell = mesh.begin_active(),
259 for (; cell != endc; ++cell)
261 for (
const unsigned int v : cell->vertex_indices())
262 if (cell->vertex_index(v) == vertex)
274 for (
const auto face :
275 cell->reference_cell().faces_for_given_vertex(v))
276 if (!cell->at_boundary(face) &&
277 cell->neighbor(face)->is_active())
299 for (
unsigned int e = 0;
e < cell->n_lines(); ++
e)
300 if (cell->line(
e)->has_children())
304 if (cell->line(
e)->child(0)->vertex_index(1) == vertex)
327 return std::vector<typename ::internal::
328 ActiveCellIterator<dim, spacedim, MeshType>::type>(
336 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
339 void find_active_cell_around_point_internal(
340 const MeshType<dim, spacedim> &mesh,
342 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>
344 std::set<typename MeshType<dim, spacedim>::active_cell_iterator>
348 typename ::internal::
349 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
352 typename ::internal::
353 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
358 using cell_iterator =
359 typename MeshType<dim, spacedim>::active_cell_iterator;
361 using cell_iterator = typename ::internal::
362 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
370 std::set<cell_iterator> adjacent_cells_new;
372 typename std::set<cell_iterator>::const_iterator cell =
376 for (; cell != endc; ++cell)
378 std::vector<cell_iterator> active_neighbors;
379 get_active_neighbors<MeshType<dim, spacedim>>(*cell,
381 for (
unsigned int i = 0; i < active_neighbors.size(); ++i)
382 if (searched_cells.find(active_neighbors[i]) ==
383 searched_cells.end())
384 adjacent_cells_new.insert(active_neighbors[i]);
388 adjacent_cells_new.end());
397 cell_iterator it = mesh.begin_active();
398 for (; it != mesh.end(); ++it)
399 if (searched_cells.find(it) == searched_cells.end())
410 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
414 typename MeshType<dim, spacedim>::active_cell_iterator
416 typename ::internal::
417 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
421 const std::vector<bool> &marked_vertices,
422 const double tolerance)
424 return find_active_cell_around_point<dim, MeshType, spacedim>(
435 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
439 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
441 std::pair<typename ::internal::
442 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
446 const MeshType<dim, spacedim> &mesh,
448 const std::vector<bool> &marked_vertices,
449 const double tolerance)
451 using active_cell_iterator = typename ::internal::
452 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
458 double best_distance = tolerance;
460 std::pair<active_cell_iterator, Point<dim>> best_cell;
463 best_cell.first = mesh.end();
467 std::vector<active_cell_iterator> adjacent_cells_tmp =
476 std::set<active_cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
477 adjacent_cells_tmp.end());
478 std::set<active_cell_iterator> searched_cells;
487 mesh.get_triangulation().n_active_cells();
489 unsigned int cells_searched = 0;
492 typename std::set<active_cell_iterator>::const_iterator
495 for (; cell != endc; ++cell)
497 if ((*cell)->is_artificial() ==
false)
500 if (marked_vertices.size() > 0)
502 bool any_vertex_marked =
false;
503 for (
const auto &v : (*cell)->vertex_indices())
505 if (marked_vertices[(*cell)->vertex_index(v)])
507 any_vertex_marked =
true;
511 if (!any_vertex_marked)
529 if ((dist < best_distance) ||
530 ((dist == best_distance) &&
531 ((*cell)->level() > best_level)))
534 best_distance = dist;
535 best_level = (*cell)->level();
536 best_cell = std::make_pair(*cell, p_cell);
569 find_active_cell_around_point_internal<dim, MeshType, spacedim>(
579 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
583 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
586 std::vector<std::pair<
587 typename ::internal::
588 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
592 const MeshType<dim, spacedim> &mesh,
594 const double tolerance,
595 const std::vector<bool> &marked_vertices)
598 mapping, mesh, p, marked_vertices, tolerance);
600 if (cell_and_point.first == mesh.end())
604 mapping, mesh, p, tolerance, cell_and_point);
609 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
613 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
616 std::vector<std::pair<
617 typename ::internal::
618 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
623 const MeshType<dim, spacedim> &mesh,
625 const double tolerance,
626 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
630 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
635 cells_and_points.push_back(first_cell);
639 const Point<dim> unit_point = cells_and_points.front().second;
640 const auto my_cell = cells_and_points.front().first;
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;
655 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
658 if (n_dirs_at_threshold == 1)
660 unsigned int neighbor_index =
661 2 * last_point_at_threshold +
662 (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
663 if (!my_cell->at_boundary(neighbor_index))
665 const auto neighbor_cell = my_cell->neighbor(neighbor_index);
667 if (neighbor_cell->is_active())
668 cells_to_add.push_back(neighbor_cell);
670 for (
const auto &child_cell : neighbor_cell->child_iterators())
672 if (child_cell->is_active())
673 cells_to_add.push_back(child_cell);
678 else if (n_dirs_at_threshold == dim)
680 unsigned int local_vertex_index = 0;
681 for (
unsigned int d = 0;
d < dim; ++
d)
682 local_vertex_index += (unit_point[
d] > 0.5 ? 1 : 0) <<
d;
683 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
685 mesh, my_cell->vertex_index(local_vertex_index));
686 for (
const auto &cell : cells)
689 cells_to_add.push_back(cell);
697 else if (n_dirs_at_threshold == 2)
700 unsigned int count_vertex_indices = 0;
702 for (
unsigned int d = 0;
d < dim; ++
d)
704 if (distance_to_center[
d] > 0.5 - tolerance)
708 unit_point[
d] > 0.5 ? 1 : 0;
709 count_vertex_indices++;
719 const unsigned int first_vertex =
722 for (
unsigned int d = 0;
d < 2; ++
d)
726 my_cell->vertex_index(first_vertex + (
d << free_direction)));
727 for (
const auto &cell : tentative_cells)
729 bool cell_not_yet_present =
true;
730 for (
const auto &other_cell : cells_to_add)
731 if (cell == other_cell)
733 cell_not_yet_present =
false;
736 if (cell_not_yet_present)
737 cells_to_add.push_back(cell);
742 const double original_distance_to_unit_cell =
744 for (
const auto &cell : cells_to_add)
752 original_distance_to_unit_cell + tolerance)
753 cells_and_points.emplace_back(cell, p_unit);
760 cells_and_points.begin(),
761 cells_and_points.end(),
762 [](
const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
764 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
767 return cells_and_points;
772 template <
class MeshType>
776 const MeshType &mesh,
777 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
780 std::vector<typename MeshType::active_cell_iterator> active_halo_layer;
781 std::vector<bool> locally_active_vertices_on_subdomain(
782 mesh.get_triangulation().n_vertices(),
false);
787 for (
const auto &cell : mesh.active_cell_iterators())
789 for (
const auto v : cell->vertex_indices())
790 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
true;
795 for (
const auto &cell : mesh.active_cell_iterators())
796 if (!predicate(cell))
797 for (
const auto v : cell->vertex_indices())
798 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
801 active_halo_layer.push_back(cell);
805 return active_halo_layer;
810 template <
class MeshType>
814 const MeshType &mesh,
815 const std::function<
bool(
const typename MeshType::cell_iterator &)>
817 const unsigned int level)
819 std::vector<typename MeshType::cell_iterator> level_halo_layer;
820 std::vector<bool> locally_active_vertices_on_level_subdomain(
821 mesh.get_triangulation().n_vertices(),
false);
826 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
827 cell != mesh.end(
level);
830 for (
const unsigned int v : cell->vertex_indices())
831 locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
837 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
838 cell != mesh.end(
level);
840 if (!predicate(cell))
841 for (
const unsigned int v : cell->vertex_indices())
842 if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
845 level_halo_layer.push_back(cell);
849 return level_halo_layer;
855 template <
class MeshType>
857 bool contains_locally_owned_cells(
858 const std::vector<typename MeshType::active_cell_iterator> &cells)
860 for (
typename std::vector<
861 typename MeshType::active_cell_iterator>::const_iterator it =
866 if ((*it)->is_locally_owned())
872 template <
class MeshType>
874 bool contains_artificial_cells(
875 const std::vector<typename MeshType::active_cell_iterator> &cells)
877 for (
typename std::vector<
878 typename MeshType::active_cell_iterator>::const_iterator it =
883 if ((*it)->is_artificial())
892 template <
class MeshType>
898 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
901 const std::vector<typename MeshType::active_cell_iterator>
906 Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) ==
false,
907 ExcMessage(
"Halo layer contains locally owned cells"));
908 Assert(contains_artificial_cells<MeshType>(active_halo_layer) ==
false,
909 ExcMessage(
"Halo layer contains artificial cells"));
911 return active_halo_layer;
916 template <
class MeshType>
920 const MeshType &mesh,
921 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
923 const double layer_thickness)
925 std::vector<typename MeshType::active_cell_iterator>
926 subdomain_boundary_cells, active_cell_layer_within_distance;
927 std::vector<bool> vertices_outside_subdomain(
928 mesh.get_triangulation().n_vertices(),
false);
930 const unsigned int spacedim = MeshType::space_dimension;
932 unsigned int n_non_predicate_cells = 0;
940 for (
const auto &cell : mesh.active_cell_iterators())
941 if (!predicate(cell))
943 for (
const unsigned int v : cell->vertex_indices())
944 vertices_outside_subdomain[cell->vertex_index(v)] =
true;
945 n_non_predicate_cells++;
952 if (n_non_predicate_cells == 0 ||
953 n_non_predicate_cells == mesh.get_triangulation().n_active_cells())
954 return std::vector<typename MeshType::active_cell_iterator>();
958 for (
const auto &cell : mesh.active_cell_iterators())
961 for (
const unsigned int v : cell->vertex_indices())
962 if (vertices_outside_subdomain[cell->vertex_index(v)] ==
true)
964 subdomain_boundary_cells.push_back(cell);
977 for (
unsigned int d = 0;
d < spacedim; ++
d)
979 bounding_box.first[
d] -= (layer_thickness + DOUBLE_EPSILON);
980 bounding_box.second[
d] += (layer_thickness + DOUBLE_EPSILON);
983 std::vector<Point<spacedim>>
984 subdomain_boundary_cells_centers;
987 subdomain_boundary_cells_radii;
989 subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size());
990 subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size());
992 for (
typename std::vector<typename MeshType::active_cell_iterator>::
993 const_iterator subdomain_boundary_cell_iterator =
994 subdomain_boundary_cells.begin();
995 subdomain_boundary_cell_iterator != subdomain_boundary_cells.end();
996 ++subdomain_boundary_cell_iterator)
998 const std::pair<Point<spacedim>,
double>
999 &subdomain_boundary_cell_enclosing_ball =
1000 (*subdomain_boundary_cell_iterator)->enclosing_ball();
1002 subdomain_boundary_cells_centers.push_back(
1003 subdomain_boundary_cell_enclosing_ball.first);
1004 subdomain_boundary_cells_radii.push_back(
1005 subdomain_boundary_cell_enclosing_ball.second);
1007 AssertThrow(subdomain_boundary_cells_radii.size() ==
1008 subdomain_boundary_cells_centers.size(),
1017 for (
const auto &cell : mesh.active_cell_iterators())
1020 if (predicate(cell))
1023 const std::pair<Point<spacedim>,
double> &cell_enclosing_ball =
1024 cell->enclosing_ball();
1027 cell_enclosing_ball.first;
1028 const double cell_enclosing_ball_radius = cell_enclosing_ball.second;
1030 bool cell_inside =
true;
1032 for (
unsigned int d = 0;
d < spacedim; ++
d)
1034 (cell_enclosing_ball_center[
d] + cell_enclosing_ball_radius >
1035 bounding_box.first[
d]) &&
1036 (cell_enclosing_ball_center[
d] - cell_enclosing_ball_radius <
1037 bounding_box.second[
d]);
1043 for (
unsigned int i = 0; i < subdomain_boundary_cells_radii.size();
1046 subdomain_boundary_cells_centers[i]) <
1047 Utilities::fixed_power<2>(cell_enclosing_ball_radius +
1048 subdomain_boundary_cells_radii[i] +
1049 layer_thickness + DOUBLE_EPSILON))
1051 active_cell_layer_within_distance.push_back(cell);
1056 return active_cell_layer_within_distance;
1061 template <
class MeshType>
1071 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1072 predicate(locally_owned_cell_predicate);
1074 const std::vector<typename MeshType::active_cell_iterator>
1075 ghost_cell_layer_within_distance =
1083 contains_locally_owned_cells<MeshType>(
1084 ghost_cell_layer_within_distance) ==
false,
1086 "Ghost cells within layer_thickness contains locally owned cells."));
1088 contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) ==
1091 "Ghost cells within layer_thickness contains artificial cells. "
1092 "The function compute_ghost_cell_layer_within_distance "
1093 "is probably called while using parallel::distributed::Triangulation. "
1094 "In such case please refer to the description of this function."));
1096 return ghost_cell_layer_within_distance;
1101 template <
class MeshType>
1107 const std::function<
bool(
1108 const typename MeshType::
1109 active_cell_iterator &)>
1112 std::vector<bool> locally_active_vertices_on_subdomain(
1113 mesh.get_triangulation().n_vertices(),
false);
1115 const unsigned int spacedim = MeshType::space_dimension;
1122 for (
const auto &cell : mesh.active_cell_iterators())
1123 if (predicate(cell))
1125 minp = cell->center();
1126 maxp = cell->center();
1132 for (
const auto &cell : mesh.active_cell_iterators())
1133 if (predicate(cell))
1134 for (
const unsigned int v : cell->vertex_indices())
1135 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
1138 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
1140 for (
unsigned int d = 0;
d < spacedim; ++
d)
1147 return std::make_pair(minp, maxp);
1152 template <
typename MeshType>
1154 std::list<std::pair<
1155 typename MeshType::cell_iterator,
1162 ExcMessage(
"The two meshes must be represent triangulations that "
1163 "have the same coarse meshes"));
1170 bool remove_ghost_cells =
false;
1171 #ifdef DEAL_II_WITH_MPI
1173 constexpr
int dim = MeshType::dimension;
1174 constexpr
int spacedim = MeshType::space_dimension;
1176 *
>(&mesh_1.get_triangulation()) !=
nullptr ||
1178 *
>(&mesh_2.get_triangulation()) !=
nullptr)
1180 Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(),
1181 ExcMessage(
"This function can only be used with meshes "
1182 "corresponding to distributed Triangulations when "
1183 "both Triangulations are equal."));
1184 remove_ghost_cells =
true;
1196 using CellList = std::list<std::pair<
typename MeshType::cell_iterator,
1197 typename MeshType::cell_iterator>>;
1201 typename MeshType::cell_iterator cell_1 = mesh_1.begin(0),
1202 cell_2 = mesh_2.begin(0);
1203 for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
1204 cell_list.emplace_back(cell_1, cell_2);
1207 typename CellList::iterator cell_pair = cell_list.begin();
1208 while (cell_pair != cell_list.end())
1212 if (cell_pair->first->has_children() &&
1213 cell_pair->second->has_children())
1215 Assert(cell_pair->first->refinement_case() ==
1216 cell_pair->second->refinement_case(),
1218 for (
unsigned int c = 0; c < cell_pair->first->n_children(); ++c)
1219 cell_list.emplace_back(cell_pair->first->child(c),
1220 cell_pair->second->child(c));
1225 const auto previous_cell_pair = cell_pair;
1227 cell_list.erase(previous_cell_pair);
1232 if (remove_ghost_cells &&
1233 ((cell_pair->first->is_active() &&
1234 !cell_pair->first->is_locally_owned()) ||
1235 (cell_pair->second->is_active() &&
1236 !cell_pair->second->is_locally_owned())))
1239 const auto previous_cell_pair = cell_pair;
1241 cell_list.erase(previous_cell_pair);
1250 for (cell_pair = cell_list.begin(); cell_pair != cell_list.end();
1252 Assert(cell_pair->first->is_active() || cell_pair->second->is_active() ||
1253 (cell_pair->first->refinement_case() !=
1254 cell_pair->second->refinement_case()),
1262 template <
int dim,
int spacedim>
1279 endc = mesh_1.
end(0);
1280 for (; cell_1 != endc; ++cell_1, ++cell_2)
1282 if (cell_1->
n_vertices() != cell_2->n_vertices())
1284 for (
const unsigned int v : cell_1->vertex_indices())
1285 if (cell_1->vertex(v) != cell_2->vertex(v))
1298 template <
typename MeshType>
1303 mesh_2.get_triangulation());
1308 template <
int dim,
int spacedim>
1309 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1315 const double tolerance)
1319 ExcMessage(
"Mapping collection needs to have either size 1 "
1320 "or size equal to the number of elements in "
1321 "the FECollection."));
1323 using cell_iterator =
1326 std::pair<cell_iterator, Point<dim>> best_cell;
1330 if (mapping.
size() == 1)
1332 const std::vector<bool> marked_vertices = {};
1334 mapping[0], mesh, p, marked_vertices, tolerance);
1341 double best_distance = tolerance;
1342 int best_level = -1;
1349 std::vector<cell_iterator> adjacent_cells_tmp =
1357 std::set<cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
1358 adjacent_cells_tmp.end());
1359 std::set<cell_iterator> searched_cells;
1369 unsigned int cells_searched = 0;
1370 while (!found && cells_searched <
n_cells)
1372 typename std::set<cell_iterator>::const_iterator
1375 for (; cell != endc; ++cell)
1380 mapping[(*cell)->active_fe_index()]
1381 .transform_real_to_unit_cell(*cell, p);
1393 if (dist < best_distance || (dist == best_distance &&
1394 (*cell)->level() > best_level))
1397 best_distance = dist;
1398 best_level = (*cell)->level();
1399 best_cell = std::make_pair(*cell, p_cell);
1425 if (!found && cells_searched <
n_cells)
1427 find_active_cell_around_point_internal<dim,
1439 template <
class MeshType>
1442 const typename MeshType::active_cell_iterator &cell)
1444 Assert(cell->is_locally_owned(),
1445 ExcMessage(
"This function only makes sense if the cell for "
1446 "which you are asking for a patch, is locally "
1449 std::vector<typename MeshType::active_cell_iterator> patch;
1450 patch.push_back(cell);
1451 for (
const unsigned int face_number : cell->face_indices())
1452 if (cell->face(face_number)->at_boundary() ==
false)
1454 if (cell->neighbor(face_number)->has_children() ==
false)
1455 patch.push_back(cell->neighbor(face_number));
1460 if (MeshType::dimension > 1)
1462 for (
unsigned int subface = 0;
1463 subface < cell->face(face_number)->n_children();
1466 cell->neighbor_child_on_subface(face_number, subface));
1472 typename MeshType::cell_iterator neighbor =
1473 cell->neighbor(face_number);
1474 while (neighbor->has_children())
1475 neighbor = neighbor->child(1 - face_number);
1477 Assert(neighbor->neighbor(1 - face_number) == cell,
1479 patch.push_back(neighbor);
1487 template <
class Container>
1488 std::vector<typename Container::cell_iterator>
1490 const std::vector<typename Container::active_cell_iterator> &patch)
1494 "Vector containing patch cells should not be an empty vector!"));
1498 int min_level = patch[0]->level();
1500 for (
unsigned int i = 0; i < patch.size(); ++i)
1502 std::set<typename Container::cell_iterator> uniform_cells;
1503 typename std::vector<
1504 typename Container::active_cell_iterator>::const_iterator patch_cell;
1506 for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell)
1511 if ((*patch_cell)->level() == min_level)
1512 uniform_cells.insert(*patch_cell);
1519 typename Container::cell_iterator parent = *patch_cell;
1521 while (parent->level() > min_level)
1522 parent = parent->parent();
1523 uniform_cells.insert(parent);
1527 return std::vector<typename Container::cell_iterator>(uniform_cells.begin(),
1528 uniform_cells.end());
1533 template <
class Container>
1536 const std::vector<typename Container::active_cell_iterator> &patch,
1538 &local_triangulation,
1541 Container::space_dimension>::active_cell_iterator,
1542 typename Container::active_cell_iterator> &patch_to_global_tria_map)
1545 const std::vector<typename Container::cell_iterator> uniform_cells =
1546 get_cells_at_coarsest_common_level<Container>(patch);
1548 local_triangulation.
clear();
1549 std::vector<Point<Container::space_dimension>>
vertices;
1550 const unsigned int n_uniform_cells = uniform_cells.size();
1551 std::vector<CellData<Container::dimension>> cells(n_uniform_cells);
1554 typename std::vector<typename Container::cell_iterator>::const_iterator
1556 for (uniform_cell = uniform_cells.begin();
1557 uniform_cell != uniform_cells.end();
1560 for (
const unsigned int v : (*uniform_cell)->vertex_indices())
1563 (*uniform_cell)->vertex(v);
1564 bool repeat_vertex =
false;
1566 for (
unsigned int m = 0; m < i; ++m)
1570 repeat_vertex =
true;
1571 cells[k].vertices[v] = m;
1575 if (repeat_vertex ==
false)
1578 cells[k].vertices[v] = i;
1589 unsigned int index = 0;
1592 Container::space_dimension>::cell_iterator,
1593 typename Container::cell_iterator>
1594 patch_to_global_tria_map_tmp;
1596 Container::space_dimension>::cell_iterator
1597 coarse_cell = local_triangulation.
begin();
1598 coarse_cell != local_triangulation.
end();
1599 ++coarse_cell, ++index)
1601 patch_to_global_tria_map_tmp.insert(
1602 std::make_pair(coarse_cell, uniform_cells[index]));
1606 Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <=
1607 1
e-15 * coarse_cell->diameter(),
1610 bool refinement_necessary;
1615 refinement_necessary =
false;
1616 for (
const auto &active_tria_cell :
1619 if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
1621 active_tria_cell->set_refine_flag();
1622 refinement_necessary =
true;
1625 for (
unsigned int i = 0; i < patch.size(); ++i)
1630 if (patch_to_global_tria_map_tmp[active_tria_cell] ==
1635 for (
const unsigned int v :
1636 active_tria_cell->vertex_indices())
1637 active_tria_cell->vertex(v) = patch[i]->vertex(v);
1639 Assert(active_tria_cell->center().distance(
1640 patch_to_global_tria_map_tmp[active_tria_cell]
1642 1
e-15 * active_tria_cell->diameter(),
1645 active_tria_cell->set_user_flag();
1651 if (refinement_necessary)
1656 Container::dimension,
1657 Container::space_dimension>::cell_iterator cell =
1658 local_triangulation.
begin();
1659 cell != local_triangulation.
end();
1662 if (patch_to_global_tria_map_tmp.find(cell) !=
1663 patch_to_global_tria_map_tmp.end())
1665 if (cell->has_children())
1672 for (
unsigned int c = 0; c < cell->n_children(); ++c)
1674 if (patch_to_global_tria_map_tmp.find(cell->child(
1675 c)) == patch_to_global_tria_map_tmp.end())
1677 patch_to_global_tria_map_tmp.insert(
1680 patch_to_global_tria_map_tmp[cell]->child(
1702 patch_to_global_tria_map_tmp.erase(cell);
1708 while (refinement_necessary);
1714 Container::space_dimension>::cell_iterator
1715 cell = local_triangulation.
begin();
1716 cell != local_triangulation.
end();
1719 if (cell->user_flag_set())
1721 Assert(patch_to_global_tria_map_tmp.find(cell) !=
1722 patch_to_global_tria_map_tmp.end(),
1725 Assert(cell->center().distance(
1726 patch_to_global_tria_map_tmp[cell]->center()) <=
1727 1
e-15 * cell->diameter(),
1735 Container::space_dimension>::cell_iterator,
1736 typename Container::cell_iterator>::iterator
1737 map_tmp_it = patch_to_global_tria_map_tmp.
begin(),
1738 map_tmp_end = patch_to_global_tria_map_tmp.end();
1742 for (; map_tmp_it != map_tmp_end; ++map_tmp_it)
1743 patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
1748 template <
int dim,
int spacedim>
1751 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1765 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1766 dof_to_set_of_cells_map;
1768 std::vector<types::global_dof_index> local_dof_indices;
1769 std::vector<types::global_dof_index> local_face_dof_indices;
1770 std::vector<types::global_dof_index> local_line_dof_indices;
1774 std::vector<bool> user_flags;
1779 std::map<typename DoFHandler<dim, spacedim>::active_line_iterator,
1781 lines_to_parent_lines_map;
1788 .clear_user_flags();
1793 endc = dof_handler.
end();
1794 for (; cell != endc; ++cell)
1800 if (cell->is_artificial() ==
false)
1802 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
1803 if (cell->line(
l)->has_children())
1804 for (
unsigned int c = 0; c < cell->line(
l)->n_children();
1807 lines_to_parent_lines_map[cell->line(
l)->child(c)] =
1811 cell->line(
l)->child(c)->set_user_flag();
1826 endc = dof_handler.
end();
1827 for (; cell != endc; ++cell)
1832 if (cell->is_artificial() ==
false)
1834 const unsigned int n_dofs_per_cell =
1836 local_dof_indices.resize(n_dofs_per_cell);
1840 cell->get_dof_indices(local_dof_indices);
1841 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i)
1842 dof_to_set_of_cells_map[local_dof_indices[i]].
insert(cell);
1852 for (
const unsigned int f : cell->face_indices())
1854 if (cell->face(f)->has_children())
1856 for (
unsigned int c = 0; c < cell->face(f)->n_children();
1871 Assert(cell->face(f)->child(c)->has_children() ==
false,
1874 const unsigned int n_dofs_per_face =
1876 local_face_dof_indices.resize(n_dofs_per_face);
1878 cell->face(f)->child(c)->get_dof_indices(
1879 local_face_dof_indices);
1880 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1881 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1885 else if ((cell->face(f)->at_boundary() ==
false) &&
1886 (cell->neighbor_is_coarser(f)))
1903 std::pair<unsigned int, unsigned int>
1904 neighbor_face_no_subface_no =
1905 cell->neighbor_of_coarser_neighbor(f);
1906 unsigned int face_no = neighbor_face_no_subface_no.first;
1907 unsigned int subface = neighbor_face_no_subface_no.second;
1909 const unsigned int n_dofs_per_face =
1911 local_face_dof_indices.resize(n_dofs_per_face);
1913 cell->neighbor(f)->face(face_no)->get_dof_indices(
1914 local_face_dof_indices);
1915 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1916 dof_to_set_of_cells_map[local_face_dof_indices[i]].
insert(
1921 for (
unsigned int c = 0;
1922 c < cell->neighbor(f)->face(face_no)->n_children();
1928 const unsigned int n_dofs_per_face =
1930 local_face_dof_indices.resize(n_dofs_per_face);
1935 ->has_children() ==
false,
1940 ->get_dof_indices(local_face_dof_indices);
1941 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1942 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1959 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
1961 if (cell->line(
l)->has_children())
1963 for (
unsigned int c = 0;
1964 c < cell->line(
l)->n_children();
1967 Assert(cell->line(
l)->child(c)->has_children() ==
1973 const unsigned int n_dofs_per_line =
1976 local_line_dof_indices.resize(n_dofs_per_line);
1978 cell->line(
l)->child(c)->get_dof_indices(
1979 local_line_dof_indices);
1980 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1981 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1989 else if (cell->line(
l)->user_flag_set() ==
true)
1993 lines_to_parent_lines_map[cell->line(
l)];
1998 const unsigned int n_dofs_per_line =
2001 local_line_dof_indices.resize(n_dofs_per_line);
2003 parent_line->get_dof_indices(local_line_dof_indices);
2004 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2005 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2008 for (
unsigned int c = 0; c < parent_line->n_children();
2011 Assert(parent_line->child(c)->has_children() ==
2015 const unsigned int n_dofs_per_line =
2018 local_line_dof_indices.resize(n_dofs_per_line);
2020 parent_line->child(c)->get_dof_indices(
2021 local_line_dof_indices);
2022 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2023 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2040 .load_user_flags(user_flags);
2047 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2048 dof_to_cell_patches;
2052 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>::
2053 iterator it = dof_to_set_of_cells_map.begin(),
2054 it_end = dof_to_set_of_cells_map.end();
2055 for (; it != it_end; ++it)
2056 dof_to_cell_patches[it->first].assign(it->second.begin(),
2059 return dof_to_cell_patches;
2065 template <
typename CellIterator>
2068 std::set<std::pair<CellIterator, unsigned int>> &pairs1,
2071 const unsigned int direction,
2073 const ::Tensor<1, CellIterator::AccessorType::space_dimension>
2077 static const int space_dim = CellIterator::AccessorType::space_dimension;
2083 constexpr
int dim = CellIterator::AccessorType::dimension;
2084 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 Assert(pairs1.size() == pairs2.size(),
2098 ExcMessage(
"Unmatched faces on periodic boundaries"));
2102 unsigned int n_matches = 0;
2105 std::bitset<3> orientation;
2106 using PairIterator =
2107 typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator;
2108 for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
2110 for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
2112 const CellIterator cell1 = it1->first;
2113 const CellIterator cell2 = it2->first;
2114 const unsigned int face_idx1 = it1->second;
2115 const unsigned int face_idx2 = it2->second;
2117 cell1->face(face_idx1),
2118 cell2->face(face_idx2),
2127 {cell1, cell2}, {face_idx1, face_idx2}, orientation,
matrix};
2128 matched_pairs.push_back(matched_face);
2142 constexpr
int dim = CellIterator::AccessorType::dimension;
2143 constexpr
int spacedim = CellIterator::AccessorType::space_dimension;
2144 if (!(((pairs1.size() > 0) &&
2145 (
dynamic_cast<const parallel::fullydistributed::
2146 Triangulation<dim, spacedim> *
>(
2147 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2148 ((pairs2.size() > 0) &&
2151 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2152 AssertThrow(n_matches == pairs1.size() && pairs2.size() == 0,
2153 ExcMessage(
"Unmatched faces on periodic boundaries"));
2159 template <
typename MeshType>
2162 const MeshType & mesh,
2164 const
unsigned int direction,
2167 const
Tensor<1, MeshType::space_dimension> &offset,
2170 static const int dim = MeshType::dimension;
2171 static const int space_dim = MeshType::space_dimension;
2181 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2182 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2184 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2185 cell != mesh.end(0);
2188 const typename MeshType::face_iterator face_1 =
2189 cell->face(2 * direction);
2190 const typename MeshType::face_iterator face_2 =
2191 cell->face(2 * direction + 1);
2193 if (face_1->at_boundary() && face_1->boundary_id() == b_id)
2195 const std::pair<typename MeshType::cell_iterator, unsigned int>
2196 pair1 = std::make_pair(cell, 2 * direction);
2197 pairs1.insert(pair1);
2200 if (face_2->at_boundary() && face_2->boundary_id() == b_id)
2202 const std::pair<typename MeshType::cell_iterator, unsigned int>
2203 pair2 = std::make_pair(cell, 2 * direction + 1);
2204 pairs2.insert(pair2);
2208 Assert(pairs1.size() == pairs2.size(),
2209 ExcMessage(
"Unmatched faces on periodic boundaries"));
2211 Assert(pairs1.size() > 0,
2212 ExcMessage(
"No new periodic face pairs have been found. "
2213 "Are you sure that you've selected the correct boundary "
2214 "id's and that the coarsest level mesh is colorized?"));
2217 const unsigned int size_old = matched_pairs.size();
2222 pairs1, pairs2, direction, matched_pairs, offset,
matrix);
2226 const unsigned int size_new = matched_pairs.size();
2227 for (
unsigned int i = size_old; i < size_new; ++i)
2229 Assert(matched_pairs[i].orientation == 1,
2231 "Found a face match with non standard orientation. "
2232 "This function is only suitable for meshes with cells "
2233 "in default orientation"));
2240 template <
typename MeshType>
2243 const MeshType & mesh,
2246 const
unsigned int direction,
2249 const
Tensor<1, MeshType::space_dimension> &offset,
2252 static const int dim = MeshType::dimension;
2253 static const int space_dim = MeshType::space_dimension;
2261 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2262 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2264 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2265 cell != mesh.end(0);
2268 for (
unsigned int i : cell->face_indices())
2270 const typename MeshType::face_iterator face = cell->face(i);
2271 if (face->at_boundary() && face->boundary_id() == b_id1)
2273 const std::pair<typename MeshType::cell_iterator, unsigned int>
2274 pair1 = std::make_pair(cell, i);
2275 pairs1.insert(pair1);
2278 if (face->at_boundary() && face->boundary_id() == b_id2)
2280 const std::pair<typename MeshType::cell_iterator, unsigned int>
2281 pair2 = std::make_pair(cell, i);
2282 pairs2.insert(pair2);
2293 if (!(((pairs1.size() > 0) &&
2296 *
>(&pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2297 ((pairs2.size() > 0) &&
2300 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2301 Assert(pairs1.size() == pairs2.size(),
2302 ExcMessage(
"Unmatched faces on periodic boundaries"));
2305 (pairs1.size() > 0 ||
2308 &mesh.begin()->get_triangulation()) !=
nullptr)),
2309 ExcMessage(
"No new periodic face pairs have been found. "
2310 "Are you sure that you've selected the correct boundary "
2311 "id's and that the coarsest level mesh is colorized?"));
2315 pairs1, pairs2, direction, matched_pairs, offset,
matrix);
2329 template <
int spacedim>
2333 const unsigned int direction,
2343 if (
matrix.m() == spacedim)
2344 for (
unsigned int i = 0; i < spacedim; ++i)
2345 for (
unsigned int j = 0; j < spacedim; ++j)
2346 distance(i) +=
matrix(i, j) * point1(j);
2350 distance += offset - point2;
2352 for (
unsigned int i = 0; i < spacedim; ++i)
2358 if (std::abs(distance(i)) > 1.e-10)
2384 std::array<unsigned int, GeometryInfo<1>::vertices_per_face>;
2385 static inline std::bitset<3>
2397 std::array<unsigned int, GeometryInfo<2>::vertices_per_face>;
2398 static inline std::bitset<3>
2406 static const MATCH_T m_tff = {{0, 1}};
2407 if (matching == m_tff)
2409 static const MATCH_T m_ttf = {{1, 0}};
2410 if (matching == m_ttf)
2423 std::array<unsigned int, GeometryInfo<3>::vertices_per_face>;
2425 static inline std::bitset<3>
2433 static const MATCH_T m_tff = {{0, 1, 2, 3}};
2434 if (matching == m_tff)
2436 static const MATCH_T m_tft = {{1, 3, 0, 2}};
2437 if (matching == m_tft)
2439 static const MATCH_T m_ttf = {{3, 2, 1, 0}};
2440 if (matching == m_ttf)
2442 static const MATCH_T m_ttt = {{2, 0, 3, 1}};
2443 if (matching == m_ttt)
2445 static const MATCH_T m_fff = {{0, 2, 1, 3}};
2446 if (matching == m_fff)
2448 static const MATCH_T m_fft = {{2, 3, 0, 1}};
2449 if (matching == m_fft)
2451 static const MATCH_T m_ftf = {{3, 1, 2, 0}};
2452 if (matching == m_ftf)
2454 static const MATCH_T m_ftt = {{1, 0, 3, 2}};
2455 if (matching == m_ftt)
2466 template <
typename FaceIterator>
2469 std::bitset<3> & orientation,
2470 const FaceIterator & face1,
2471 const FaceIterator & face2,
2472 const unsigned int direction,
2477 ExcMessage(
"The supplied matrix must be a square matrix"));
2479 static const int dim = FaceIterator::AccessorType::dimension;
2483 std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
2487 std::set<unsigned int> face2_vertices;
2488 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2489 face2_vertices.insert(i);
2491 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2493 for (std::set<unsigned int>::iterator it = face2_vertices.begin();
2494 it != face2_vertices.end();
2504 face2_vertices.erase(it);
2511 if (face2_vertices.empty())
2514 return face2_vertices.empty();
2519 template <
typename FaceIterator>
2522 const FaceIterator & face1,
2523 const FaceIterator & face2,
2524 const unsigned int direction,
2529 std::bitset<3> dummy;
2538 #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
static double distance_to_unit_cell(const Point< dim > &p)
const ::Triangulation< dim, spacedim > & tria