54 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
58 const std::vector<bool> & marked_vertices)
69 marked_vertices.size() == 0,
71 marked_vertices.size()));
81 marked_vertices.size() == 0 ||
82 std::equal(marked_vertices.begin(),
83 marked_vertices.end(),
85 [](
bool p,
bool q) { return !p || q; }),
87 "marked_vertices should be a subset of used vertices in the triangulation "
88 "but marked_vertices contains one or more vertices that are not used vertices!"));
96 const std::vector<bool> &used = (marked_vertices.size() == 0) ?
102 std::vector<bool>::const_iterator
first =
103 std::find(used.begin(), used.end(),
true);
109 unsigned int best_vertex = std::distance(used.begin(),
first);
110 double best_dist = (p -
vertices[best_vertex]).norm_square();
114 for (
unsigned int j = best_vertex + 1; j <
vertices.size(); ++j)
117 double dist = (p -
vertices[j]).norm_square();
118 if (dist < best_dist)
130 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
133 const MeshType<dim, spacedim> &mesh,
135 const std::vector<bool> & marked_vertices)
150 marked_vertices.size() == 0,
152 marked_vertices.size()));
162 marked_vertices.size() == 0 ||
163 std::equal(marked_vertices.begin(),
164 marked_vertices.end(),
166 [](
bool p,
bool q) { return !p || q; }),
168 "marked_vertices should be a subset of used vertices in the triangulation "
169 "but marked_vertices contains one or more vertices that are not used vertices!"));
172 if (marked_vertices.size())
175 if (marked_vertices[it->first] ==
false)
190 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
192 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
195 typename ::internal::
196 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
199 const unsigned int vertex)
205 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
212 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
215 typename ::internal::
216 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
217 cell = mesh.begin_active(),
253 for (; cell != endc; ++cell)
255 for (
const unsigned int v : cell->vertex_indices())
256 if (cell->vertex_index(v) == vertex)
268 for (
const auto face :
269 cell->reference_cell().faces_for_given_vertex(v))
270 if (!cell->at_boundary(face) &&
271 cell->neighbor(face)->is_active())
293 for (
unsigned int e = 0;
e < cell->n_lines(); ++
e)
294 if (cell->line(
e)->has_children())
298 if (cell->line(
e)->child(0)->vertex_index(1) == vertex)
322 typename ::internal::
323 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
331 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
333 find_active_cell_around_point_internal(
334 const MeshType<dim, spacedim> &mesh,
336 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>
338 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>
342 typename ::internal::
343 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
346 typename ::internal::
347 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
352 using cell_iterator =
353 typename MeshType<dim, spacedim>::active_cell_iterator;
355 using cell_iterator = typename ::internal::
356 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
364 std::set<cell_iterator> adjacent_cells_new;
366 typename std::set<cell_iterator>::const_iterator cell =
370 for (; cell != endc; ++cell)
372 std::vector<cell_iterator> active_neighbors;
373 get_active_neighbors<MeshType<dim, spacedim>>(*cell,
375 for (
unsigned int i = 0; i < active_neighbors.size(); ++i)
376 if (searched_cells.find(active_neighbors[i]) ==
377 searched_cells.end())
378 adjacent_cells_new.insert(active_neighbors[i]);
382 adjacent_cells_new.end());
391 cell_iterator it = mesh.begin_active();
392 for (; it != mesh.end(); ++it)
393 if (searched_cells.find(it) == searched_cells.end())
404 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
406 typename MeshType<dim, spacedim>::active_cell_iterator
408 typename ::internal::
409 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
413 const std::vector<bool> & marked_vertices,
414 const double tolerance)
416 return find_active_cell_around_point<dim, MeshType, spacedim>(
427 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
429 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
431 std::pair<typename ::internal::
432 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
436 const MeshType<dim, spacedim> &mesh,
438 const std::vector<bool> & marked_vertices,
439 const double tolerance)
441 using active_cell_iterator = typename ::internal::
442 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type;
448 double best_distance = tolerance;
450 std::pair<active_cell_iterator, Point<dim>> best_cell;
453 best_cell.first = mesh.end();
457 std::vector<active_cell_iterator> adjacent_cells_tmp =
466 std::set<active_cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
467 adjacent_cells_tmp.end());
468 std::set<active_cell_iterator> searched_cells;
477 mesh.get_triangulation().n_active_cells();
479 unsigned int cells_searched = 0;
482 typename std::set<active_cell_iterator>::const_iterator
485 for (; cell != endc; ++cell)
487 if ((*cell)->is_artificial() ==
false)
503 if ((dist < best_distance) ||
504 ((dist == best_distance) &&
505 ((*cell)->level() > best_level)))
508 best_distance = dist;
509 best_level = (*cell)->level();
510 best_cell = std::make_pair(*cell, p_cell);
537 if (marked_vertices.size() > 0)
549 find_active_cell_around_point_internal<dim, MeshType, spacedim>(
559 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
561 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
564 std::vector<std::pair<
565 typename ::internal::
566 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
570 const MeshType<dim, spacedim> &mesh,
572 const double tolerance,
573 const std::vector<bool> &marked_vertices)
576 mapping, mesh, p, marked_vertices, tolerance);
578 if (cell_and_point.first == mesh.end())
582 mapping, mesh, p, tolerance, cell_and_point);
587 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
589 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
592 std::vector<std::pair<
593 typename ::internal::
594 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
599 const MeshType<dim, spacedim> &mesh,
601 const double tolerance,
602 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
606 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
611 cells_and_points.push_back(first_cell);
615 const Point<dim> unit_point = cells_and_points.front().second;
616 const auto my_cell = cells_and_points.front().first;
618 unsigned int n_dirs_at_threshold = 0;
620 for (
unsigned int d = 0;
d < dim; ++
d)
622 distance_to_center[
d] = std::abs(unit_point[
d] - 0.5);
623 if (distance_to_center[
d] > 0.5 - tolerance)
625 ++n_dirs_at_threshold;
626 last_point_at_threshold =
d;
630 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
633 if (n_dirs_at_threshold == 1)
635 unsigned int neighbor_index =
636 2 * last_point_at_threshold +
637 (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
638 if (!my_cell->at_boundary(neighbor_index))
639 cells_to_add.push_back(my_cell->neighbor(neighbor_index));
642 else if (n_dirs_at_threshold == dim)
644 unsigned int local_vertex_index = 0;
645 for (
unsigned int d = 0;
d < dim; ++
d)
646 local_vertex_index += (unit_point[
d] > 0.5 ? 1 : 0) <<
d;
647 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
649 mesh, my_cell->vertex_index(local_vertex_index));
650 for (
const auto &cell : cells)
652 cells_to_add.push_back(cell);
659 else if (n_dirs_at_threshold == 2)
662 unsigned int count_vertex_indices = 0;
664 for (
unsigned int d = 0;
d < dim; ++
d)
666 if (distance_to_center[
d] > 0.5 - tolerance)
670 unit_point[
d] > 0.5 ? 1 : 0;
671 count_vertex_indices++;
681 const unsigned int first_vertex =
684 for (
unsigned int d = 0;
d < 2; ++
d)
688 my_cell->vertex_index(first_vertex + (
d << free_direction)));
689 for (
const auto &cell : tentative_cells)
691 bool cell_not_yet_present =
true;
692 for (
const auto &other_cell : cells_to_add)
693 if (cell == other_cell)
695 cell_not_yet_present =
false;
698 if (cell_not_yet_present)
699 cells_to_add.push_back(cell);
704 const double original_distance_to_unit_cell =
706 for (
const auto &cell : cells_to_add)
714 original_distance_to_unit_cell + tolerance)
715 cells_and_points.emplace_back(cell, p_unit);
722 cells_and_points.begin(),
723 cells_and_points.end(),
724 [](
const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
726 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
729 return cells_and_points;
734 template <
class MeshType>
735 std::vector<typename MeshType::active_cell_iterator>
737 const MeshType &mesh,
738 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
741 std::vector<typename MeshType::active_cell_iterator> active_halo_layer;
742 std::vector<bool> locally_active_vertices_on_subdomain(
743 mesh.get_triangulation().n_vertices(),
false);
748 for (
const auto &cell : mesh.active_cell_iterators())
750 for (
const auto v : cell->vertex_indices())
751 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
true;
756 for (
const auto &cell : mesh.active_cell_iterators())
757 if (!predicate(cell))
758 for (
const auto v : cell->vertex_indices())
759 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
762 active_halo_layer.push_back(cell);
766 return active_halo_layer;
771 template <
class MeshType>
772 std::vector<typename MeshType::cell_iterator>
774 const MeshType &mesh,
775 const std::function<
bool(
const typename MeshType::cell_iterator &)>
777 const unsigned int level)
779 std::vector<typename MeshType::cell_iterator> level_halo_layer;
780 std::vector<bool> locally_active_vertices_on_level_subdomain(
781 mesh.get_triangulation().n_vertices(),
false);
786 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
787 cell != mesh.end(
level);
790 for (
const unsigned int v : cell->vertex_indices())
791 locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] =
797 for (
typename MeshType::cell_iterator cell = mesh.begin(
level);
798 cell != mesh.end(
level);
800 if (!predicate(cell))
801 for (
const unsigned int v : cell->vertex_indices())
802 if (locally_active_vertices_on_level_subdomain[cell->vertex_index(
805 level_halo_layer.push_back(cell);
809 return level_halo_layer;
815 template <
class MeshType>
817 contains_locally_owned_cells(
818 const std::vector<typename MeshType::active_cell_iterator> &cells)
820 for (
typename std::vector<
821 typename MeshType::active_cell_iterator>::const_iterator it =
826 if ((*it)->is_locally_owned())
832 template <
class MeshType>
834 contains_artificial_cells(
835 const std::vector<typename MeshType::active_cell_iterator> &cells)
837 for (
typename std::vector<
838 typename MeshType::active_cell_iterator>::const_iterator it =
843 if ((*it)->is_artificial())
852 template <
class MeshType>
853 std::vector<typename MeshType::active_cell_iterator>
856 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
859 const std::vector<typename MeshType::active_cell_iterator>
864 Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) ==
false,
865 ExcMessage(
"Halo layer contains locally owned cells"));
866 Assert(contains_artificial_cells<MeshType>(active_halo_layer) ==
false,
867 ExcMessage(
"Halo layer contains artificial cells"));
869 return active_halo_layer;
874 template <
class MeshType>
875 std::vector<typename MeshType::active_cell_iterator>
877 const MeshType &mesh,
878 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
880 const double layer_thickness)
882 std::vector<typename MeshType::active_cell_iterator>
883 subdomain_boundary_cells, active_cell_layer_within_distance;
884 std::vector<bool> vertices_outside_subdomain(
885 mesh.get_triangulation().n_vertices(),
false);
887 const unsigned int spacedim = MeshType::space_dimension;
889 unsigned int n_non_predicate_cells = 0;
897 for (
const auto &cell : mesh.active_cell_iterators())
898 if (!predicate(cell))
900 for (
const unsigned int v : cell->vertex_indices())
901 vertices_outside_subdomain[cell->vertex_index(v)] =
true;
902 n_non_predicate_cells++;
909 if (n_non_predicate_cells == 0 ||
910 n_non_predicate_cells == mesh.get_triangulation().n_active_cells())
911 return std::vector<typename MeshType::active_cell_iterator>();
915 for (
const auto &cell : mesh.active_cell_iterators())
918 for (
const unsigned int v : cell->vertex_indices())
919 if (vertices_outside_subdomain[cell->vertex_index(v)] ==
true)
921 subdomain_boundary_cells.push_back(cell);
934 for (
unsigned int d = 0;
d < spacedim; ++
d)
936 bounding_box.first[
d] -= (layer_thickness + DOUBLE_EPSILON);
937 bounding_box.second[
d] += (layer_thickness + DOUBLE_EPSILON);
940 std::vector<Point<spacedim>>
941 subdomain_boundary_cells_centers;
944 subdomain_boundary_cells_radii;
946 subdomain_boundary_cells_centers.reserve(subdomain_boundary_cells.size());
947 subdomain_boundary_cells_radii.reserve(subdomain_boundary_cells.size());
949 for (
typename std::vector<typename MeshType::active_cell_iterator>::
950 const_iterator subdomain_boundary_cell_iterator =
951 subdomain_boundary_cells.begin();
952 subdomain_boundary_cell_iterator != subdomain_boundary_cells.end();
953 ++subdomain_boundary_cell_iterator)
955 const std::pair<Point<spacedim>,
double>
956 &subdomain_boundary_cell_enclosing_ball =
957 (*subdomain_boundary_cell_iterator)->enclosing_ball();
959 subdomain_boundary_cells_centers.push_back(
960 subdomain_boundary_cell_enclosing_ball.first);
961 subdomain_boundary_cells_radii.push_back(
962 subdomain_boundary_cell_enclosing_ball.second);
964 AssertThrow(subdomain_boundary_cells_radii.size() ==
965 subdomain_boundary_cells_centers.size(),
974 for (
const auto &cell : mesh.active_cell_iterators())
980 const std::pair<Point<spacedim>,
double> &cell_enclosing_ball =
981 cell->enclosing_ball();
984 cell_enclosing_ball.first;
985 const double cell_enclosing_ball_radius = cell_enclosing_ball.second;
987 bool cell_inside =
true;
989 for (
unsigned int d = 0;
d < spacedim; ++
d)
991 (cell_enclosing_ball_center[
d] + cell_enclosing_ball_radius >
992 bounding_box.first[
d]) &&
993 (cell_enclosing_ball_center[
d] - cell_enclosing_ball_radius <
994 bounding_box.second[
d]);
1000 for (
unsigned int i = 0; i < subdomain_boundary_cells_radii.size();
1003 subdomain_boundary_cells_centers[i]) <
1004 Utilities::fixed_power<2>(cell_enclosing_ball_radius +
1005 subdomain_boundary_cells_radii[i] +
1006 layer_thickness + DOUBLE_EPSILON))
1008 active_cell_layer_within_distance.push_back(cell);
1013 return active_cell_layer_within_distance;
1018 template <
class MeshType>
1019 std::vector<typename MeshType::active_cell_iterator>
1021 const double layer_thickness)
1024 std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1025 predicate(locally_owned_cell_predicate);
1027 const std::vector<typename MeshType::active_cell_iterator>
1028 ghost_cell_layer_within_distance =
1036 contains_locally_owned_cells<MeshType>(
1037 ghost_cell_layer_within_distance) ==
false,
1039 "Ghost cells within layer_thickness contains locally owned cells."));
1041 contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) ==
1044 "Ghost cells within layer_thickness contains artificial cells. "
1045 "The function compute_ghost_cell_layer_within_distance "
1046 "is probably called while using parallel::distributed::Triangulation. "
1047 "In such case please refer to the description of this function."));
1049 return ghost_cell_layer_within_distance;
1054 template <
class MeshType>
1057 const MeshType &mesh,
1058 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1061 std::vector<bool> locally_active_vertices_on_subdomain(
1062 mesh.get_triangulation().n_vertices(),
false);
1064 const unsigned int spacedim = MeshType::space_dimension;
1071 for (
const auto &cell : mesh.active_cell_iterators())
1072 if (predicate(cell))
1074 minp = cell->center();
1075 maxp = cell->center();
1081 for (
const auto &cell : mesh.active_cell_iterators())
1082 if (predicate(cell))
1083 for (
const unsigned int v : cell->vertex_indices())
1084 if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
1087 locally_active_vertices_on_subdomain[cell->vertex_index(v)] =
1089 for (
unsigned int d = 0;
d < spacedim; ++
d)
1096 return std::make_pair(minp, maxp);
1101 template <
typename MeshType>
1102 std::list<std::pair<
typename MeshType::cell_iterator,
1103 typename MeshType::cell_iterator>>
1107 ExcMessage(
"The two meshes must be represent triangulations that "
1108 "have the same coarse meshes"));
1115 bool remove_ghost_cells =
false;
1116 #ifdef DEAL_II_WITH_MPI
1118 constexpr
int dim = MeshType::dimension;
1119 constexpr
int spacedim = MeshType::space_dimension;
1121 *
>(&mesh_1.get_triangulation()) !=
nullptr ||
1123 *
>(&mesh_2.get_triangulation()) !=
nullptr)
1125 Assert(&mesh_1.get_triangulation() == &mesh_2.get_triangulation(),
1126 ExcMessage(
"This function can only be used with meshes "
1127 "corresponding to distributed Triangulations when "
1128 "both Triangulations are equal."));
1129 remove_ghost_cells =
true;
1141 using CellList = std::list<std::pair<
typename MeshType::cell_iterator,
1142 typename MeshType::cell_iterator>>;
1146 typename MeshType::cell_iterator cell_1 = mesh_1.begin(0),
1147 cell_2 = mesh_2.begin(0);
1148 for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
1149 cell_list.emplace_back(cell_1, cell_2);
1152 typename CellList::iterator cell_pair = cell_list.begin();
1153 while (cell_pair != cell_list.end())
1157 if (cell_pair->first->has_children() &&
1158 cell_pair->second->has_children())
1160 Assert(cell_pair->first->refinement_case() ==
1161 cell_pair->second->refinement_case(),
1163 for (
unsigned int c = 0; c < cell_pair->first->n_children(); ++c)
1164 cell_list.emplace_back(cell_pair->first->child(c),
1165 cell_pair->second->child(c));
1170 const auto previous_cell_pair = cell_pair;
1172 cell_list.erase(previous_cell_pair);
1177 if (remove_ghost_cells &&
1178 ((cell_pair->first->is_active() &&
1179 !cell_pair->first->is_locally_owned()) ||
1180 (cell_pair->second->is_active() &&
1181 !cell_pair->second->is_locally_owned())))
1184 const auto previous_cell_pair = cell_pair;
1186 cell_list.erase(previous_cell_pair);
1195 for (cell_pair = cell_list.begin(); cell_pair != cell_list.end();
1197 Assert(cell_pair->first->is_active() || cell_pair->second->is_active() ||
1198 (cell_pair->first->refinement_case() !=
1199 cell_pair->second->refinement_case()),
1207 template <
int dim,
int spacedim>
1224 endc = mesh_1.
end(0);
1225 for (; cell_1 != endc; ++cell_1, ++cell_2)
1227 if (cell_1->
n_vertices() != cell_2->n_vertices())
1229 for (
const unsigned int v : cell_1->vertex_indices())
1230 if (cell_1->vertex(v) != cell_2->vertex(v))
1243 template <
typename MeshType>
1248 mesh_2.get_triangulation());
1253 template <
int dim,
int spacedim>
1254 std::pair<typename ::DoFHandler<dim, spacedim>::active_cell_iterator,
1260 const double tolerance)
1264 ExcMessage(
"Mapping collection needs to have either size 1 "
1265 "or size equal to the number of elements in "
1266 "the FECollection."));
1268 using cell_iterator =
1269 typename ::DoFHandler<dim, spacedim>::active_cell_iterator;
1271 std::pair<cell_iterator, Point<dim>> best_cell;
1275 if (mapping.
size() == 1)
1277 const std::vector<bool> marked_vertices = {};
1279 mapping[0], mesh, p, marked_vertices, tolerance);
1286 double best_distance = tolerance;
1287 int best_level = -1;
1294 std::vector<cell_iterator> adjacent_cells_tmp =
1302 std::set<cell_iterator>
adjacent_cells(adjacent_cells_tmp.begin(),
1303 adjacent_cells_tmp.end());
1304 std::set<cell_iterator> searched_cells;
1314 unsigned int cells_searched = 0;
1315 while (!found && cells_searched <
n_cells)
1317 typename std::set<cell_iterator>::const_iterator
1320 for (; cell != endc; ++cell)
1325 mapping[(*cell)->active_fe_index()]
1326 .transform_real_to_unit_cell(*cell, p);
1338 if (dist < best_distance || (dist == best_distance &&
1339 (*cell)->level() > best_level))
1342 best_distance = dist;
1343 best_level = (*cell)->level();
1344 best_cell = std::make_pair(*cell, p_cell);
1370 if (!found && cells_searched <
n_cells)
1372 find_active_cell_around_point_internal<dim,
1384 template <
class MeshType>
1385 std::vector<typename MeshType::active_cell_iterator>
1388 Assert(cell->is_locally_owned(),
1389 ExcMessage(
"This function only makes sense if the cell for "
1390 "which you are asking for a patch, is locally "
1393 std::vector<typename MeshType::active_cell_iterator> patch;
1394 patch.push_back(cell);
1395 for (
const unsigned int face_number : cell->face_indices())
1396 if (cell->face(face_number)->at_boundary() ==
false)
1398 if (cell->neighbor(face_number)->has_children() ==
false)
1399 patch.push_back(cell->neighbor(face_number));
1404 if (MeshType::dimension > 1)
1406 for (
unsigned int subface = 0;
1407 subface < cell->face(face_number)->n_children();
1410 cell->neighbor_child_on_subface(face_number, subface));
1416 typename MeshType::cell_iterator neighbor =
1417 cell->neighbor(face_number);
1418 while (neighbor->has_children())
1419 neighbor = neighbor->child(1 - face_number);
1421 Assert(neighbor->neighbor(1 - face_number) == cell,
1423 patch.push_back(neighbor);
1431 template <
class Container>
1432 std::vector<typename Container::cell_iterator>
1434 const std::vector<typename Container::active_cell_iterator> &patch)
1438 "Vector containing patch cells should not be an empty vector!"));
1442 int min_level = patch[0]->level();
1444 for (
unsigned int i = 0; i < patch.size(); ++i)
1446 std::set<typename Container::cell_iterator> uniform_cells;
1447 typename std::vector<
1448 typename Container::active_cell_iterator>::const_iterator patch_cell;
1450 for (patch_cell = patch.begin(); patch_cell != patch.end(); ++patch_cell)
1455 if ((*patch_cell)->level() == min_level)
1456 uniform_cells.insert(*patch_cell);
1463 typename Container::cell_iterator parent = *patch_cell;
1465 while (parent->level() > min_level)
1466 parent = parent->parent();
1467 uniform_cells.insert(parent);
1471 return std::vector<typename Container::cell_iterator>(uniform_cells.begin(),
1472 uniform_cells.end());
1477 template <
class Container>
1480 const std::vector<typename Container::active_cell_iterator> &patch,
1482 &local_triangulation,
1485 Container::space_dimension>::active_cell_iterator,
1486 typename Container::active_cell_iterator> &patch_to_global_tria_map)
1489 const std::vector<typename Container::cell_iterator> uniform_cells =
1490 get_cells_at_coarsest_common_level<Container>(patch);
1492 local_triangulation.
clear();
1493 std::vector<Point<Container::space_dimension>>
vertices;
1494 const unsigned int n_uniform_cells = uniform_cells.size();
1495 std::vector<CellData<Container::dimension>> cells(n_uniform_cells);
1498 typename std::vector<typename Container::cell_iterator>::const_iterator
1500 for (uniform_cell = uniform_cells.begin();
1501 uniform_cell != uniform_cells.end();
1504 for (
const unsigned int v : (*uniform_cell)->vertex_indices())
1507 (*uniform_cell)->vertex(v);
1508 bool repeat_vertex =
false;
1510 for (
unsigned int m = 0; m < i; ++m)
1514 repeat_vertex =
true;
1515 cells[k].vertices[v] = m;
1519 if (repeat_vertex ==
false)
1522 cells[k].vertices[v] = i;
1533 unsigned int index = 0;
1536 Container::space_dimension>::cell_iterator,
1537 typename Container::cell_iterator>
1538 patch_to_global_tria_map_tmp;
1540 Container::space_dimension>::cell_iterator
1541 coarse_cell = local_triangulation.
begin();
1542 coarse_cell != local_triangulation.
end();
1543 ++coarse_cell, ++index)
1545 patch_to_global_tria_map_tmp.insert(
1546 std::make_pair(coarse_cell, uniform_cells[index]));
1550 Assert(coarse_cell->center().distance(uniform_cells[index]->center()) <=
1551 1
e-15 * coarse_cell->diameter(),
1554 bool refinement_necessary;
1559 refinement_necessary =
false;
1560 for (
const auto &active_tria_cell :
1563 if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
1565 active_tria_cell->set_refine_flag();
1566 refinement_necessary =
true;
1569 for (
unsigned int i = 0; i < patch.size(); ++i)
1574 if (patch_to_global_tria_map_tmp[active_tria_cell] ==
1579 for (
const unsigned int v :
1580 active_tria_cell->vertex_indices())
1581 active_tria_cell->vertex(v) = patch[i]->vertex(v);
1583 Assert(active_tria_cell->center().distance(
1584 patch_to_global_tria_map_tmp[active_tria_cell]
1586 1
e-15 * active_tria_cell->diameter(),
1589 active_tria_cell->set_user_flag();
1595 if (refinement_necessary)
1600 Container::dimension,
1601 Container::space_dimension>::cell_iterator cell =
1602 local_triangulation.
begin();
1603 cell != local_triangulation.
end();
1606 if (patch_to_global_tria_map_tmp.find(cell) !=
1607 patch_to_global_tria_map_tmp.end())
1609 if (cell->has_children())
1616 for (
unsigned int c = 0; c < cell->n_children(); ++c)
1618 if (patch_to_global_tria_map_tmp.find(cell->child(
1619 c)) == patch_to_global_tria_map_tmp.end())
1621 patch_to_global_tria_map_tmp.insert(
1624 patch_to_global_tria_map_tmp[cell]->child(
1646 patch_to_global_tria_map_tmp.erase(cell);
1652 while (refinement_necessary);
1658 Container::space_dimension>::cell_iterator
1659 cell = local_triangulation.
begin();
1660 cell != local_triangulation.
end();
1663 if (cell->user_flag_set())
1665 Assert(patch_to_global_tria_map_tmp.find(cell) !=
1666 patch_to_global_tria_map_tmp.end(),
1669 Assert(cell->center().distance(
1670 patch_to_global_tria_map_tmp[cell]->center()) <=
1671 1
e-15 * cell->diameter(),
1679 Container::space_dimension>::cell_iterator,
1680 typename Container::cell_iterator>::iterator
1681 map_tmp_it = patch_to_global_tria_map_tmp.
begin(),
1682 map_tmp_end = patch_to_global_tria_map_tmp.end();
1686 for (; map_tmp_it != map_tmp_end; ++map_tmp_it)
1687 patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
1692 template <
int dim,
int spacedim>
1695 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1709 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1710 dof_to_set_of_cells_map;
1712 std::vector<types::global_dof_index> local_dof_indices;
1713 std::vector<types::global_dof_index> local_face_dof_indices;
1714 std::vector<types::global_dof_index> local_line_dof_indices;
1718 std::vector<bool> user_flags;
1723 std::map<typename DoFHandler<dim, spacedim>::active_line_iterator,
1725 lines_to_parent_lines_map;
1732 .clear_user_flags();
1737 endc = dof_handler.
end();
1738 for (; cell != endc; ++cell)
1744 if (cell->is_artificial() ==
false)
1746 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
1747 if (cell->line(
l)->has_children())
1748 for (
unsigned int c = 0; c < cell->line(
l)->n_children();
1751 lines_to_parent_lines_map[cell->line(
l)->child(c)] =
1755 cell->line(
l)->child(c)->set_user_flag();
1770 endc = dof_handler.
end();
1771 for (; cell != endc; ++cell)
1776 if (cell->is_artificial() ==
false)
1778 const unsigned int n_dofs_per_cell =
1780 local_dof_indices.resize(n_dofs_per_cell);
1784 cell->get_dof_indices(local_dof_indices);
1785 for (
unsigned int i = 0; i < n_dofs_per_cell; ++i)
1786 dof_to_set_of_cells_map[local_dof_indices[i]].
insert(cell);
1796 for (
const unsigned int f : cell->face_indices())
1798 if (cell->face(f)->has_children())
1800 for (
unsigned int c = 0; c < cell->face(f)->n_children();
1815 Assert(cell->face(f)->child(c)->has_children() ==
false,
1818 const unsigned int n_dofs_per_face =
1820 local_face_dof_indices.resize(n_dofs_per_face);
1822 cell->face(f)->child(c)->get_dof_indices(
1823 local_face_dof_indices);
1824 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1825 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1829 else if ((cell->face(f)->at_boundary() ==
false) &&
1830 (cell->neighbor_is_coarser(f)))
1847 std::pair<unsigned int, unsigned int>
1848 neighbor_face_no_subface_no =
1849 cell->neighbor_of_coarser_neighbor(f);
1850 unsigned int face_no = neighbor_face_no_subface_no.first;
1851 unsigned int subface = neighbor_face_no_subface_no.second;
1853 const unsigned int n_dofs_per_face =
1855 local_face_dof_indices.resize(n_dofs_per_face);
1857 cell->neighbor(f)->face(face_no)->get_dof_indices(
1858 local_face_dof_indices);
1859 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1860 dof_to_set_of_cells_map[local_face_dof_indices[i]].
insert(
1865 for (
unsigned int c = 0;
1866 c < cell->neighbor(f)->face(face_no)->n_children();
1872 const unsigned int n_dofs_per_face =
1874 local_face_dof_indices.resize(n_dofs_per_face);
1879 ->has_children() ==
false,
1884 ->get_dof_indices(local_face_dof_indices);
1885 for (
unsigned int i = 0; i < n_dofs_per_face; ++i)
1886 dof_to_set_of_cells_map[local_face_dof_indices[i]]
1903 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
1905 if (cell->line(
l)->has_children())
1907 for (
unsigned int c = 0;
1908 c < cell->line(
l)->n_children();
1911 Assert(cell->line(
l)->child(c)->has_children() ==
1917 const unsigned int n_dofs_per_line =
1920 local_line_dof_indices.resize(n_dofs_per_line);
1922 cell->line(
l)->child(c)->get_dof_indices(
1923 local_line_dof_indices);
1924 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1925 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1933 else if (cell->line(
l)->user_flag_set() ==
true)
1937 lines_to_parent_lines_map[cell->line(
l)];
1942 const unsigned int n_dofs_per_line =
1945 local_line_dof_indices.resize(n_dofs_per_line);
1947 parent_line->get_dof_indices(local_line_dof_indices);
1948 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1949 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1952 for (
unsigned int c = 0; c < parent_line->n_children();
1955 Assert(parent_line->child(c)->has_children() ==
1959 const unsigned int n_dofs_per_line =
1962 local_line_dof_indices.resize(n_dofs_per_line);
1964 parent_line->child(c)->get_dof_indices(
1965 local_line_dof_indices);
1966 for (
unsigned int i = 0; i < n_dofs_per_line; ++i)
1967 dof_to_set_of_cells_map[local_line_dof_indices[i]]
1984 .load_user_flags(user_flags);
1991 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
1992 dof_to_cell_patches;
1996 std::set<typename DoFHandler<dim, spacedim>::active_cell_iterator>>::
1997 iterator it = dof_to_set_of_cells_map.begin(),
1998 it_end = dof_to_set_of_cells_map.end();
1999 for (; it != it_end; ++it)
2000 dof_to_cell_patches[it->first].assign(it->second.begin(),
2003 return dof_to_cell_patches;
2009 template <
typename CellIterator>
2012 std::set<std::pair<CellIterator, unsigned int>> &pairs1,
2015 const unsigned int direction,
2017 const ::Tensor<1, CellIterator::AccessorType::space_dimension>
2021 static const int space_dim = CellIterator::AccessorType::space_dimension;
2027 constexpr
int dim = CellIterator::AccessorType::dimension;
2028 constexpr
int spacedim = CellIterator::AccessorType::space_dimension;
2033 if (!(((pairs1.size() > 0) &&
2034 (
dynamic_cast<const parallel::fullydistributed::
2035 Triangulation<dim, spacedim> *
>(
2036 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2037 ((pairs2.size() > 0) &&
2040 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2041 Assert(pairs1.size() == pairs2.size(),
2042 ExcMessage(
"Unmatched faces on periodic boundaries"));
2046 unsigned int n_matches = 0;
2049 std::bitset<3> orientation;
2050 using PairIterator =
2051 typename std::set<std::pair<CellIterator, unsigned int>>::const_iterator;
2052 for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
2054 for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
2056 const CellIterator cell1 = it1->first;
2057 const CellIterator cell2 = it2->first;
2058 const unsigned int face_idx1 = it1->second;
2059 const unsigned int face_idx2 = it2->second;
2061 cell1->face(face_idx1),
2062 cell2->face(face_idx2),
2071 {cell1, cell2}, {face_idx1, face_idx2}, orientation,
matrix};
2072 matched_pairs.push_back(matched_face);
2086 constexpr
int dim = CellIterator::AccessorType::dimension;
2087 constexpr
int spacedim = CellIterator::AccessorType::space_dimension;
2088 if (!(((pairs1.size() > 0) &&
2089 (
dynamic_cast<const parallel::fullydistributed::
2090 Triangulation<dim, spacedim> *
>(
2091 &pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2092 ((pairs2.size() > 0) &&
2095 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2096 AssertThrow(n_matches == pairs1.size() && pairs2.size() == 0,
2097 ExcMessage(
"Unmatched faces on periodic boundaries"));
2103 template <
typename MeshType>
2106 const MeshType & mesh,
2108 const unsigned int direction,
2114 static const int dim = MeshType::dimension;
2115 static const int space_dim = MeshType::space_dimension;
2125 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2126 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2128 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2129 cell != mesh.end(0);
2132 const typename MeshType::face_iterator face_1 =
2133 cell->face(2 * direction);
2134 const typename MeshType::face_iterator face_2 =
2135 cell->face(2 * direction + 1);
2137 if (face_1->at_boundary() && face_1->boundary_id() == b_id)
2139 const std::pair<typename MeshType::cell_iterator, unsigned int>
2140 pair1 = std::make_pair(cell, 2 * direction);
2141 pairs1.insert(pair1);
2144 if (face_2->at_boundary() && face_2->boundary_id() == b_id)
2146 const std::pair<typename MeshType::cell_iterator, unsigned int>
2147 pair2 = std::make_pair(cell, 2 * direction + 1);
2148 pairs2.insert(pair2);
2152 Assert(pairs1.size() == pairs2.size(),
2153 ExcMessage(
"Unmatched faces on periodic boundaries"));
2155 Assert(pairs1.size() > 0,
2156 ExcMessage(
"No new periodic face pairs have been found. "
2157 "Are you sure that you've selected the correct boundary "
2158 "id's and that the coarsest level mesh is colorized?"));
2161 const unsigned int size_old = matched_pairs.size();
2166 pairs1, pairs2, direction, matched_pairs, offset,
matrix);
2170 const unsigned int size_new = matched_pairs.size();
2171 for (
unsigned int i = size_old; i < size_new; ++i)
2173 Assert(matched_pairs[i].orientation == 1,
2175 "Found a face match with non standard orientation. "
2176 "This function is only suitable for meshes with cells "
2177 "in default orientation"));
2184 template <
typename MeshType>
2187 const MeshType & mesh,
2190 const unsigned int direction,
2196 static const int dim = MeshType::dimension;
2197 static const int space_dim = MeshType::space_dimension;
2205 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs1;
2206 std::set<std::pair<typename MeshType::cell_iterator, unsigned int>> pairs2;
2208 for (
typename MeshType::cell_iterator cell = mesh.begin(0);
2209 cell != mesh.end(0);
2212 for (
unsigned int i : cell->face_indices())
2214 const typename MeshType::face_iterator face = cell->face(i);
2215 if (face->at_boundary() && face->boundary_id() == b_id1)
2217 const std::pair<typename MeshType::cell_iterator, unsigned int>
2218 pair1 = std::make_pair(cell, i);
2219 pairs1.insert(pair1);
2222 if (face->at_boundary() && face->boundary_id() == b_id2)
2224 const std::pair<typename MeshType::cell_iterator, unsigned int>
2225 pair2 = std::make_pair(cell, i);
2226 pairs2.insert(pair2);
2237 if (!(((pairs1.size() > 0) &&
2240 *
>(&pairs1.begin()->first->get_triangulation()) !=
nullptr)) ||
2241 ((pairs2.size() > 0) &&
2244 *
>(&pairs2.begin()->first->get_triangulation()) !=
nullptr))))
2245 Assert(pairs1.size() == pairs2.size(),
2246 ExcMessage(
"Unmatched faces on periodic boundaries"));
2249 (pairs1.size() > 0 ||
2252 &mesh.begin()->get_triangulation()) !=
nullptr)),
2253 ExcMessage(
"No new periodic face pairs have been found. "
2254 "Are you sure that you've selected the correct boundary "
2255 "id's and that the coarsest level mesh is colorized?"));
2259 pairs1, pairs2, direction, matched_pairs, offset,
matrix);
2273 template <
int spacedim>
2277 const unsigned int direction,
2287 if (
matrix.m() == spacedim)
2288 for (
unsigned int i = 0; i < spacedim; ++i)
2289 for (
unsigned int j = 0; j < spacedim; ++j)
2290 distance(i) +=
matrix(i, j) * point1(j);
2294 distance += offset - point2;
2296 for (
unsigned int i = 0; i < spacedim; ++i)
2302 if (std::abs(distance(i)) > 1.e-10)
2328 std::array<unsigned int, GeometryInfo<1>::vertices_per_face>;
2329 static inline std::bitset<3>
2341 std::array<unsigned int, GeometryInfo<2>::vertices_per_face>;
2342 static inline std::bitset<3>
2350 static const MATCH_T m_tff = {{0, 1}};
2351 if (matching == m_tff)
2353 static const MATCH_T m_ttf = {{1, 0}};
2354 if (matching == m_ttf)
2367 std::array<unsigned int, GeometryInfo<3>::vertices_per_face>;
2369 static inline std::bitset<3>
2377 static const MATCH_T m_tff = {{0, 1, 2, 3}};
2378 if (matching == m_tff)
2380 static const MATCH_T m_tft = {{1, 3, 0, 2}};
2381 if (matching == m_tft)
2383 static const MATCH_T m_ttf = {{3, 2, 1, 0}};
2384 if (matching == m_ttf)
2386 static const MATCH_T m_ttt = {{2, 0, 3, 1}};
2387 if (matching == m_ttt)
2389 static const MATCH_T m_fff = {{0, 2, 1, 3}};
2390 if (matching == m_fff)
2392 static const MATCH_T m_fft = {{2, 3, 0, 1}};
2393 if (matching == m_fft)
2395 static const MATCH_T m_ftf = {{3, 1, 2, 0}};
2396 if (matching == m_ftf)
2398 static const MATCH_T m_ftt = {{1, 0, 3, 2}};
2399 if (matching == m_ftt)
2410 template <
typename FaceIterator>
2413 std::bitset<3> & orientation,
2414 const FaceIterator & face1,
2415 const FaceIterator & face2,
2416 const unsigned int direction,
2421 ExcMessage(
"The supplied matrix must be a square matrix"));
2423 static const int dim = FaceIterator::AccessorType::dimension;
2427 std::array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;
2431 std::set<unsigned int> face2_vertices;
2432 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2433 face2_vertices.insert(i);
2435 for (
unsigned int i = 0; i < face1->n_vertices(); ++i)
2437 for (std::set<unsigned int>::iterator it = face2_vertices.begin();
2438 it != face2_vertices.end();
2448 face2_vertices.erase(it);
2455 if (face2_vertices.empty())
2458 return face2_vertices.empty();
2463 template <
typename FaceIterator>
2466 const FaceIterator & face1,
2467 const FaceIterator & face2,
2468 const unsigned int direction,
2473 std::bitset<3> dummy;
2482 #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 unsigned int 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
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
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_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)
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
unsigned int global_dof_index
static double distance_to_unit_cell(const Point< dim > &p)
const ::Triangulation< dim, spacedim > & tria