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 <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
196 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
199 typename ::internal::
200 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
203 const unsigned int vertex)
209 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
216 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
219 typename ::internal::
220 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
221 cell = mesh.begin_active(),
257 for (; cell != endc; ++cell)
259 for (
const unsigned int v : cell->vertex_indices())
260 if (cell->vertex_index(v) == vertex)
272 for (
const auto face :
273 cell->reference_cell().faces_for_given_vertex(v))
274 if (!cell->at_boundary(face) &&
275 cell->neighbor(face)->is_active())
297 for (
unsigned int e = 0;
e < cell->n_lines(); ++
e)
298 if (cell->line(
e)->has_children())
302 if (cell->line(
e)->child(0)->vertex_index(1) == vertex)
326 typename ::internal::
327 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::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;
371 typename std::set<cell_iterator>::const_iterator cell =
375 for (; cell != endc; ++cell)
377 std::vector<cell_iterator> active_neighbors;
378 get_active_neighbors<MeshType<dim, spacedim>>(*cell,
380 for (
unsigned int i = 0; i < active_neighbors.size(); ++i)
381 if (searched_cells.find(active_neighbors[i]) ==
382 searched_cells.end())
383 adjacent_cells_new.insert(active_neighbors[i]);
387 adjacent_cells_new.end());
396 cell_iterator it = mesh.begin_active();
397 for (; it != mesh.end(); ++it)
398 if (searched_cells.find(it) == searched_cells.end())
409 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
413 typename MeshType<dim, spacedim>::active_cell_iterator
415 typename ::internal::
416 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
420 const std::vector<bool> &marked_vertices,
421 const double tolerance)
423 return find_active_cell_around_point<dim, MeshType, spacedim>(
434 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
438 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
440 std::pair<typename ::internal::
441 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
445 const MeshType<dim, spacedim> &mesh,
447 const std::vector<bool> &marked_vertices,
448 const double tolerance)
450 using active_cell_iterator = typename ::internal::
451 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
457 double best_distance = tolerance;
459 std::pair<active_cell_iterator, Point<dim>> best_cell;
462 best_cell.first = mesh.end();
466 std::vector<active_cell_iterator> adjacent_cells_tmp =
475 std::set<active_cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
476 adjacent_cells_tmp.end());
477 std::set<active_cell_iterator> searched_cells;
486 mesh.get_triangulation().n_active_cells();
488 unsigned int cells_searched = 0;
491 typename std::set<active_cell_iterator>::const_iterator
494 for (; cell != endc; ++cell)
496 if ((*cell)->is_artificial() ==
false)
512 if ((dist < best_distance) ||
513 ((dist == best_distance) &&
514 ((*cell)->level() > best_level)))
517 best_distance = dist;
518 best_level = (*cell)->level();
519 best_cell = std::make_pair(*cell, p_cell);
546 if (marked_vertices.size() > 0)
558 find_active_cell_around_point_internal<dim, MeshType, spacedim>(
568 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
572 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
575 std::vector<std::pair<
576 typename ::internal::
577 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
581 const MeshType<dim, spacedim> &mesh,
583 const double tolerance,
584 const std::vector<bool> &marked_vertices)
587 mapping, mesh, p, marked_vertices, tolerance);
589 if (cell_and_point.first == mesh.end())
593 mapping, mesh, p, tolerance, cell_and_point);
598 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
602 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
605 std::vector<std::pair<
606 typename ::internal::
607 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
612 const MeshType<dim, spacedim> &mesh,
614 const double tolerance,
615 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
619 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
624 cells_and_points.push_back(first_cell);
628 const Point<dim> unit_point = cells_and_points.front().second;
629 const auto my_cell = cells_and_points.front().first;
632 unsigned int n_dirs_at_threshold = 0;
634 for (
unsigned int d = 0;
d < dim; ++
d)
636 distance_to_center[
d] = std::abs(unit_point[
d] - 0.5);
637 if (distance_to_center[
d] > 0.5 - tolerance)
639 ++n_dirs_at_threshold;
640 last_point_at_threshold =
d;
644 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
647 if (n_dirs_at_threshold == 1)
649 unsigned int neighbor_index =
650 2 * last_point_at_threshold +
651 (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
652 if (!my_cell->at_boundary(neighbor_index))
654 const auto neighbor_cell = my_cell->neighbor(neighbor_index);
656 if (neighbor_cell->is_active())
657 cells_to_add.push_back(neighbor_cell);
659 for (
const auto &child_cell : neighbor_cell->child_iterators())
661 if (child_cell->is_active())
662 cells_to_add.push_back(child_cell);
667 else if (n_dirs_at_threshold == dim)
669 unsigned int local_vertex_index = 0;
670 for (
unsigned int d = 0;
d < dim; ++
d)
671 local_vertex_index += (unit_point[
d] > 0.5 ? 1 : 0) <<
d;
672 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
674 mesh, my_cell->vertex_index(local_vertex_index));
675 for (
const auto &cell : cells)
678 cells_to_add.push_back(cell);
686 else if (n_dirs_at_threshold == 2)
689 unsigned int count_vertex_indices = 0;
691 for (
unsigned int d = 0;
d < dim; ++
d)
693 if (distance_to_center[
d] > 0.5 - tolerance)
697 unit_point[
d] > 0.5 ? 1 : 0;
698 count_vertex_indices++;
708 const unsigned int first_vertex =
711 for (
unsigned int d = 0;
d < 2; ++
d)
715 my_cell->vertex_index(first_vertex + (
d << free_direction)));
716 for (
const auto &cell : tentative_cells)
718 bool cell_not_yet_present =
true;
719 for (
const auto &other_cell : cells_to_add)
720 if (cell == other_cell)
722 cell_not_yet_present =
false;
725 if (cell_not_yet_present)
726 cells_to_add.push_back(cell);
731 const double original_distance_to_unit_cell =
733 for (
const auto &cell : cells_to_add)
741 original_distance_to_unit_cell + tolerance)
742 cells_and_points.emplace_back(cell, p_unit);
749 cells_and_points.begin(),
750 cells_and_points.end(),
751 [](
const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
753 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
756 return cells_and_points;
761 template <
class MeshType>
765 const MeshType &mesh,
766 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
769 std::vector<typename MeshType::active_cell_iterator> active_halo_layer;
770 std::vector<bool> locally_active_vertices_on_subdomain(
771 mesh.get_triangulation().n_vertices(),
false);
776 for (
const auto &cell : mesh.active_cell_iterators())
778 for (
const auto v : cell->vertex_indices())
779 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
true;
784 for (
const auto &cell : mesh.active_cell_iterators())
785 if (!predicate(cell))
786 for (
const auto v : cell->vertex_indices())
787 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
790 active_halo_layer.push_back(cell);
794 return active_halo_layer;
799 template <
class MeshType>
803 const MeshType &mesh,
804 const std::function<
bool(
const typename MeshType::cell_iterator &)>
806 const unsigned int level)
808 std::vector<typename MeshType::cell_iterator> level_halo_layer;
809 std::vector<bool> locally_active_vertices_on_level_subdomain(
810 mesh.get_triangulation().n_vertices(),
false);
815 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
816 cell != mesh.end(
level);
819 for (
const unsigned int v : cell->vertex_indices())
820 locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
826 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
827 cell != mesh.end(
level);
829 if (!predicate(cell))
830 for (
const unsigned int v : cell->vertex_indices())
831 if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
834 level_halo_layer.push_back(cell);
838 return level_halo_layer;
844 template <
class MeshType>
846 bool contains_locally_owned_cells(
847 const std::vector<typename MeshType::active_cell_iterator> &cells)
849 for (
typename std::vector<
850 typename MeshType::active_cell_iterator>::const_iterator it =
855 if ((*it)->is_locally_owned())
861 template <
class MeshType>
863 bool contains_artificial_cells(
864 const std::vector<typename MeshType::active_cell_iterator> &cells)
866 for (
typename std::vector<
867 typename MeshType::active_cell_iterator>::const_iterator it =
872 if ((*it)->is_artificial())
881 template <
class MeshType>
887 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
890 const std::vector<typename MeshType::active_cell_iterator>
895 Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) ==
false,
896 ExcMessage(
"Halo layer contains locally owned cells"));
897 Assert(contains_artificial_cells<MeshType>(active_halo_layer) ==
false,
898 ExcMessage(
"Halo layer contains artificial cells"));
900 return active_halo_layer;
905 template <
class MeshType>
909 const MeshType &mesh,
910 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
912 const double layer_thickness)
914 std::vector<typename MeshType::active_cell_iterator>
915 subdomain_boundary_cells, active_cell_layer_within_distance;
916 std::vector<bool> vertices_outside_subdomain(
917 mesh.get_triangulation().n_vertices(),
false);
919 const unsigned int spacedim = MeshType::space_dimension;
921 unsigned int n_non_predicate_cells = 0;
929 for (
const auto &cell : mesh.active_cell_iterators())
930 if (!predicate(cell))
932 for (
const unsigned int v : cell->vertex_indices())
933 vertices_outside_subdomain[cell->vertex_index(v)] =
true;
934 n_non_predicate_cells++;
941 if (n_non_predicate_cells == 0 ||
942 n_non_predicate_cells == mesh.get_triangulation().n_active_cells())
943 return std::vector<typename MeshType::active_cell_iterator>();
947 for (
const auto &cell : mesh.active_cell_iterators())
950 for (
const unsigned int v : cell->vertex_indices())
951 if (vertices_outside_subdomain[cell->vertex_index(v)] ==
true)
953 subdomain_boundary_cells.push_back(cell);
966 for (
unsigned int d = 0;
d < spacedim; ++
d)
968 bounding_box.first[
d] -= (layer_thickness + DOUBLE_EPSILON);
969 bounding_box.second[
d] += (layer_thickness + DOUBLE_EPSILON);
972 std::vector<Point<spacedim>>
973 subdomain_boundary_cells_centers;
976 subdomain_boundary_cells_radii;
978 subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size());
979 subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size());
981 for (
typename std::vector<typename MeshType::active_cell_iterator>::
982 const_iterator subdomain_boundary_cell_iterator =
983 subdomain_boundary_cells.begin();
984 subdomain_boundary_cell_iterator != subdomain_boundary_cells.end();
985 ++subdomain_boundary_cell_iterator)
987 const std::pair<Point<spacedim>,
double>
988 &subdomain_boundary_cell_enclosing_ball =
989 (*subdomain_boundary_cell_iterator)->enclosing_ball();
991 subdomain_boundary_cells_centers.push_back(
992 subdomain_boundary_cell_enclosing_ball.first);
993 subdomain_boundary_cells_radii.push_back(
994 subdomain_boundary_cell_enclosing_ball.second);
996 AssertThrow(subdomain_boundary_cells_radii.size() ==
997 subdomain_boundary_cells_centers.size(),
1006 for (
const auto &cell : mesh.active_cell_iterators())
1009 if (predicate(cell))
1012 const std::pair<Point<spacedim>,
double> &cell_enclosing_ball =
1013 cell->enclosing_ball();
1016 cell_enclosing_ball.first;
1017 const double cell_enclosing_ball_radius = cell_enclosing_ball.second;
1019 bool cell_inside =
true;
1021 for (
unsigned int d = 0;
d < spacedim; ++
d)
1023 (cell_enclosing_ball_center[
d] + cell_enclosing_ball_radius >
1024 bounding_box.first[
d]) &&
1025 (cell_enclosing_ball_center[
d] - cell_enclosing_ball_radius <
1026 bounding_box.second[
d]);
1032 for (
unsigned int i = 0; i < subdomain_boundary_cells_radii.size();
1035 subdomain_boundary_cells_centers[i]) <
1036 Utilities::fixed_power<2>(cell_enclosing_ball_radius +
1037 subdomain_boundary_cells_radii[i] +
1038 layer_thickness + DOUBLE_EPSILON))
1040 active_cell_layer_within_distance.push_back(cell);
1045 return active_cell_layer_within_distance;
1050 template <
class MeshType>
1060 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1061 predicate(locally_owned_cell_predicate);
1063 const std::vector<typename MeshType::active_cell_iterator>
1064 ghost_cell_layer_within_distance =
1072 contains_locally_owned_cells<MeshType>(
1073 ghost_cell_layer_within_distance) ==
false,
1075 "Ghost cells within layer_thickness contains locally owned cells."));
1077 contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) ==
1080 "Ghost cells within layer_thickness contains artificial cells. "
1081 "The function compute_ghost_cell_layer_within_distance "
1082 "is probably called while using parallel::distributed::Triangulation. "
1083 "In such case please refer to the description of this function."));
1085 return ghost_cell_layer_within_distance;
1090 template <
class MeshType>
1096 const std::function<
bool(
1097 const typename MeshType::
1098 active_cell_iterator &)>
1101 std::vector<bool> locally_active_vertices_on_subdomain(
1102 mesh.get_triangulation().n_vertices(),
false);
1104 const unsigned int spacedim = MeshType::space_dimension;
1111 for (
const auto &cell : mesh.active_cell_iterators())
1112 if (predicate(cell))
1114 minp = cell->center();
1115 maxp = cell->center();
1121 for (
const auto &cell : mesh.active_cell_iterators())
1122 if (predicate(cell))
1123 for (
const unsigned int v : cell->vertex_indices())
1124 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
1127 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
1129 for (
unsigned int d = 0;
d < spacedim; ++
d)
1136 return std::make_pair(minp, maxp);
1141 template <
typename MeshType>
1143 std::list<std::pair<
1144 typename MeshType::cell_iterator,
1151 ExcMessage(
"The two meshes must be represent triangulations that "
1152 "have the same coarse meshes"));
1159 bool remove_ghost_cells =
false;
1160 #ifdef DEAL_II_WITH_MPI
1162 constexpr
int dim = MeshType::dimension;
1163 constexpr
int spacedim = MeshType::space_dimension;
1165 *
>(&mesh_1.get_triangulation()) !=
nullptr ||
1167 *
>(&mesh_2.get_triangulation()) !=
nullptr)
1169 Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(),
1170 ExcMessage(
"This function can only be used with meshes "
1171 "corresponding to distributed Triangulations when "
1172 "both Triangulations are equal."));
1173 remove_ghost_cells =
true;
1185 using CellList = std::list<std::pair<
typename MeshType::cell_iterator,
1186 typename MeshType::cell_iterator>>;
1190 typename MeshType::cell_iterator cell_1 = mesh_1.begin(0),
1191 cell_2 = mesh_2.begin(0);
1192 for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
1193 cell_list.emplace_back(cell_1, cell_2);
1196 typename CellList::iterator cell_pair = cell_list.begin();
1197 while (cell_pair != cell_list.end())
1201 if (cell_pair->first->has_children() &&
1202 cell_pair->second->has_children())
1204 Assert(cell_pair->first->refinement_case() ==
1205 cell_pair->second->refinement_case(),
1207 for (
unsigned int c = 0; c < cell_pair->first->n_children(); ++c)
1208 cell_list.emplace_back(cell_pair->first->child(c),
1209 cell_pair->second->child(c));
1214 const auto previous_cell_pair = cell_pair;
1216 cell_list.erase(previous_cell_pair);
1221 if (remove_ghost_cells &&
1222 ((cell_pair->first->is_active() &&
1223 !cell_pair->first->is_locally_owned()) ||
1224 (cell_pair->second->is_active() &&
1225 !cell_pair->second->is_locally_owned())))
1228 const auto previous_cell_pair = cell_pair;
1230 cell_list.erase(previous_cell_pair);
1239 for (cell_pair = cell_list.begin(); cell_pair != cell_list.end();
1241 Assert(cell_pair->first->is_active() || cell_pair->second->is_active() ||
1242 (cell_pair->first->refinement_case() !=
1243 cell_pair->second->refinement_case()),
1251 template <
int dim,
int spacedim>
1268 endc = mesh_1.
end(0);
1269 for (; cell_1 != endc; ++cell_1, ++cell_2)
1271 if (cell_1->
n_vertices() != cell_2->n_vertices())
1273 for (
const unsigned int v : cell_1->vertex_indices())
1274 if (cell_1->vertex(v) != cell_2->vertex(v))
1287 template <
typename MeshType>
1292 mesh_2.get_triangulation());
1297 template <
int dim,
int spacedim>
1298 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1304 const double tolerance)
1308 ExcMessage(
"Mapping collection needs to have either size 1 "
1309 "or size equal to the number of elements in "
1310 "the FECollection."));
1312 using cell_iterator =
1315 std::pair<cell_iterator, Point<dim>> best_cell;
1319 if (mapping.
size() == 1)
1321 const std::vector<bool> marked_vertices = {};
1323 mapping[0], mesh, p, marked_vertices, tolerance);
1330 double best_distance = tolerance;
1331 int best_level = -1;
1338 std::vector<cell_iterator> adjacent_cells_tmp =
1346 std::set<cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
1347 adjacent_cells_tmp.end());
1348 std::set<cell_iterator> searched_cells;
1358 unsigned int cells_searched = 0;
1359 while (!found && cells_searched <
n_cells)
1361 typename std::set<cell_iterator>::const_iterator
1364 for (; cell != endc; ++cell)
1369 mapping[(*cell)->active_fe_index()]
1370 .transform_real_to_unit_cell(*cell, p);
1382 if (dist < best_distance || (dist == best_distance &&
1383 (*cell)->level() > best_level))
1386 best_distance = dist;
1387 best_level = (*cell)->level();
1388 best_cell = std::make_pair(*cell, p_cell);
1414 if (!found && cells_searched <
n_cells)
1416 find_active_cell_around_point_internal<dim,
1428 template <
class MeshType>
1431 const typename MeshType::active_cell_iterator &cell)
1433 Assert(cell->is_locally_owned(),
1434 ExcMessage(
"This function only makes sense if the cell for "
1435 "which you are asking for a patch, is locally "
1438 std::vector<typename MeshType::active_cell_iterator> patch;
1439 patch.push_back(cell);
1440 for (
const unsigned int face_number : cell->face_indices())
1441 if (cell->face(face_number)->at_boundary() ==
false)
1443 if (cell->neighbor(face_number)->has_children() ==
false)
1444 patch.push_back(cell->neighbor(face_number));
1449 if (MeshType::dimension > 1)
1451 for (
unsigned int subface = 0;
1452 subface < cell->face(face_number)->n_children();
1455 cell->neighbor_child_on_subface(face_number, subface));
1461 typename MeshType::cell_iterator neighbor =
1462 cell->neighbor(face_number);
1463 while (neighbor->has_children())
1464 neighbor = neighbor->child(1 - face_number);
1466 Assert(neighbor->neighbor(1 - face_number) == cell,
1468 patch.push_back(neighbor);
1476 template <
class Container>
1477 std::vector<typename Container::cell_iterator>
1479 const std::vector<typename Container::active_cell_iterator> &patch)
1483 "Vector containing patch cells should not be an empty vector!"));
1487 int min_level = patch[0]->level();
1489 for (
unsigned int i = 0; i < patch.size(); ++i)
1491 std::set<typename Container::cell_iterator> uniform_cells;
1492 typename std::vector<
1493 typename Container::active_cell_iterator>::const_iterator patch_cell;
1495 for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell)
1500 if ((*patch_cell)->level() == min_level)
1501 uniform_cells.insert(*patch_cell);
1508 typename Container::cell_iterator parent = *patch_cell;
1510 while (parent->level() > min_level)
1511 parent = parent->parent();
1512 uniform_cells.insert(parent);
1516 return std::vector<typename Container::cell_iterator>(uniform_cells.begin(),
1517 uniform_cells.end());
1522 template <
class Container>
1525 const std::vector<typename Container::active_cell_iterator> &patch,
1527 &local_triangulation,
1530 Container::space_dimension>::active_cell_iterator,
1531 typename Container::active_cell_iterator> &patch_to_global_tria_map)
1534 const std::vector<typename Container::cell_iterator> uniform_cells =
1535 get_cells_at_coarsest_common_level<Container>(patch);
1537 local_triangulation.
clear();
1538 std::vector<Point<Container::space_dimension>>
vertices;
1539 const unsigned int n_uniform_cells = uniform_cells.size();
1540 std::vector<CellData<Container::dimension>> cells(n_uniform_cells);
1543 typename std::vector<typename Container::cell_iterator>::const_iterator
1545 for (uniform_cell = uniform_cells.begin();
1546 uniform_cell != uniform_cells.end();
1549 for (
const unsigned int v : (*uniform_cell)->vertex_indices())
1552 (*uniform_cell)->vertex(v);
1553 bool repeat_vertex =
false;
1555 for (
unsigned int m = 0; m < i; ++m)
1559 repeat_vertex =
true;
1560 cells[k].vertices[v] = m;
1564 if (repeat_vertex ==
false)
1567 cells[k].vertices[v] = i;
1578 unsigned int index = 0;
1581 Container::space_dimension>::cell_iterator,
1582 typename Container::cell_iterator>
1583 patch_to_global_tria_map_tmp;
1585 Container::space_dimension>::cell_iterator
1586 coarse_cell = local_triangulation.
begin();
1587 coarse_cell != local_triangulation.
end();
1588 ++coarse_cell, ++index)
1590 patch_to_global_tria_map_tmp.insert(
1591 std::make_pair(coarse_cell, uniform_cells[index]));
1595 Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <=
1596 1
e-15 * coarse_cell->diameter(),
1599 bool refinement_necessary;
1604 refinement_necessary =
false;
1605 for (
const auto &active_tria_cell :
1608 if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
1610 active_tria_cell->set_refine_flag();
1611 refinement_necessary =
true;
1614 for (
unsigned int i = 0; i < patch.size(); ++i)
1619 if (patch_to_global_tria_map_tmp[active_tria_cell] ==
1624 for (
const unsigned int v :
1625 active_tria_cell->vertex_indices())
1626 active_tria_cell->vertex(v) = patch[i]->vertex(v);
1628 Assert(active_tria_cell->center().distance(
1629 patch_to_global_tria_map_tmp[active_tria_cell]
1631 1
e-15 * active_tria_cell->diameter(),
1634 active_tria_cell->set_user_flag();
1640 if (refinement_necessary)
1645 Container::dimension,
1646 Container::space_dimension>::cell_iterator cell =
1647 local_triangulation.
begin();
1648 cell != local_triangulation.
end();
1651 if (patch_to_global_tria_map_tmp.find(cell) !=
1652 patch_to_global_tria_map_tmp.end())
1654 if (cell->has_children())
1661 for (
unsigned int c = 0; c < cell->n_children(); ++c)
1663 if (patch_to_global_tria_map_tmp.find(cell->child(
1664 c)) == patch_to_global_tria_map_tmp.end())
1666 patch_to_global_tria_map_tmp.insert(
1669 patch_to_global_tria_map_tmp[cell]->child(
1691 patch_to_global_tria_map_tmp.erase(cell);
1697 while (refinement_necessary);
1703 Container::space_dimension>::cell_iterator
1704 cell = local_triangulation.
begin();
1705 cell != local_triangulation.
end();
1708 if (cell->user_flag_set())
1710 Assert(patch_to_global_tria_map_tmp.find(cell) !=
1711 patch_to_global_tria_map_tmp.end(),
1714 Assert(cell->center().distance(
1715 patch_to_global_tria_map_tmp[cell]->center()) <=
1716 1
e-15 * cell->diameter(),
1724 Container::space_dimension>::cell_iterator,
1725 typename Container::cell_iterator>::iterator
1726 map_tmp_it = patch_to_global_tria_map_tmp.
begin(),
1727 map_tmp_end = patch_to_global_tria_map_tmp.end();
1731 for (; map_tmp_it != map_tmp_end; ++map_tmp_it)
1732 patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
1737 template <
int dim,
int spacedim>
1740 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1754 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1755 dof_to_set_of_cells_map;
1757 std::vector<types::global_dof_index> local_dof_indices;
1758 std::vector<types::global_dof_index> local_face_dof_indices;
1759 std::vector<types::global_dof_index> local_line_dof_indices;
1763 std::vector<bool> user_flags;
1768 std::map<typename DoFHandler<dim, spacedim>::active_line_iterator,
1770 lines_to_parent_lines_map;
1777 .clear_user_flags();
1782 endc = dof_handler.
end();
1783 for (; cell != endc; ++cell)
1789 if (cell->is_artificial() ==
false)
1791 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
1792 if (cell->line(
l)->has_children())
1793 for (
unsigned int c = 0; c < cell->line(
l)->n_children();
1796 lines_to_parent_lines_map[cell->line(
l)->child(c)] =
1800 cell->line(
l)->child(c)->set_user_flag();
1815 endc = dof_handler.
end();
1816 for (; cell != endc; ++cell)
1821 if (cell->is_artificial() ==
false)
1823 const unsigned int n_dofs_per_cell =
1825 local_dof_indices.resize(n_dofs_per_cell);
1829 cell->get_dof_indices(local_dof_indices);
1830 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i)
1831 dof_to_set_of_cells_map[local_dof_indices[i]].
insert(cell);
1841 for (
const unsigned int f : cell->face_indices())
1843 if (cell->face(f)->has_children())
1845 for (
unsigned int c = 0; c < cell->face(f)->n_children();
1860 Assert(cell->face(f)->child(c)->has_children() ==
false,
1863 const unsigned int n_dofs_per_face =
1865 local_face_dof_indices.resize(n_dofs_per_face);
1867 cell->face(f)->child(c)->get_dof_indices(
1868 local_face_dof_indices);
1869 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1870 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1874 else if ((cell->face(f)->at_boundary() ==
false) &&
1875 (cell->neighbor_is_coarser(f)))
1892 std::pair<unsigned int, unsigned int>
1893 neighbor_face_no_subface_no =
1894 cell->neighbor_of_coarser_neighbor(f);
1895 unsigned int face_no = neighbor_face_no_subface_no.first;
1896 unsigned int subface = neighbor_face_no_subface_no.second;
1898 const unsigned int n_dofs_per_face =
1900 local_face_dof_indices.resize(n_dofs_per_face);
1902 cell->neighbor(f)->face(face_no)->get_dof_indices(
1903 local_face_dof_indices);
1904 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1905 dof_to_set_of_cells_map[local_face_dof_indices[i]].
insert(
1910 for (
unsigned int c = 0;
1911 c < cell->neighbor(f)->face(face_no)->n_children();
1917 const unsigned int n_dofs_per_face =
1919 local_face_dof_indices.resize(n_dofs_per_face);
1924 ->has_children() ==
false,
1929 ->get_dof_indices(local_face_dof_indices);
1930 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1931 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1948 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
1950 if (cell->line(
l)->has_children())
1952 for (
unsigned int c = 0;
1953 c < cell->line(
l)->n_children();
1956 Assert(cell->line(
l)->child(c)->has_children() ==
1962 const unsigned int n_dofs_per_line =
1965 local_line_dof_indices.resize(n_dofs_per_line);
1967 cell->line(
l)->child(c)->get_dof_indices(
1968 local_line_dof_indices);
1969 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1970 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1978 else if (cell->line(
l)->user_flag_set() ==
true)
1982 lines_to_parent_lines_map[cell->line(
l)];
1987 const unsigned int n_dofs_per_line =
1990 local_line_dof_indices.resize(n_dofs_per_line);
1992 parent_line->get_dof_indices(local_line_dof_indices);
1993 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1994 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1997 for (
unsigned int c = 0; c < parent_line->n_children();
2000 Assert(parent_line->child(c)->has_children() ==
2004 const unsigned int n_dofs_per_line =
2007 local_line_dof_indices.resize(n_dofs_per_line);
2009 parent_line->child(c)->get_dof_indices(
2010 local_line_dof_indices);
2011 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
2012 dof_to_set_of_cells_map[local_line_dof_indices[i]]
2029 .load_user_flags(user_flags);
2036 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2037 dof_to_cell_patches;
2041 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>::
2042 iterator it = dof_to_set_of_cells_map.begin(),
2043 it_end = dof_to_set_of_cells_map.end();
2044 for (; it != it_end; ++it)
2045 dof_to_cell_patches[it->first].assign(it->second.begin(),
2048 return dof_to_cell_patches;
2054 template <
typename CellIterator>
2057 std::set<std::pair<CellIterator, unsigned int>> &pairs1,
2060 const unsigned int direction,
2062 const ::Tensor<1, CellIterator::AccessorType::space_dimension>
2066 static const int space_dim = CellIterator::AccessorType::space_dimension;
2072 constexpr
int dim = CellIterator::AccessorType::dimension;
2073 constexpr
int spacedim = CellIterator::AccessorType::space_dimension;
2078 if (!(((pairs1.size() > 0) &&
2079 (
dynamic_cast<const parallel::fullydistributed::
2080 Triangulation<dim, spacedim> *
>(
2081 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2082 ((pairs2.size() > 0) &&
2085 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2086 Assert(pairs1.size() == pairs2.size(),
2087 ExcMessage(
"Unmatched faces on periodic boundaries"));
2091 unsigned int n_matches = 0;
2094 std::bitset<3> orientation;
2095 using PairIterator =
2096 typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator;
2097 for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
2099 for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
2101 const CellIterator cell1 = it1->first;
2102 const CellIterator cell2 = it2->first;
2103 const unsigned int face_idx1 = it1->second;
2104 const unsigned int face_idx2 = it2->second;
2106 cell1->face(face_idx1),
2107 cell2->face(face_idx2),
2116 {cell1, cell2}, {face_idx1, face_idx2}, orientation,
matrix};
2117 matched_pairs.push_back(matched_face);
2131 constexpr
int dim = CellIterator::AccessorType::dimension;
2132 constexpr
int spacedim = CellIterator::AccessorType::space_dimension;
2133 if (!(((pairs1.size() > 0) &&
2134 (
dynamic_cast<const parallel::fullydistributed::
2135 Triangulation<dim, spacedim> *
>(
2136 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2137 ((pairs2.size() > 0) &&
2140 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2141 AssertThrow(n_matches == pairs1.size() && pairs2.size() == 0,
2142 ExcMessage(
"Unmatched faces on periodic boundaries"));
2148 template <
typename MeshType>
2151 const MeshType & mesh,
2153 const
unsigned int direction,
2156 const
Tensor<1, MeshType::space_dimension> &offset,
2159 static const int dim = MeshType::dimension;
2160 static const int space_dim = MeshType::space_dimension;
2170 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2171 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2173 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2174 cell != mesh.end(0);
2177 const typename MeshType::face_iterator face_1 =
2178 cell->face(2 * direction);
2179 const typename MeshType::face_iterator face_2 =
2180 cell->face(2 * direction + 1);
2182 if (face_1->at_boundary() && face_1->boundary_id() == b_id)
2184 const std::pair<typename MeshType::cell_iterator, unsigned int>
2185 pair1 = std::make_pair(cell, 2 * direction);
2186 pairs1.insert(pair1);
2189 if (face_2->at_boundary() && face_2->boundary_id() == b_id)
2191 const std::pair<typename MeshType::cell_iterator, unsigned int>
2192 pair2 = std::make_pair(cell, 2 * direction + 1);
2193 pairs2.insert(pair2);
2197 Assert(pairs1.size() == pairs2.size(),
2198 ExcMessage(
"Unmatched faces on periodic boundaries"));
2200 Assert(pairs1.size() > 0,
2201 ExcMessage(
"No new periodic face pairs have been found. "
2202 "Are you sure that you've selected the correct boundary "
2203 "id's and that the coarsest level mesh is colorized?"));
2206 const unsigned int size_old = matched_pairs.size();
2211 pairs1, pairs2, direction, matched_pairs, offset,
matrix);
2215 const unsigned int size_new = matched_pairs.size();
2216 for (
unsigned int i = size_old; i < size_new; ++i)
2218 Assert(matched_pairs[i].orientation == 1,
2220 "Found a face match with non standard orientation. "
2221 "This function is only suitable for meshes with cells "
2222 "in default orientation"));
2229 template <
typename MeshType>
2232 const MeshType & mesh,
2235 const
unsigned int direction,
2238 const
Tensor<1, MeshType::space_dimension> &offset,
2241 static const int dim = MeshType::dimension;
2242 static const int space_dim = MeshType::space_dimension;
2250 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2251 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2253 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2254 cell != mesh.end(0);
2257 for (
unsigned int i : cell->face_indices())
2259 const typename MeshType::face_iterator face = cell->face(i);
2260 if (face->at_boundary() && face->boundary_id() == b_id1)
2262 const std::pair<typename MeshType::cell_iterator, unsigned int>
2263 pair1 = std::make_pair(cell, i);
2264 pairs1.insert(pair1);
2267 if (face->at_boundary() && face->boundary_id() == b_id2)
2269 const std::pair<typename MeshType::cell_iterator, unsigned int>
2270 pair2 = std::make_pair(cell, i);
2271 pairs2.insert(pair2);
2282 if (!(((pairs1.size() > 0) &&
2285 *
>(&pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2286 ((pairs2.size() > 0) &&
2289 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2290 Assert(pairs1.size() == pairs2.size(),
2291 ExcMessage(
"Unmatched faces on periodic boundaries"));
2294 (pairs1.size() > 0 ||
2297 &mesh.begin()->get_triangulation()) !=
nullptr)),
2298 ExcMessage(
"No new periodic face pairs have been found. "
2299 "Are you sure that you've selected the correct boundary "
2300 "id's and that the coarsest level mesh is colorized?"));
2304 pairs1, pairs2, direction, matched_pairs, offset,
matrix);
2318 template <
int spacedim>
2322 const unsigned int direction,
2332 if (
matrix.m() == spacedim)
2333 for (
unsigned int i = 0; i < spacedim; ++i)
2334 for (
unsigned int j = 0; j < spacedim; ++j)
2335 distance(i) +=
matrix(i, j) * point1(j);
2339 distance += offset - point2;
2341 for (
unsigned int i = 0; i < spacedim; ++i)
2347 if (std::abs(distance(i)) > 1.e-10)
2373 std::array<unsigned int, GeometryInfo<1>::vertices_per_face>;
2374 static inline std::bitset<3>
2386 std::array<unsigned int, GeometryInfo<2>::vertices_per_face>;
2387 static inline std::bitset<3>
2395 static const MATCH_T m_tff = {{0, 1}};
2396 if (matching == m_tff)
2398 static const MATCH_T m_ttf = {{1, 0}};
2399 if (matching == m_ttf)
2412 std::array<unsigned int, GeometryInfo<3>::vertices_per_face>;
2414 static inline std::bitset<3>
2422 static const MATCH_T m_tff = {{0, 1, 2, 3}};
2423 if (matching == m_tff)
2425 static const MATCH_T m_tft = {{1, 3, 0, 2}};
2426 if (matching == m_tft)
2428 static const MATCH_T m_ttf = {{3, 2, 1, 0}};
2429 if (matching == m_ttf)
2431 static const MATCH_T m_ttt = {{2, 0, 3, 1}};
2432 if (matching == m_ttt)
2434 static const MATCH_T m_fff = {{0, 2, 1, 3}};
2435 if (matching == m_fff)
2437 static const MATCH_T m_fft = {{2, 3, 0, 1}};
2438 if (matching == m_fft)
2440 static const MATCH_T m_ftf = {{3, 1, 2, 0}};
2441 if (matching == m_ftf)
2443 static const MATCH_T m_ftt = {{1, 0, 3, 2}};
2444 if (matching == m_ftt)
2455 template <
typename FaceIterator>
2458 std::bitset<3> & orientation,
2459 const FaceIterator & face1,
2460 const FaceIterator & face2,
2461 const unsigned int direction,
2466 ExcMessage(
"The supplied matrix must be a square matrix"));
2468 static const int dim = FaceIterator::AccessorType::dimension;
2472 std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
2476 std::set<unsigned int> face2_vertices;
2477 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2478 face2_vertices.insert(i);
2480 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2482 for (std::set<unsigned int>::iterator it = face2_vertices.begin();
2483 it != face2_vertices.end();
2493 face2_vertices.erase(it);
2500 if (face2_vertices.empty())
2503 return face2_vertices.empty();
2508 template <
typename FaceIterator>
2511 const FaceIterator & face1,
2512 const FaceIterator & face2,
2513 const unsigned int direction,
2518 std::bitset<3> dummy;
2527 #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