18 #include <deal.II/base/mpi.templates.h>
23 #ifdef DEAL_II_WITH_ARBORX
27 template <
int dim,
typename Number>
35 template <
int dim,
typename Number>
38 template <
typename QueryType>
39 std::pair<std::vector<std::pair<int, int>>, std::vector<int>>
40 query(
const QueryType &queries);
42 class BoundingBoxIntersectPredicate
47 #ifdef DEAL_II_WITH_CGAL
91 #include <boost/random/mersenne_twister.hpp>
92 #include <boost/random/uniform_real_distribution.hpp>
102 #include <unordered_map>
112 template <
int spacedim>
170 template <
int spacedim>
189 template <
int dim,
int spacedim>
199 template <
int dim,
int spacedim>
208 "GridTools::rotate() is only available for spacedim = 2."));
243 const unsigned int axis,
255 template <
int dim,
int spacedim>
277 const unsigned int n_dofs = S.
n();
289 const auto constrained_rhs =
291 solver.
solve(SF, u, constrained_rhs, prec);
304 const bool solve_for_absolute_positions)
339 std::array<AffineConstraints<double>, dim> constraints;
340 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
350 for (
const unsigned int vertex_no : cell->vertex_indices())
352 const unsigned int vertex_index = cell->vertex_index(vertex_no);
353 const Point<dim> &vertex_point = cell->vertex(vertex_no);
355 const typename std::map<unsigned int, Point<dim>>::const_iterator
356 map_iter = new_points.find(vertex_index);
358 if (map_iter != map_end)
359 for (
unsigned int i = 0; i < dim; ++i)
360 if (constraints[i].is_constrained(
361 cell->vertex_dof_index(vertex_no, 0)) ==
false)
363 constraints[i].add_constraint(
364 cell->vertex_dof_index(vertex_no, 0),
366 (solve_for_absolute_positions ?
367 map_iter->second(i) :
368 map_iter->second(i) - vertex_point[i]));
373 for (
unsigned int i = 0; i < dim; ++i)
374 constraints[i].close();
378 for (
unsigned int i = 0; i < dim; ++i)
383 for (
unsigned int i = 0; i < dim; ++i)
390 std::vector<bool> vertex_touched(
triangulation.n_vertices(),
false);
392 for (
const unsigned int vertex_no : cell->vertex_indices())
393 if (vertex_touched[cell->vertex_index(vertex_no)] ==
false)
398 cell->vertex_dof_index(vertex_no, 0);
399 for (
unsigned int i = 0; i < dim; ++i)
400 if (solve_for_absolute_positions)
401 v(i) = us[i](dof_index);
403 v(i) += us[i](dof_index);
405 vertex_touched[cell->vertex_index(vertex_no)] =
true;
413 template <
int dim,
int spacedim>
417 const bool keep_boundary,
418 const unsigned int seed)
435 double almost_infinite_length = 0;
440 almost_infinite_length += cell->diameter();
442 std::vector<double> minimal_length(
triangulation.n_vertices(),
443 almost_infinite_length);
446 std::vector<bool> at_boundary(keep_boundary ?
triangulation.n_vertices() :
451 const bool is_parallel_shared =
454 for (
const auto &cell :
triangulation.active_cell_iterators())
455 if (is_parallel_shared || cell->is_locally_owned())
459 for (
unsigned int i = 0; i < cell->n_lines(); ++i)
462 line = cell->line(i);
464 if (keep_boundary && line->at_boundary())
466 at_boundary[line->vertex_index(0)] =
true;
467 at_boundary[line->vertex_index(1)] =
true;
470 minimal_length[line->vertex_index(0)] =
472 minimal_length[line->vertex_index(0)]);
473 minimal_length[line->vertex_index(1)] =
475 minimal_length[line->vertex_index(1)]);
481 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
482 if (cell->at_boundary(vertex) ==
true)
483 at_boundary[cell->vertex_index(vertex)] =
true;
485 minimal_length[cell->vertex_index(0)] =
487 minimal_length[cell->vertex_index(0)]);
488 minimal_length[cell->vertex_index(1)] =
490 minimal_length[cell->vertex_index(1)]);
495 boost::random::mt19937 rng(seed);
496 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
500 if (
auto distributed_triangulation =
504 const std::vector<bool> locally_owned_vertices =
506 std::vector<bool> vertex_moved(
triangulation.n_vertices(),
false);
509 for (
const auto &cell :
triangulation.active_cell_iterators())
510 if (cell->is_locally_owned())
512 for (
const unsigned int vertex_no : cell->vertex_indices())
514 const unsigned global_vertex_no =
515 cell->vertex_index(vertex_no);
520 if ((keep_boundary && at_boundary[global_vertex_no]) ||
521 vertex_moved[global_vertex_no] ||
522 !locally_owned_vertices[global_vertex_no])
527 for (
unsigned int d = 0;
d < spacedim; ++
d)
528 shift_vector(
d) = uniform_distribution(rng);
530 shift_vector *= factor * minimal_length[global_vertex_no] /
531 std::sqrt(shift_vector.
square());
534 cell->vertex(vertex_no) += shift_vector;
535 vertex_moved[global_vertex_no] =
true;
539 distributed_triangulation->communicate_locally_moved_vertices(
540 locally_owned_vertices);
550 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
551 const std::vector<Point<spacedim>> &old_vertex_locations =
554 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
558 if (keep_boundary && at_boundary[vertex])
559 new_vertex_locations[vertex] = old_vertex_locations[vertex];
564 for (
unsigned int d = 0;
d < spacedim; ++
d)
565 shift_vector(
d) = uniform_distribution(rng);
567 shift_vector *= factor * minimal_length[vertex] /
568 std::sqrt(shift_vector.
square());
571 new_vertex_locations[vertex] =
572 old_vertex_locations[vertex] + shift_vector;
577 for (
const auto &cell :
triangulation.active_cell_iterators())
578 for (
const unsigned int vertex_no : cell->vertex_indices())
579 cell->vertex(vertex_no) =
580 new_vertex_locations[cell->vertex_index(vertex_no)];
596 for (; cell != endc; ++cell)
597 if (!cell->is_artificial())
598 for (
const unsigned int face : cell->face_indices())
599 if (cell->face(face)->has_children() &&
600 !cell->face(face)->at_boundary())
604 cell->face(face)->child(0)->vertex(1) =
605 (cell->face(face)->vertex(0) +
606 cell->face(face)->vertex(1)) /
610 cell->face(face)->child(0)->vertex(1) =
611 .5 * (cell->face(face)->vertex(0) +
612 cell->face(face)->vertex(1));
613 cell->face(face)->child(0)->vertex(2) =
614 .5 * (cell->face(face)->vertex(0) +
615 cell->face(face)->vertex(2));
616 cell->face(face)->child(1)->vertex(3) =
617 .5 * (cell->face(face)->vertex(1) +
618 cell->face(face)->vertex(3));
619 cell->face(face)->child(2)->vertex(3) =
620 .5 * (cell->face(face)->vertex(2) +
621 cell->face(face)->vertex(3));
624 cell->face(face)->child(0)->vertex(3) =
625 .25 * (cell->face(face)->vertex(0) +
626 cell->face(face)->vertex(1) +
627 cell->face(face)->vertex(2) +
628 cell->face(face)->vertex(3));
636 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
640 const
Point<spacedim> &p,
641 const std::vector<
bool> &marked_vertices)
650 marked_vertices.empty(),
652 marked_vertices.size()));
659 marked_vertices.empty() ||
660 std::equal(marked_vertices.begin(),
661 marked_vertices.end(),
663 [](
bool p,
bool q) { return !p || q; }),
665 "marked_vertices should be a subset of used vertices in the triangulation "
666 "but marked_vertices contains one or more vertices that are not used vertices!"));
670 const std::vector<bool> &vertices_to_use =
675 std::vector<bool>::const_iterator
first =
676 std::find(vertices_to_use.begin(), vertices_to_use.end(),
true);
681 unsigned int best_vertex = std::distance(vertices_to_use.begin(),
first);
682 double best_dist = (p -
vertices[best_vertex]).norm_square();
686 for (
unsigned int j = best_vertex + 1; j <
vertices.size(); ++j)
687 if (vertices_to_use[j])
689 const double dist = (p -
vertices[j]).norm_square();
690 if (dist < best_dist)
702 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
706 const MeshType<dim, spacedim> &mesh,
707 const
Point<spacedim> &p,
708 const std::vector<
bool> &marked_vertices)
711 if (mapping.preserves_vertex_locations() ==
true)
721 marked_vertices.empty(),
723 marked_vertices.size()));
731 marked_vertices.empty() ||
732 std::equal(marked_vertices.begin(),
733 marked_vertices.end(),
735 [](
bool p,
bool q) { return !p || q; }),
737 "marked_vertices should be a subset of used vertices in the triangulation "
738 "but marked_vertices contains one or more vertices that are not used vertices!"));
741 if (marked_vertices.size() != 0)
744 if (marked_vertices[it->first] ==
false)
759 template <
int dim,
int spacedim>
760 std::vector<std::vector<Tensor<1, spacedim>>>
768 const unsigned int n_vertices = vertex_to_cells.size();
773 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
775 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
778 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
779 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
781 typename std::set<typename Triangulation<dim, spacedim>::
782 active_cell_iterator>::iterator it =
783 vertex_to_cells[vertex].begin();
784 for (
unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
786 vertex_to_cell_centers[vertex][cell] =
788 vertex_to_cell_centers[vertex][cell] /=
789 vertex_to_cell_centers[vertex][cell].
norm();
792 return vertex_to_cell_centers;
798 template <
int spacedim>
801 const unsigned int a,
802 const unsigned int b,
806 const double scalar_product_a = center_directions[a] * point_direction;
807 const double scalar_product_b = center_directions[
b] * point_direction;
812 return (scalar_product_a > scalar_product_b);
816 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
820 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
822 std::pair<typename ::internal::
823 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
828 const MeshType<dim, spacedim> &mesh,
831 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
834 &vertex_to_cell_centers,
835 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
836 const std::vector<bool> &marked_vertices,
838 &used_vertices_rtree,
839 const double tolerance,
843 *relevant_cell_bounding_boxes_rtree)
845 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
848 cell_and_position.first = mesh.end();
852 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
854 cell_and_position_approx;
856 if (relevant_cell_bounding_boxes_rtree !=
nullptr &&
857 !relevant_cell_bounding_boxes_rtree->empty())
862 if (relevant_cell_bounding_boxes_rtree->qbegin(
863 boost::geometry::index::intersects(bb)) ==
864 relevant_cell_bounding_boxes_rtree->qend())
865 return cell_and_position;
868 bool found_cell =
false;
869 bool approx_cell =
false;
871 unsigned int closest_vertex_index = 0;
876 if (marked_vertices.size() && !used_vertices_rtree.empty())
879 std::find(marked_vertices.begin(), marked_vertices.end(),
true);
880 Assert(itr != marked_vertices.end(),
882 closest_vertex_index = std::distance(marked_vertices.begin(), itr);
886 auto current_cell = cell_hint;
889 const auto cell_marked = [&mesh, &marked_vertices](
const auto &cell) {
890 if (marked_vertices.empty())
893 if (cell != mesh.active_cell_iterators().end())
894 for (
unsigned int i = 0; i < cell->n_vertices(); ++i)
895 if (marked_vertices[cell->vertex_index(i)])
902 const auto any_cell_marked = [&cell_marked](
const auto &cells) {
903 return std::any_of(cells.begin(),
905 [&cell_marked](
const auto &cell) {
906 return cell_marked(cell);
909 (void)any_cell_marked;
911 while (found_cell ==
false)
916 cell_marked(cell_hint))
918 const auto cell_vertices = mapping.
get_vertices(current_cell);
919 const unsigned int closest_vertex =
920 find_closest_vertex_of_cell<dim, spacedim>(current_cell,
923 vertex_to_point = p - cell_vertices[closest_vertex];
924 closest_vertex_index = current_cell->vertex_index(closest_vertex);
933 #if defined(DEAL_II_WITH_BOOST_BUNDLED) || \
934 !(defined(__clang_major__) && __clang_major__ >= 16) || \
935 BOOST_VERSION >= 108100
936 if (!used_vertices_rtree.empty())
939 using ValueType = std::pair<Point<spacedim>,
unsigned int>;
940 std::function<
bool(
const ValueType &)> marked;
941 if (marked_vertices.size() == mesh.n_vertices())
942 marked = [&marked_vertices](
const ValueType &value) ->
bool {
943 return marked_vertices[value.second];
946 marked = [](
const ValueType &) ->
bool {
return true; };
948 std::vector<std::pair<Point<spacedim>,
unsigned int>> res;
949 used_vertices_rtree.query(
950 boost::geometry::index::nearest(p, 1) &&
951 boost::geometry::index::satisfies(marked),
952 std::back_inserter(res));
957 ::
ExcMessage(
"There can not be multiple results"));
960 if (any_cell_marked(vertex_to_cells[res[0].
second]))
961 closest_vertex_index = res[0].second;
967 mapping, mesh, p, marked_vertices);
969 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
975 Assert(any_cell_marked(vertex_to_cells[closest_vertex_index]),
980 const double vertex_point_norm = vertex_to_point.
norm();
981 if (vertex_point_norm > 0)
982 vertex_to_point /= vertex_point_norm;
984 const unsigned int n_neighbor_cells =
985 vertex_to_cells[closest_vertex_index].size();
988 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
990 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
991 neighbor_permutation[i] = i;
993 auto comp = [&](
const unsigned int a,
const unsigned int b) ->
bool {
994 return internal::compare_point_association<spacedim>(
998 vertex_to_cell_centers[closest_vertex_index]);
1001 std::sort(neighbor_permutation.begin(),
1002 neighbor_permutation.end(),
1007 double best_distance = tolerance;
1011 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
1015 auto cell = vertex_to_cells[closest_vertex_index].begin();
1018 if (!(*cell)->is_artificial())
1022 if ((*cell)->reference_cell().contains_point(p_unit,
1025 cell_and_position.first = *cell;
1026 cell_and_position.second = p_unit;
1028 approx_cell =
false;
1035 const double dist = p_unit.
distance(
1036 (*cell)->reference_cell().closest_point(p_unit));
1037 if (dist < best_distance)
1039 best_distance = dist;
1040 cell_and_position_approx.first = *cell;
1041 cell_and_position_approx.second = p_unit;
1050 if (found_cell ==
true)
1051 return cell_and_position;
1052 else if (approx_cell ==
true)
1053 return cell_and_position_approx;
1069 mapping, mesh, p, marked_vertices, tolerance);
1071 current_cell =
typename MeshType<dim, spacedim>::active_cell_iterator();
1073 return cell_and_position;
1078 template <
int dim,
int spacedim>
1087 unsigned int closest_vertex = 0;
1089 for (
unsigned int v = 1; v < cell->
n_vertices(); ++v)
1092 if (vertex_distance < minimum_distance)
1095 minimum_distance = vertex_distance;
1098 return closest_vertex;
1105 namespace BoundingBoxPredicate
1107 template <
typename MeshType>
1109 concepts::is_triangulation_or_dof_handler<MeshType>)
1113 cell_iterator &parent_cell,
1114 const std::function<
bool(
1115 const typename MeshType::
1116 active_cell_iterator &)>
1119 bool has_predicate =
1121 std::vector<typename MeshType::active_cell_iterator> active_cells;
1122 if (parent_cell->is_active())
1123 active_cells = {parent_cell};
1127 active_cells = get_active_child_cells<MeshType>(parent_cell);
1129 const unsigned int spacedim = MeshType::space_dimension;
1133 while (i < active_cells.size() && !predicate(active_cells[i]))
1137 if (active_cells.empty() || i == active_cells.size())
1140 return std::make_tuple(bbox, has_predicate);
1147 for (; i < active_cells.size(); ++i)
1148 if (predicate(active_cells[i]))
1150 for (
unsigned int d = 0;
d < spacedim; ++
d)
1152 minp[
d] =
std::min(minp[
d], active_cells[i]->vertex(v)[
d]);
1153 maxp[
d] =
std::max(maxp[
d], active_cells[i]->vertex(v)[
d]);
1156 has_predicate =
true;
1158 return std::make_tuple(bbox, has_predicate);
1165 template <
typename MeshType>
1169 const MeshType &mesh,
1170 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1172 const unsigned int refinement_level,
1173 const bool allow_merge,
1174 const unsigned int max_boxes)
1181 refinement_level <= mesh.n_levels(),
1183 "Error: refinement level is higher then total levels in the triangulation!"));
1185 const unsigned int spacedim = MeshType::space_dimension;
1186 std::vector<BoundingBox<spacedim>> bounding_boxes;
1190 for (
unsigned int i = 0; i < refinement_level; ++i)
1191 for (
const typename MeshType::cell_iterator &cell :
1192 mesh.active_cell_iterators_on_level(i))
1194 bool has_predicate =
false;
1196 std::tie(bbox, has_predicate) =
1198 MeshType>(cell, predicate);
1200 bounding_boxes.push_back(bbox);
1204 for (
const typename MeshType::cell_iterator &cell :
1205 mesh.cell_iterators_on_level(refinement_level))
1207 bool has_predicate =
false;
1209 std::tie(bbox, has_predicate) =
1211 MeshType>(cell, predicate);
1213 bounding_boxes.push_back(bbox);
1218 return bounding_boxes;
1224 std::vector<unsigned int> merged_boxes_idx;
1225 bool found_neighbors =
true;
1230 while (found_neighbors)
1232 found_neighbors =
false;
1233 for (
unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
1235 if (std::find(merged_boxes_idx.begin(),
1236 merged_boxes_idx.end(),
1237 i) == merged_boxes_idx.end())
1238 for (
unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
1239 if (std::find(merged_boxes_idx.begin(),
1240 merged_boxes_idx.end(),
1241 j) == merged_boxes_idx.end() &&
1242 bounding_boxes[i].get_neighbor_type(
1243 bounding_boxes[j]) ==
1246 bounding_boxes[i].merge_with(bounding_boxes[j]);
1247 merged_boxes_idx.push_back(j);
1248 found_neighbors =
true;
1254 std::vector<BoundingBox<spacedim>> merged_b_boxes;
1255 for (
unsigned int i = 0; i < bounding_boxes.size(); ++i)
1256 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
1257 merged_boxes_idx.end())
1258 merged_b_boxes.push_back(bounding_boxes[i]);
1263 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
1265 std::vector<double> volumes;
1266 for (
unsigned int i = 0; i < merged_b_boxes.size(); ++i)
1267 volumes.push_back(merged_b_boxes[i].volume());
1269 while (merged_b_boxes.size() > max_boxes)
1271 unsigned int min_idx =
1272 std::min_element(volumes.begin(), volumes.end()) -
1274 volumes.erase(volumes.begin() + min_idx);
1276 bool not_removed =
true;
1277 for (
unsigned int i = 0;
1278 i < merged_b_boxes.size() && not_removed;
1283 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
1284 merged_b_boxes[min_idx]) ==
1286 merged_b_boxes[i].get_neighbor_type(
1287 merged_b_boxes[min_idx]) ==
1290 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
1291 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
1292 not_removed =
false;
1295 ExcMessage(
"Error: couldn't merge bounding boxes!"));
1298 Assert(merged_b_boxes.size() <= max_boxes,
1300 "Error: couldn't reach target number of bounding boxes!"));
1301 return merged_b_boxes;
1307 template <
int spacedim>
1309 std::tuple<std::vector<std::vector<unsigned int>>,
1310 std::map<unsigned int, unsigned int>,
1311 std::map<unsigned int, std::vector<unsigned int>>>
1319 unsigned int n_procs = global_bboxes.size();
1320 std::vector<std::vector<unsigned int>> point_owners(n_procs);
1321 std::map<unsigned int, unsigned int> map_owners_found;
1322 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
1324 unsigned int n_points = points.size();
1325 for (
unsigned int pt = 0; pt < n_points; ++pt)
1328 std::vector<unsigned int> owners_found;
1330 for (
unsigned int rk = 0; rk < n_procs; ++rk)
1333 if (bbox.point_inside(points[pt]))
1335 point_owners[rk].emplace_back(pt);
1336 owners_found.emplace_back(rk);
1340 Assert(owners_found.size() > 0,
1341 ExcMessage(
"No owners found for the point " +
1343 if (owners_found.size() == 1)
1344 map_owners_found[pt] = owners_found[0];
1347 map_owners_guessed[pt] = owners_found;
1350 return std::make_tuple(std::move(point_owners),
1351 std::move(map_owners_found),
1352 std::move(map_owners_guessed));
1355 template <
int spacedim>
1357 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
1358 std::map<unsigned int, unsigned int>,
1359 std::map<unsigned int, std::vector<unsigned int>>>
1367 std::map<unsigned int, std::vector<unsigned int>> point_owners;
1368 std::map<unsigned int, unsigned int> map_owners_found;
1369 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
1370 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> search_result;
1372 unsigned int n_points = points.size();
1373 for (
unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
1375 search_result.clear();
1378 covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
1379 std::back_inserter(search_result));
1382 std::set<unsigned int> owners_found;
1384 for (
const auto &rank_bbox : search_result)
1388 const bool pt_inserted = owners_found.insert(pt_n).second;
1390 point_owners[rank_bbox.second].emplace_back(pt_n);
1392 Assert(owners_found.size() > 0,
1393 ExcMessage(
"No owners found for the point " +
1395 if (owners_found.size() == 1)
1396 map_owners_found[pt_n] = *owners_found.begin();
1401 std::back_inserter(map_owners_guessed[pt_n]));
1404 return std::make_tuple(std::move(point_owners),
1405 std::move(map_owners_found),
1406 std::move(map_owners_guessed));
1411 template <
int dim,
int spacedim>
1412 std::map<unsigned int, types::global_vertex_index>
1416 std::map<unsigned int, types::global_vertex_index>
1417 local_to_global_vertex_index;
1419 #ifndef DEAL_II_WITH_MPI
1426 ExcMessage(
"This function does not make any sense "
1427 "for parallel::distributed::Triangulation "
1428 "objects if you do not have MPI enabled."));
1432 using active_cell_iterator =
1434 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
1439 unsigned int max_cellid_size = 0;
1440 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
1442 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
1448 std::set<active_cell_iterator> missing_vert_cells;
1449 std::set<unsigned int> used_vertex_index;
1450 for (
const auto &cell :
triangulation.active_cell_iterators())
1452 if (cell->is_locally_owned())
1454 for (
const unsigned int i : cell->vertex_indices())
1457 for (
const auto &adjacent_cell :
1458 vertex_to_cell[cell->vertex_index(i)])
1459 lowest_subdomain_id =
std::min(lowest_subdomain_id,
1460 adjacent_cell->subdomain_id());
1463 if (lowest_subdomain_id == cell->subdomain_id())
1467 if (used_vertex_index.find(cell->vertex_index(i)) ==
1468 used_vertex_index.end())
1471 local_to_global_vertex_index[cell->vertex_index(i)] =
1476 for (
const auto &adjacent_cell :
1477 vertex_to_cell[cell->vertex_index(i)])
1478 if (adjacent_cell->subdomain_id() !=
1479 cell->subdomain_id())
1483 tmp(adjacent_cell->subdomain_id(),
1484 cell->vertex_index(i));
1485 if (vertices_added.find(tmp) ==
1486 vertices_added.end())
1488 vertices_to_send[adjacent_cell
1491 cell->vertex_index(i),
1492 cell->id().to_string());
1493 if (cell->id().to_string().size() >
1496 cell->id().to_string().size();
1497 vertices_added.insert(tmp);
1500 used_vertex_index.insert(cell->vertex_index(i));
1507 vertices_to_recv[lowest_subdomain_id].insert(
1508 cell->vertex_index(i));
1509 missing_vert_cells.insert(cell);
1516 if (cell->is_ghost())
1518 for (
const unsigned int i : cell->face_indices())
1520 if (cell->at_boundary(i) ==
false)
1522 if (cell->neighbor(i)->is_active())
1525 spacedim>::active_cell_iterator
1526 adjacent_cell = cell->neighbor(i);
1527 if ((adjacent_cell->is_locally_owned()))
1530 adjacent_cell->subdomain_id();
1531 if (cell->subdomain_id() < adj_subdomain_id)
1532 for (
unsigned int j = 0;
1533 j < cell->face(i)->n_vertices();
1536 vertices_to_recv[cell->subdomain_id()].insert(
1537 cell->face(i)->vertex_index(j));
1538 missing_vert_cells.insert(cell);
1554 int ierr = MPI_Exscan(&next_index,
1562 for (
auto &global_index_it : local_to_global_vertex_index)
1563 global_index_it.second +=
shift;
1577 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
1578 vertices_to_send.size());
1579 std::vector<MPI_Request> first_requests(vertices_to_send.size());
1583 std::string>>>::iterator
1584 vert_to_send_it = vertices_to_send.begin(),
1585 vert_to_send_end = vertices_to_send.end();
1586 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
1587 ++vert_to_send_it, ++i)
1589 int destination = vert_to_send_it->first;
1590 const unsigned int n_vertices = vert_to_send_it->second.size();
1591 const int buffer_size = 2 * n_vertices;
1592 vertices_send_buffers[i].resize(buffer_size);
1595 for (
unsigned int j = 0; j < n_vertices; ++j)
1597 vertices_send_buffers[i][2 * j] =
1598 std::get<0>(vert_to_send_it->second[j]);
1599 vertices_send_buffers[i][2 * j + 1] =
1600 local_to_global_vertex_index[std::get<1>(
1601 vert_to_send_it->second[j])];
1605 ierr = MPI_Isend(vertices_send_buffers[i].data(),
1611 &first_requests[i]);
1616 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
1617 vertices_to_recv.size());
1618 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
1619 vert_to_recv_it = vertices_to_recv.begin(),
1620 vert_to_recv_end = vertices_to_recv.end();
1621 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1622 ++vert_to_recv_it, ++i)
1624 int source = vert_to_recv_it->first;
1625 const unsigned int n_vertices = vert_to_recv_it->second.size();
1626 const int buffer_size = 2 * n_vertices;
1627 vertices_recv_buffers[i].resize(buffer_size);
1630 ierr = MPI_Recv(vertices_recv_buffers[i].data(),
1641 MPI_Waitall(first_requests.size(),
1642 first_requests.data(),
1643 MPI_STATUSES_IGNORE);
1647 std::vector<std::vector<char>> cellids_send_buffers(
1648 vertices_to_send.size());
1649 std::vector<MPI_Request> second_requests(vertices_to_send.size());
1650 vert_to_send_it = vertices_to_send.begin();
1651 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
1652 ++vert_to_send_it, ++i)
1654 int destination = vert_to_send_it->first;
1655 const unsigned int n_vertices = vert_to_send_it->second.size();
1656 const int buffer_size = max_cellid_size * n_vertices;
1657 cellids_send_buffers[i].resize(buffer_size);
1660 unsigned int pos = 0;
1661 for (
unsigned int j = 0; j < n_vertices; ++j)
1663 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
1664 for (
unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
1666 if (k < cell_id.size())
1667 cellids_send_buffers[i][pos] = cell_id[k];
1671 cellids_send_buffers[i][pos] =
'-';
1676 ierr = MPI_Isend(cellids_send_buffers[i].data(),
1682 &second_requests[i]);
1687 std::vector<std::vector<char>> cellids_recv_buffers(
1688 vertices_to_recv.size());
1689 vert_to_recv_it = vertices_to_recv.begin();
1690 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1691 ++vert_to_recv_it, ++i)
1693 int source = vert_to_recv_it->first;
1694 const unsigned int n_vertices = vert_to_recv_it->second.size();
1695 const int buffer_size = max_cellid_size * n_vertices;
1696 cellids_recv_buffers[i].resize(buffer_size);
1699 ierr = MPI_Recv(cellids_recv_buffers[i].data(),
1711 vert_to_recv_it = vertices_to_recv.begin();
1712 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
1713 ++i, ++vert_to_recv_it)
1715 for (
unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
1717 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
1719 vertices_recv_buffers[i][2 * j + 1];
1720 const std::string cellid_recv(
1721 &cellids_recv_buffers[i][max_cellid_size * j],
1722 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
1724 typename std::set<active_cell_iterator>::iterator
1725 cell_set_it = missing_vert_cells.begin(),
1726 end_cell_set = missing_vert_cells.end();
1727 for (; (found ==
false) && (cell_set_it != end_cell_set);
1730 typename std::set<active_cell_iterator>::iterator
1732 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
1734 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
1735 for (; candidate_cell != end_cell; ++candidate_cell)
1737 std::string current_cellid =
1738 (*candidate_cell)->id().to_string();
1739 current_cellid.resize(max_cellid_size,
'-');
1740 if (current_cellid.compare(cellid_recv) == 0)
1742 local_to_global_vertex_index
1743 [(*candidate_cell)->vertex_index(local_pos_recv)] =
1756 MPI_Waitall(second_requests.size(),
1757 second_requests.data(),
1758 MPI_STATUSES_IGNORE);
1761 return local_to_global_vertex_index;
1766 template <
int dim,
int spacedim>
1774 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
1775 "are already partitioned implicitly and can not be "
1776 "partitioned again explicitly."));
1778 std::vector<unsigned int> cell_weights;
1788 for (
const auto &cell :
triangulation.active_cell_iterators())
1789 if (cell->is_locally_owned())
1790 cell_weights[cell->active_cell_index()] =
1799 if (
const auto shared_tria =
1803 shared_tria->get_communicator(),
1807 Assert(std::accumulate(cell_weights.begin(),
1809 std::uint64_t(0)) > 0,
1810 ExcMessage(
"The global sum of weights over all active cells "
1811 "is zero. Please verify how you generate weights."));
1823 template <
int dim,
int spacedim>
1826 const std::vector<unsigned int> &cell_weights,
1832 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
1833 "are already partitioned implicitly and can not be "
1834 "partitioned again explicitly."));
1838 if (n_partitions == 1)
1840 for (
const auto &cell :
triangulation.active_cell_iterators())
1841 cell->set_subdomain_id(0);
1855 sp_cell_connectivity.
copy_from(cell_connectivity);
1858 sp_cell_connectivity,
1865 template <
int dim,
int spacedim>
1874 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
1875 "are already partitioned implicitly and can not be "
1876 "partitioned again explicitly."));
1878 std::vector<unsigned int> cell_weights;
1888 for (
const auto &cell :
triangulation.active_cell_iterators() |
1890 cell_weights[cell->active_cell_index()] =
1899 if (
const auto shared_tria =
1903 shared_tria->get_communicator(),
1907 Assert(std::accumulate(cell_weights.begin(),
1909 std::uint64_t(0)) > 0,
1910 ExcMessage(
"The global sum of weights over all active cells "
1911 "is zero. Please verify how you generate weights."));
1917 cell_connection_graph,
1924 template <
int dim,
int spacedim>
1927 const std::vector<unsigned int> &cell_weights,
1934 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
1935 "are already partitioned implicitly and can not be "
1936 "partitioned again explicitly."));
1939 ExcMessage(
"Connectivity graph has wrong size"));
1941 ExcMessage(
"Connectivity graph has wrong size"));
1947 if (n_partitions == 1)
1949 for (
const auto &cell :
triangulation.active_cell_iterators())
1950 cell->set_subdomain_id(0);
1958 std::vector<unsigned int> partition_indices(
triangulation.n_active_cells());
1966 for (
const auto &cell :
triangulation.active_cell_iterators())
1967 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
1979 unsigned int ¤t_proc_idx,
1980 unsigned int ¤t_cell_idx,
1982 const unsigned int n_partitions)
1984 if (cell->is_active())
1986 while (current_cell_idx >=
1988 (current_proc_idx + 1) / n_partitions))
1990 cell->set_subdomain_id(current_proc_idx);
1995 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2005 template <
int dim,
int spacedim>
2009 const bool group_siblings)
2013 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
2014 "are already partitioned implicitly and can not be "
2015 "partitioned again explicitly."));
2023 if (n_partitions == 1)
2025 for (
const auto &cell :
triangulation.active_cell_iterators())
2026 cell->set_subdomain_id(0);
2032 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
2033 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
2039 coarse_cell_to_p4est_tree_permutation.resize(
triangulation.n_cells(0));
2041 coarse_cell_to_p4est_tree_permutation);
2043 p4est_tree_to_coarse_cell_permutation =
2046 unsigned int current_proc_idx = 0;
2047 unsigned int current_cell_idx = 0;
2052 for (
unsigned int idx = 0; idx <
triangulation.n_cells(0); ++idx)
2054 const unsigned int coarse_cell_idx =
2055 p4est_tree_to_coarse_cell_permutation[idx];
2078 for (; cell != endc; ++cell)
2080 if (cell->is_active())
2082 bool all_children_active =
true;
2083 std::map<unsigned int, unsigned int> map_cpu_n_cells;
2084 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2085 if (!cell->child(n)->is_active())
2087 all_children_active =
false;
2091 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
2093 if (!all_children_active)
2096 unsigned int new_owner = cell->child(0)->subdomain_id();
2097 for (std::map<unsigned int, unsigned int>::iterator it =
2098 map_cpu_n_cells.begin();
2099 it != map_cpu_n_cells.end();
2101 if (it->second > map_cpu_n_cells[new_owner])
2102 new_owner = it->first;
2104 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2105 cell->child(n)->set_subdomain_id(new_owner);
2111 template <
int dim,
int spacedim>
2116 for (
int lvl = n_levels - 1; lvl >= 0; --lvl)
2118 for (
const auto &cell :
triangulation.cell_iterators_on_level(lvl))
2120 if (cell->is_active())
2121 cell->set_level_subdomain_id(cell->subdomain_id());
2124 Assert(cell->child(0)->level_subdomain_id() !=
2127 cell->set_level_subdomain_id(
2128 cell->child(0)->level_subdomain_id());
2140 template <
int dim,
int spacedim>
2145 const std::vector<CellId> &cell_ids,
2146 std::vector<types::subdomain_id> &subdomain_ids)
2148 #ifndef DEAL_II_WITH_P4EST
2151 (void)subdomain_ids;
2155 "You are attempting to use a functionality that is only available "
2156 "if deal.II was configured to use p4est, but cmake did not find a "
2157 "valid p4est library."));
2162 for (
const auto &cell_id : cell_ids)
2165 typename ::internal::p4est::types<dim>::quadrant p4est_cell,
2168 ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
2169 for (
const auto &child_index : cell_id.get_child_indices())
2171 ::internal::p4est::init_quadrant_children<dim>(
2172 p4est_cell, p4est_children);
2174 p4est_children[
static_cast<unsigned int>(child_index)];
2180 const_cast<typename ::internal::p4est::types<dim>::forest
2182 cell_id.get_coarse_cell_id(),
2189 subdomain_ids.push_back(owner);
2196 template <
int spacedim>
2200 const std::vector<CellId> &,
2201 std::vector<types::subdomain_id> &)
2210 template <
int dim,
int spacedim>
2211 std::vector<types::subdomain_id>
2213 const std::vector<CellId> &cell_ids)
2215 std::vector<types::subdomain_id> subdomain_ids;
2216 subdomain_ids.reserve(cell_ids.size());
2225 *parallel_tria =
dynamic_cast<
2239 const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
2240 shared_tria->get_true_subdomain_ids_of_cells();
2242 for (
const auto &cell_id : cell_ids)
2244 const unsigned int active_cell_index =
2245 shared_tria->create_cell_iterator(cell_id)->active_cell_index();
2246 subdomain_ids.push_back(
2247 true_subdomain_ids_of_cells[active_cell_index]);
2254 for (
const auto &cell_id : cell_ids)
2256 subdomain_ids.push_back(
2257 triangulation.create_cell_iterator(cell_id)->subdomain_id());
2261 return subdomain_ids;
2266 template <
int dim,
int spacedim>
2269 std::vector<types::subdomain_id> &subdomain)
2274 for (
const auto &cell :
triangulation.active_cell_iterators())
2275 subdomain[cell->active_cell_index()] = cell->subdomain_id();
2280 template <
int dim,
int spacedim>
2286 unsigned int count = 0;
2287 for (
const auto &cell :
triangulation.active_cell_iterators())
2288 if (cell->subdomain_id() == subdomain)
2296 template <
int dim,
int spacedim>
2301 std::vector<bool> locally_owned_vertices =
2308 if (
const auto *tr =
dynamic_cast<
2311 for (
const auto &cell :
triangulation.active_cell_iterators())
2312 if (cell->is_artificial() ||
2313 (cell->is_ghost() &&
2314 (cell->subdomain_id() < tr->locally_owned_subdomain())))
2315 for (
const unsigned int v : cell->vertex_indices())
2316 locally_owned_vertices[cell->vertex_index(v)] =
false;
2318 return locally_owned_vertices;
2325 namespace FixUpDistortedChildCells
2348 template <
typename Iterator,
int spacedim>
2353 const unsigned int structdim =
2354 Iterator::AccessorType::structure_dimension;
2355 Assert(spacedim == Iterator::AccessorType::dimension,
2361 Assert(object->refinement_case() ==
2369 Tensor<spacedim - structdim, spacedim>
2372 for (
const unsigned int i : object->vertex_indices())
2373 parent_vertices[i] =
object->vertex(i);
2376 parent_vertices, parent_alternating_forms);
2378 const Tensor<spacedim - structdim, spacedim>
2379 average_parent_alternating_form =
2380 std::accumulate(parent_alternating_forms,
2381 parent_alternating_forms +
2395 Tensor<spacedim - structdim, spacedim> child_alternating_forms
2399 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2400 for (
const unsigned int i : object->child(c)->vertex_indices())
2401 child_vertices[c][i] =
object->child(c)->vertex(i);
2409 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2411 1] = object_mid_point;
2413 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2415 child_vertices[c], child_alternating_forms[c]);
2427 double objective = 0;
2428 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2429 for (
const unsigned int i : object->child(c)->vertex_indices())
2431 (child_alternating_forms[c][i] -
2432 average_parent_alternating_form / std::pow(2., 1. * structdim))
2444 template <
typename Iterator>
2447 const unsigned int f,
2448 std::integral_constant<int, 1>)
2450 return object->vertex(f);
2460 template <
typename Iterator>
2463 const unsigned int f,
2464 std::integral_constant<int, 2>)
2466 return object->line(f)->center();
2476 template <
typename Iterator>
2479 const unsigned int f,
2480 std::integral_constant<int, 3>)
2482 return object->face(f)->center();
2509 template <
typename Iterator>
2513 const unsigned int structdim =
2514 Iterator::AccessorType::structure_dimension;
2516 double diameter =
object->diameter();
2517 for (
const unsigned int f : object->face_indices())
2518 for (
unsigned int e = f + 1;
e <
object->n_faces(); ++
e)
2523 std::integral_constant<int, structdim>())
2525 object,
e, std::integral_constant<int, structdim>())));
2536 template <
typename Iterator>
2540 const unsigned int structdim =
2541 Iterator::AccessorType::structure_dimension;
2542 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
2549 Assert(object->refinement_case() ==
2559 unsigned int iteration = 0;
2564 double initial_delta = 0;
2571 const double step_length =
diameter / 4 / (iteration + 1);
2576 for (
unsigned int d = 0;
d < spacedim; ++
d)
2578 const double eps = step_length / 10;
2592 if (gradient.norm() == 0)
2601 std::min(2 * current_value / (gradient * gradient),
2602 step_length / gradient.norm()) *
2607 const double previous_value = current_value;
2611 initial_delta = (previous_value - current_value);
2614 if ((iteration >= 1) &&
2615 ((previous_value - current_value < 0) ||
2616 (
std::fabs(previous_value - current_value) <
2617 0.001 * initial_delta)))
2622 while (iteration < 20);
2638 double old_min_product, new_min_product;
2643 parent_vertices[i] = object->vertex(i);
2645 Tensor<spacedim - structdim, spacedim>
2648 parent_vertices, parent_alternating_forms);
2654 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2655 for (
const unsigned int i : object->child(c)->vertex_indices())
2656 child_vertices[c][i] =
object->child(c)->vertex(i);
2658 Tensor<spacedim - structdim, spacedim> child_alternating_forms
2662 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2664 child_vertices[c], child_alternating_forms[c]);
2667 child_alternating_forms[0][0] * parent_alternating_forms[0];
2668 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2669 for (
const unsigned int i : object->child(c)->vertex_indices())
2670 for (
const unsigned int j : object->vertex_indices())
2671 old_min_product = std::min<double>(old_min_product,
2672 child_alternating_forms[c][i] *
2673 parent_alternating_forms[j]);
2681 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2683 1] = object_mid_point;
2685 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2687 child_vertices[c], child_alternating_forms[c]);
2690 child_alternating_forms[0][0] * parent_alternating_forms[0];
2691 for (
unsigned int c = 0; c <
object->n_children(); ++c)
2692 for (
const unsigned int i : object->child(c)->vertex_indices())
2693 for (
const unsigned int j : object->vertex_indices())
2694 new_min_product = std::min<double>(new_min_product,
2695 child_alternating_forms[c][i] *
2696 parent_alternating_forms[j]);
2704 if (new_min_product >= old_min_product)
2705 object->child(0)->vertex(
2712 return (
std::max(new_min_product, old_min_product) > 0);
2718 template <
int dim,
int spacedim>
2721 const typename ::Triangulation<dim, spacedim>::cell_iterator
2723 std::integral_constant<int, dim>,
2724 std::integral_constant<int, spacedim>)
2734 for (
auto f : cell->face_indices())
2737 Assert(cell->face(f)->refinement_case() ==
2741 bool subface_is_more_refined =
false;
2742 for (
unsigned int g = 0;
2743 g < GeometryInfo<dim>::max_children_per_face;
2745 if (cell->face(f)->child(g)->has_children())
2747 subface_is_more_refined =
true;
2751 if (subface_is_more_refined ==
true)
2762 template <
int dim,
int spacedim>
2770 dim != 1 && spacedim != 1,
2771 "This function is only valid when dim != 1 or spacedim != 1.");
2775 for (
typename std::list<
2777 cell_ptr = distorted_cells.distorted_cells.
begin();
2778 cell_ptr != distorted_cells.distorted_cells.
end();
2784 Assert(!cell->is_active(),
2786 "This function is only valid for a list of cells that "
2787 "have children (i.e., no cell in the list may be active)."));
2791 std::integral_constant<int, dim>(),
2792 std::integral_constant<int, spacedim>());
2796 unfixable_subset.distorted_cells.push_back(cell);
2799 return unfixable_subset;
2804 template <
int dim,
int spacedim>
2807 const bool reset_boundary_ids)
2810 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
2811 auto m_it = dst_manifold_ids.begin();
2812 for (
const auto b : src_boundary_ids)
2817 const std::vector<types::boundary_id> reset_boundary_id =
2818 reset_boundary_ids ?
2819 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
2829 template <
int dim,
int spacedim>
2832 const std::vector<types::boundary_id> &src_boundary_ids,
2833 const std::vector<types::manifold_id> &dst_manifold_ids,
2835 const std::vector<types::boundary_id> &reset_boundary_ids_)
2838 const auto reset_boundary_ids =
2839 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
2848 for (
auto f : cell->face_indices())
2849 if (cell->face(f)->at_boundary())
2850 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
2852 const auto bid = cell->face(f)->line(
e)->boundary_id();
2853 const unsigned int ind = std::find(src_boundary_ids.begin(),
2854 src_boundary_ids.end(),
2856 src_boundary_ids.begin();
2857 if (ind < src_boundary_ids.size())
2858 cell->face(f)->line(
e)->set_manifold_id(
2859 dst_manifold_ids[ind]);
2864 for (
auto f : cell->face_indices())
2865 if (cell->face(f)->at_boundary())
2867 const auto bid = cell->face(f)->boundary_id();
2868 const unsigned int ind =
2869 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
2870 src_boundary_ids.begin();
2872 if (ind < src_boundary_ids.size())
2875 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
2877 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
2881 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
2883 const auto bid = cell->face(f)->line(
e)->boundary_id();
2884 const unsigned int ind = std::find(src_boundary_ids.begin(),
2885 src_boundary_ids.end(),
2887 src_boundary_ids.begin();
2888 if (ind < src_boundary_ids.size())
2889 cell->face(f)->line(
e)->set_boundary_id(
2890 reset_boundary_ids[ind]);
2896 template <
int dim,
int spacedim>
2899 const bool compute_face_ids)
2905 for (; cell != endc; ++cell)
2907 cell->set_manifold_id(cell->material_id());
2908 if (compute_face_ids ==
true)
2910 for (
auto f : cell->face_indices())
2912 if (cell->at_boundary(f) ==
false)
2913 cell->face(f)->set_manifold_id(
2915 cell->neighbor(f)->material_id()));
2917 cell->face(f)->set_manifold_id(cell->material_id());
2924 template <
int dim,
int spacedim>
2929 const std::set<types::manifold_id> &)> &disambiguation_function,
2930 bool overwrite_only_flat_manifold_ids)
2935 const unsigned int n_subobjects =
2939 std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
2940 std::vector<unsigned int> backup;
2944 unsigned next_index = 1;
2948 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
2950 if (cell->line(
l)->user_index() == 0)
2953 manifold_ids[next_index].insert(cell->manifold_id());
2954 cell->line(
l)->set_user_index(next_index++);
2957 manifold_ids[cell->line(
l)->user_index()].insert(
2958 cell->manifold_id());
2961 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
2963 if (cell->quad(
l)->user_index() == 0)
2966 manifold_ids[next_index].insert(cell->manifold_id());
2967 cell->quad(
l)->set_user_index(next_index++);
2970 manifold_ids[cell->quad(
l)->user_index()].insert(
2971 cell->manifold_id());
2977 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
2979 const auto id = cell->line(
l)->user_index();
2983 if (cell->line(
l)->manifold_id() ==
2985 overwrite_only_flat_manifold_ids ==
false)
2986 cell->line(
l)->set_manifold_id(
2987 disambiguation_function(manifold_ids[
id]));
2988 cell->line(
l)->set_user_index(0);
2992 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
2994 const auto id = cell->quad(
l)->user_index();
2998 if (cell->quad(
l)->manifold_id() ==
3000 overwrite_only_flat_manifold_ids ==
false)
3001 cell->quad(
l)->set_manifold_id(
3002 disambiguation_function(manifold_ids[
id]));
3003 cell->quad(
l)->set_user_index(0);
3012 template <
int dim,
int spacedim>
3015 const double limit_angle_fraction)
3023 "have hanging nodes."));
3027 bool has_cells_with_more_than_dim_faces_on_boundary =
true;
3028 bool has_cells_with_dim_faces_on_boundary =
false;
3030 unsigned int refinement_cycles = 0;
3032 while (has_cells_with_more_than_dim_faces_on_boundary)
3034 has_cells_with_more_than_dim_faces_on_boundary =
false;
3038 unsigned int boundary_face_counter = 0;
3039 for (
auto f : cell->face_indices())
3040 if (cell->face(f)->at_boundary())
3041 boundary_face_counter++;
3042 if (boundary_face_counter > dim)
3044 has_cells_with_more_than_dim_faces_on_boundary =
true;
3047 else if (boundary_face_counter == dim)
3048 has_cells_with_dim_faces_on_boundary =
true;
3050 if (has_cells_with_more_than_dim_faces_on_boundary)
3053 refinement_cycles++;
3057 if (has_cells_with_dim_faces_on_boundary)
3060 refinement_cycles++;
3064 while (refinement_cycles > 0)
3067 cell->set_coarsen_flag();
3069 refinement_cycles--;
3079 std::vector<CellData<dim>> cells_to_add;
3083 const unsigned int v0 = 0,
v1 = 1, v2 = (dim > 1 ? 2 : 0),
3084 v3 = (dim > 1 ? 3 : 0);
3088 double angle_fraction = 0;
3094 p0[spacedim > 1 ? 1 : 0] = 1;
3098 if (cell->face(
v0)->at_boundary() && cell->face(v3)->at_boundary())
3100 p0 = cell->vertex(
v0) - cell->vertex(v2);
3101 p1 = cell->vertex(v3) - cell->vertex(v2);
3102 vertex_at_corner = v2;
3104 else if (cell->face(v3)->at_boundary() &&
3105 cell->face(
v1)->at_boundary())
3107 p0 = cell->vertex(v2) - cell->vertex(v3);
3108 p1 = cell->vertex(
v1) - cell->vertex(v3);
3109 vertex_at_corner = v3;
3111 else if (cell->face(1)->at_boundary() &&
3112 cell->face(2)->at_boundary())
3114 p0 = cell->vertex(
v0) - cell->vertex(
v1);
3115 p1 = cell->vertex(v3) - cell->vertex(
v1);
3116 vertex_at_corner =
v1;
3118 else if (cell->face(2)->at_boundary() &&
3119 cell->face(0)->at_boundary())
3121 p0 = cell->vertex(v2) - cell->vertex(
v0);
3122 p1 = cell->vertex(
v1) - cell->vertex(
v0);
3123 vertex_at_corner =
v0;
3134 if (angle_fraction > limit_angle_fraction)
3136 auto flags_removal = [&](
unsigned int f1,
3139 unsigned int n2) ->
void {
3140 cells_to_remove[cell->active_cell_index()] =
true;
3141 cells_to_remove[cell->neighbor(n1)->active_cell_index()] =
true;
3142 cells_to_remove[cell->neighbor(n2)->active_cell_index()] =
true;
3144 faces_to_remove[cell->face(f1)->index()] =
true;
3145 faces_to_remove[cell->face(f2)->index()] =
true;
3147 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] =
true;
3148 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] =
true;
3151 auto cell_creation = [&](
const unsigned int vv0,
3152 const unsigned int vv1,
3153 const unsigned int f0,
3154 const unsigned int f1,
3156 const unsigned int n0,
3157 const unsigned int v0n0,
3158 const unsigned int v1n0,
3160 const unsigned int n1,
3161 const unsigned int v0n1,
3162 const unsigned int v1n1) {
3168 c1.
vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
3169 c1.
vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
3175 c2.
vertices[
v1] = cell->neighbor(n1)->vertex_index(v0n1);
3176 c2.
vertices[v2] = cell->vertex_index(vv1);
3177 c2.
vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
3182 l1.
vertices[0] = cell->vertex_index(vv0);
3183 l1.
vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
3189 l2.
vertices[0] = cell->vertex_index(vv0);
3190 l2.
vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
3196 cells_to_add.push_back(c1);
3197 cells_to_add.push_back(c2);
3202 switch (vertex_at_corner)
3205 flags_removal(0, 2, 3, 1);
3206 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
3209 flags_removal(1, 2, 3, 0);
3210 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
3213 flags_removal(3, 0, 1, 2);
3214 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
3217 flags_removal(3, 1, 0, 2);
3218 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
3231 if (cells_to_add.empty())
3233 while (refinement_cycles > 0)
3236 cell->set_coarsen_flag();
3238 refinement_cycles--;
3246 if (cells_to_remove[cell->active_cell_index()] ==
false)
3249 for (
const unsigned int v : cell->vertex_indices())
3250 c.
vertices[v] = cell->vertex_index(v);
3253 cells_to_add.push_back(c);
3261 for (; face != endf; ++face)
3262 if ((face->at_boundary() ||
3264 faces_to_remove[face->index()] ==
false)
3266 for (
unsigned int l = 0;
l < face->
n_lines(); ++
l)
3271 for (
const unsigned int v : face->vertex_indices())
3272 line.
vertices[v] = face->vertex_index(v);
3278 for (
const unsigned int v : face->line(
l)->vertex_indices())
3279 line.
vertices[v] = face->line(
l)->vertex_index(v);
3288 for (
const unsigned int v : face->vertex_indices())
3289 quad.
vertices[v] = face->vertex_index(v);
3297 subcelldata_to_add);
3302 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
3321 template <
int dim,
int spacedim>
3324 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3325 std::vector<std::vector<Point<dim>>>,
3326 std::vector<std::vector<unsigned int>>>
3338 auto &cells = std::get<0>(cqmp);
3339 auto &qpoints = std::get<1>(cqmp);
3340 auto &maps = std::get<2>(cqmp);
3342 return std::make_tuple(std::move(cells),
3349 template <
int dim,
int spacedim>
3352 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3353 std::vector<std::vector<Point<dim>>>,
3354 std::vector<std::vector<unsigned int>>,
3355 std::vector<unsigned int>>
3365 Assert((dim == spacedim),
3366 ExcMessage(
"Only implemented for dim==spacedim."));
3375 const unsigned int np = points.size();
3377 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
3379 std::vector<std::vector<Point<dim>>> qpoints_out;
3380 std::vector<std::vector<unsigned int>> maps_out;
3381 std::vector<unsigned int> missing_points_out;
3385 return std::make_tuple(std::move(cells_out),
3386 std::move(qpoints_out),
3387 std::move(maps_out),
3388 std::move(missing_points_out));
3396 std::vector<std::pair<Point<spacedim>,
unsigned int>> points_and_ids(np);
3397 for (
unsigned int i = 0; i < np; ++i)
3398 points_and_ids[i] = std::make_pair(points[i], i);
3399 const auto p_tree =
pack_rtree(points_and_ids);
3402 std::vector<bool> found_points(points.size(),
false);
3405 const auto already_found = [&found_points](
const auto &id) {
3407 return found_points[
id.second];
3413 const auto store_cell_point_and_id =
3417 const unsigned int &id) {
3418 const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
3419 if (it != cells_out.rend())
3421 const auto cell_id =
3422 (cells_out.size() - 1 - (it - cells_out.rbegin()));
3423 qpoints_out[cell_id].emplace_back(ref_point);
3424 maps_out[cell_id].emplace_back(
id);
3428 cells_out.emplace_back(cell);
3429 qpoints_out.emplace_back(std::vector<
Point<dim>>({ref_point}));
3430 maps_out.emplace_back(std::vector<unsigned int>({
id}));
3435 const auto check_all_points_within_box = [&](
const auto &leaf) {
3436 const double relative_tolerance = 1
e-12;
3439 const auto &cell_hint = leaf.second;
3441 for (
const auto &point_and_id :
3442 p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
3443 bgi::intersects(box)))
3445 const auto id = point_and_id.second;
3446 const auto cell_and_ref =
3450 const auto &cell = cell_and_ref.first;
3451 const auto &ref_point = cell_and_ref.second;
3454 store_cell_point_and_id(cell, ref_point,
id);
3456 missing_points_out.emplace_back(
id);
3459 found_points[id] =
true;
3465 check_all_points_within_box(
3466 std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
3469 for (
unsigned int i = 0; i < np; ++i)
3470 if (found_points[i] ==
false)
3473 const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
3475 if (leaf != b_tree.qend())
3476 check_all_points_within_box(*leaf);
3484 for (
unsigned int i = 0; i < np; ++i)
3485 if (found_points[i] ==
false)
3486 missing_points_out.emplace_back(i);
3493 unsigned int c = cells_out.size();
3494 unsigned int qps = 0;
3500 for (
unsigned int n = 0; n < c; ++n)
3503 qps += qpoints_out[n].size();
3506 Assert(qps + missing_points_out.size() == np,
3510 return std::make_tuple(std::move(cells_out),
3511 std::move(qpoints_out),
3512 std::move(maps_out),
3513 std::move(missing_points_out));
3518 template <
int dim,
int spacedim>
3521 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3522 std::vector<std::vector<Point<dim>>>,
3523 std::vector<std::vector<unsigned int>>,
3524 std::vector<std::vector<Point<spacedim>>>,
3525 std::vector<std::vector<unsigned int>>>
3533 const double tolerance,
3534 const std::vector<bool> &marked_vertices,
3535 const bool enforce_unique_mapping)
3545 enforce_unique_mapping)
3550 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3551 std::vector<std::vector<Point<dim>>>,
3552 std::vector<std::vector<unsigned int>>,
3553 std::vector<std::vector<Point<spacedim>>>,
3554 std::vector<std::vector<unsigned int>>>
3557 std::pair<int, int> dummy{-1, -1};
3559 for (
unsigned int i = 0; i < all.size(); ++i)
3561 if (dummy != std::get<0>(all[i]))
3563 std::get<0>(result).push_back(
3566 std::get<0>(all[i]).
first,
3567 std::get<0>(all[i]).
second});
3569 const unsigned int new_size = std::get<0>(result).size();
3571 std::get<1>(result).resize(new_size);
3572 std::get<2>(result).resize(new_size);
3573 std::get<3>(result).resize(new_size);
3574 std::get<4>(result).resize(new_size);
3576 dummy = std::get<0>(all[i]);
3579 std::get<1>(result).back().push_back(
3580 std::get<3>(all[i]));
3581 std::get<2>(result).back().push_back(std::get<2>(all[i]));
3582 std::get<3>(result).back().push_back(std::get<4>(all[i]));
3583 std::get<4>(result).back().push_back(std::get<1>(all[i]));
3599 template <
int spacedim,
typename T>
3600 std::tuple<std::vector<unsigned int>,
3601 std::vector<unsigned int>,
3602 std::vector<unsigned int>>
3604 const MPI_Comm
comm,
3606 const std::vector<T> &entities,
3607 const double tolerance)
3609 std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
3610 ranks_and_indices.reserve(entities.size());
3612 #if defined(DEAL_II_WITH_ARBORX)
3613 static constexpr
bool use_arborx =
true;
3615 static constexpr
bool use_arborx =
false;
3619 const auto process_bboxes = [&]() ->
void {
3620 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes_temp;
3621 auto *global_bboxes_to_be_used = &global_bboxes;
3623 if (global_bboxes.size() == 1 && use_arborx ==
false)
3625 global_bboxes_temp =
3627 global_bboxes_to_be_used = &global_bboxes_temp;
3631 const auto is_valid = [](
const auto &bb) {
3632 for (
unsigned int i = 0; i < spacedim; ++i)
3633 if (bb.get_boundary_points().first[i] >
3634 bb.get_boundary_points().second[i])
3641 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>>
3644 for (
unsigned rank = 0; rank < global_bboxes_to_be_used->size(); ++rank)
3645 for (
const auto &box : (*global_bboxes_to_be_used)[rank])
3647 boxes_and_ranks.emplace_back(box, rank);
3650 const auto tree =
pack_rtree(boxes_and_ranks);
3653 for (
unsigned int i = 0; i < entities.size(); ++i)
3660 std::set<unsigned int> my_ranks;
3662 for (
const auto &box_and_rank :
3663 tree | boost::geometry::index::adaptors::queried(
3664 boost::geometry::index::intersects(bb)))
3665 my_ranks.insert(box_and_rank.second);
3667 for (
const auto rank : my_ranks)
3668 ranks_and_indices.emplace_back(rank, i);
3672 if constexpr (use_arborx)
3674 if (global_bboxes.size() == 1)
3677 comm, global_bboxes[0]);
3678 std::vector<BoundingBox<spacedim>> query_bounding_boxes;
3679 for (
const auto &entity : entities)
3680 query_bounding_boxes.emplace_back(
3684 query_bounding_boxes);
3685 const auto &[indices_ranks, offsets] =
3686 distributed_tree.
query(bb_intersect);
3688 for (
unsigned long int i = 0; i < offsets.size() - 1; ++i)
3690 std::set<unsigned int> my_ranks;
3691 for (
int j = offsets[i]; j < offsets[i + 1]; ++j)
3692 my_ranks.insert(indices_ranks[j].second);
3694 for (
const auto rank : my_ranks)
3695 ranks_and_indices.emplace_back(rank, i);
3712 std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
3714 std::vector<unsigned int> ranks;
3715 std::vector<unsigned int> ptr;
3716 std::vector<unsigned int> indices;
3720 for (
const std::pair<unsigned int, unsigned int> &i : ranks_and_indices)
3722 if (current_rank != i.first)
3724 current_rank = i.first;
3725 ranks.push_back(current_rank);
3726 ptr.push_back(indices.size());
3729 indices.push_back(i.second);
3731 ptr.push_back(indices.size());
3733 return {std::move(ranks), std::move(ptr), std::move(indices)};
3738 template <
int dim,
int spacedim>
3740 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
3746 const std::vector<bool> &marked_vertices,
3747 const double tolerance,
3748 const bool enforce_unique_mapping)
3751 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
3753 locally_owned_active_cells_around_point;
3770 cell_hint = first_cell.first;
3773 const auto active_cells_around_point =
3782 if (enforce_unique_mapping)
3789 for (
const auto &cell : active_cells_around_point)
3790 lowes_rank =
std::min(lowes_rank, cell.first->subdomain_id());
3792 if (lowes_rank != my_rank)
3796 locally_owned_active_cells_around_point.reserve(
3797 active_cells_around_point.size());
3799 for (
const auto &cell : active_cells_around_point)
3800 if (cell.first->is_locally_owned())
3801 locally_owned_active_cells_around_point.push_back(cell);
3804 std::sort(locally_owned_active_cells_around_point.begin(),
3805 locally_owned_active_cells_around_point.end(),
3806 [](
const auto &a,
const auto &
b) { return a.first < b.first; });
3808 if (enforce_unique_mapping &&
3809 locally_owned_active_cells_around_point.size() > 1)
3811 return {locally_owned_active_cells_around_point.front()};
3813 return locally_owned_active_cells_around_point;
3816 template <
int dim,
int spacedim>
3822 template <
int dim,
int spacedim>
3830 Assert(recv_components.empty() ||
3831 std::get<1>(*std::max_element(recv_components.begin(),
3832 recv_components.end(),
3833 [](
const auto &a,
const auto &
b) {
3834 return std::get<1>(a) <
3836 })) < n_searched_points,
3848 std::sort(send_components.begin(),
3849 send_components.end(),
3850 [&](
const auto &a,
const auto &
b) {
3851 if (std::get<1>(a) != std::get<1>(b))
3852 return std::get<1>(a) < std::get<1>(b);
3854 if (std::get<2>(a) != std::get<2>(b))
3855 return std::get<2>(a) < std::get<2>(b);
3857 return std::get<0>(a) < std::get<0>(b);
3862 i < send_components.size();
3865 std::get<5>(send_components[i]) = i;
3867 if (dummy != std::get<1>(send_components[i]))
3869 dummy = std::get<1>(send_components[i]);
3870 send_ranks.push_back(dummy);
3871 send_ptrs.push_back(i);
3874 send_ptrs.push_back(send_components.size());
3878 std::sort(send_components.begin(),
3879 send_components.end(),
3880 [&](
const auto &a,
const auto &
b) {
3881 if (std::get<0>(a) != std::get<0>(b))
3882 return std::get<0>(a) < std::get<0>(b);
3884 if (std::get<1>(a) != std::get<1>(b))
3885 return std::get<1>(a) < std::get<1>(b);
3887 if (std::get<2>(a) != std::get<2>(b))
3888 return std::get<2>(a) < std::get<2>(b);
3890 return std::get<5>(a) < std::get<5>(b);
3894 if (recv_components.size() > 0)
3897 std::sort(recv_components.begin(),
3898 recv_components.end(),
3899 [&](
const auto &a,
const auto &
b) {
3900 if (std::get<0>(a) != std::get<0>(b))
3901 return std::get<0>(a) < std::get<0>(b);
3903 return std::get<1>(a) < std::get<1>(b);
3908 i < recv_components.size();
3911 std::get<2>(recv_components[i]) = i;
3913 if (dummy != std::get<0>(recv_components[i]))
3915 dummy = std::get<0>(recv_components[i]);
3916 recv_ranks.push_back(dummy);
3917 recv_ptrs.push_back(i);
3920 recv_ptrs.push_back(recv_components.size());
3924 std::sort(recv_components.begin(),
3925 recv_components.end(),
3926 [&](
const auto &a,
const auto &
b) {
3927 if (std::get<1>(a) != std::get<1>(b))
3928 return std::get<1>(a) < std::get<1>(b);
3930 if (std::get<0>(a) != std::get<0>(b))
3931 return std::get<0>(a) < std::get<0>(b);
3933 return std::get<2>(a) < std::get<2>(b);
3940 template <
int dim,
int spacedim>
3946 const std::vector<bool> &marked_vertices,
3947 const double tolerance,
3948 const bool perform_handshake,
3949 const bool enforce_unique_mapping)
3960 comm, global_bboxes, points, tolerance);
3962 const auto &potential_owners_ranks = std::get<0>(potential_owners);
3963 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
3964 const auto &potential_owners_indices = std::get<2>(potential_owners);
3968 const auto translate = [&](
const unsigned int other_rank) {
3969 const auto ptr = std::find(potential_owners_ranks.begin(),
3970 potential_owners_ranks.end(),
3975 const auto other_rank_index =
3976 std::distance(potential_owners_ranks.begin(), ptr);
3978 return other_rank_index;
3982 (marked_vertices.empty()) ||
3985 "The marked_vertices vector has to be either empty or its size has "
3986 "to equal the number of vertices of the triangulation."));
3988 using RequestType = std::vector<std::pair<unsigned int, Point<spacedim>>>;
3989 using AnswerType = std::vector<unsigned int>;
3995 const bool has_relevant_vertices =
3996 (marked_vertices.empty()) ||
3997 (std::find(marked_vertices.begin(), marked_vertices.end(),
true) !=
3998 marked_vertices.end());
4000 const auto create_request = [&](
const unsigned int other_rank) {
4001 const auto other_rank_index = translate(other_rank);
4003 RequestType request;
4004 request.reserve(potential_owners_ptrs[other_rank_index + 1] -
4005 potential_owners_ptrs[other_rank_index]);
4007 for (
unsigned int i = potential_owners_ptrs[other_rank_index];
4008 i < potential_owners_ptrs[other_rank_index + 1];
4010 request.emplace_back(potential_owners_indices[i],
4011 points[potential_owners_indices[i]]);
4016 const auto answer_request =
4017 [&](
const unsigned int &other_rank,
4018 const RequestType &request) -> AnswerType {
4019 AnswerType answer(request.size(), 0);
4021 if (has_relevant_vertices)
4025 for (
unsigned int i = 0; i < request.size(); ++i)
4027 const auto &index_and_point = request[i];
4029 const auto cells_and_reference_positions =
4032 index_and_point.second,
4036 enforce_unique_mapping);
4041 for (
const auto &cell_and_reference_position :
4042 cells_and_reference_positions)
4044 const auto cell = cell_and_reference_position.first;
4045 auto reference_position =
4046 cell_and_reference_position.second;
4048 reference_position =
4049 cell->reference_cell().closest_point(reference_position);
4051 send_components.emplace_back(
4052 std::pair<int, int>(cell->level(), cell->index()),
4054 index_and_point.first,
4056 index_and_point.second,
4060 answer[i] = cells_and_reference_positions.size();
4064 if (perform_handshake)
4070 const auto process_answer = [&](
const unsigned int other_rank,
4071 const AnswerType &answer) {
4072 if (perform_handshake)
4074 const auto other_rank_index = translate(other_rank);
4076 for (
unsigned int i = 0; i < answer.size(); ++i)
4077 for (
unsigned int j = 0; j < answer[i]; ++j)
4078 recv_components.emplace_back(
4080 potential_owners_indices
4081 [i + potential_owners_ptrs[other_rank_index]],
4086 Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
4087 potential_owners_ranks,
4100 template <
int structdim,
int spacedim>
4102 DistributedComputePointLocationsInternal<dim, spacedim>
4105 const unsigned int n_points_1D,
4108 const bool consistent_numbering_of_sender_and_receiver)
const
4110 using CellIterator =
4126 for (
const auto &recv_component : recv_components)
4133 std::get<2>(recv_component));
4135 for (
unsigned int i = 0; i < quad.
size(); ++i)
4139 result.recv_components.emplace_back(
4140 std::get<0>(recv_component),
4141 result.recv_components.size(),
4149 result.n_searched_points = result.recv_components.size();
4154 std::map<unsigned int, std::vector<unsigned int>> indices_of_rank;
4155 std::map<unsigned int, unsigned int> counter;
4156 std::set<unsigned int> send_ranks;
4157 if (consistent_numbering_of_sender_and_receiver)
4159 for (
const auto &sc : send_components)
4160 send_ranks.insert(std::get<1>(sc));
4162 for (
const auto rank : send_ranks)
4166 indices_of_rank = communicate_indices(result.recv_components,
4170 for (
const auto &send_component : send_components)
4172 const CellIterator cell(&
tria,
4173 std::get<0>(send_component).
first,
4174 std::get<0>(send_component).
second);
4178 std::get<3>(send_component));
4180 const auto rank = std::get<1>(send_component);
4182 for (
unsigned int q = 0; q < quad.
size(); ++q)
4186 result.send_components.emplace_back(std::make_tuple(
4187 std::get<0>(send_component),
4189 indices_of_rank.empty() ?
4190 result.send_components.size() :
4191 indices_of_rank.at(rank)[counter.at(rank)],
4196 if (!indices_of_rank.empty())
4201 result.finalize_setup();
4208 template <
int structdim,
int spacedim>
4209 std::map<unsigned int, std::vector<unsigned int>>
4212 const std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
4213 &point_recv_components,
4214 const MPI_Comm
comm)
const
4216 #ifndef DEAL_II_WITH_MPI
4218 (void)point_recv_components;
4224 const auto mpi_tag =
4229 std::set<unsigned int> send_ranks;
4230 for (
const auto &sc : send_components)
4231 send_ranks.insert(std::get<1>(sc));
4232 std::set<unsigned int> recv_ranks;
4233 for (
const auto &rc : recv_components)
4234 recv_ranks.insert(std::get<0>(rc));
4236 std::vector<MPI_Request> requests;
4237 requests.reserve(send_ranks.size());
4240 std::map<unsigned int, std::vector<unsigned int>> indices_of_rank;
4241 indices_of_rank[my_rank] = std::vector<unsigned int>();
4244 std::map<unsigned int, std::vector<unsigned int>> send_indices_of_rank;
4245 for (
const auto rank : recv_ranks)
4246 if (rank != my_rank)
4247 send_indices_of_rank[rank] = std::vector<unsigned int>();
4250 for (
const auto &point_recv_component : point_recv_components)
4252 const auto rank = std::get<0>(point_recv_component);
4253 const auto idx = std::get<1>(point_recv_component);
4255 if (rank == my_rank)
4256 indices_of_rank[rank].emplace_back(idx);
4258 send_indices_of_rank[rank].emplace_back(idx);
4262 for (
const auto rank : recv_ranks)
4264 if (rank == my_rank)
4269 requests.push_back(MPI_Request());
4271 const int ierr = MPI_Isend(buffer.data(),
4282 for (
const auto rank : send_ranks)
4284 if (rank == my_rank)
4288 int ierr = MPI_Probe(MPI_ANY_SOURCE, mpi_tag,
comm, &status);
4292 ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
4295 std::vector<char> buffer(message_length);
4297 ierr = MPI_Recv(buffer.data(),
4306 indices_of_rank[status.MPI_SOURCE] =
4307 Utilities::unpack<std::vector<unsigned int>>(buffer,
false);
4312 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4315 return indices_of_rank;
4321 template <
int structdim,
int dim,
int spacedim>
4325 const std::vector<std::vector<
Point<spacedim>>> &intersection_requests,
4327 const std::vector<bool> &marked_vertices,
4328 const double tolerance)
4330 using IntersectionRequest = std::vector<Point<spacedim>>;
4332 using IntersectionAnswer =
4335 spacedim>::IntersectionType;
4348 comm, global_bboxes, intersection_requests, tolerance);
4350 const auto &potential_owners_ranks = std::get<0>(potential_owners);
4351 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
4352 const auto &potential_owners_indices = std::get<2>(potential_owners);
4354 const auto translate = [&](
const unsigned int other_rank) {
4355 const auto ptr = std::find(potential_owners_ranks.begin(),
4356 potential_owners_ranks.end(),
4361 const auto other_rank_index =
4362 std::distance(potential_owners_ranks.begin(), ptr);
4364 return other_rank_index;
4368 (marked_vertices.empty()) ||
4371 "The marked_vertices vector has to be either empty or its size has "
4372 "to equal the number of vertices of the triangulation."));
4378 const bool has_relevant_vertices =
4379 (marked_vertices.empty()) ||
4380 (std::find(marked_vertices.begin(), marked_vertices.end(),
true) !=
4381 marked_vertices.end());
4387 std::vector<std::pair<unsigned int, IntersectionRequest>>;
4391 std::vector<std::pair<unsigned int, IntersectionAnswer>>;
4393 const auto create_request = [&](
const unsigned int other_rank) {
4394 const auto other_rank_index = translate(other_rank);
4396 RequestType request;
4397 request.reserve(potential_owners_ptrs[other_rank_index + 1] -
4398 potential_owners_ptrs[other_rank_index]);
4400 for (
unsigned int i = potential_owners_ptrs[other_rank_index];
4401 i < potential_owners_ptrs[other_rank_index + 1];
4403 request.emplace_back(
4404 potential_owners_indices[i],
4405 intersection_requests[potential_owners_indices[i]]);
4413 const auto construct_locally_owned_cell_bounding_boxes_rtree =
4414 [&cache](
const std::vector<bool> &marked_verts) {
4415 const auto cell_marked = [&marked_verts](
const auto &cell) {
4416 for (
const unsigned int v : cell->vertex_indices())
4417 if (marked_verts[cell->vertex_index(v)])
4422 const auto &boxes_and_cells =
4425 if (marked_verts.empty())
4426 return boxes_and_cells;
4428 std::vector<std::pair<
4431 potential_boxes_and_cells;
4433 for (
const auto &box_and_cell : boxes_and_cells)
4434 if (cell_marked(box_and_cell.second))
4435 potential_boxes_and_cells.emplace_back(box_and_cell);
4437 return pack_rtree(potential_boxes_and_cells);
4442 std::pair<BoundingBox<spacedim>,
4446 const auto answer_request =
4447 [&](
const unsigned int &other_rank,
4448 const RequestType &request) -> AnswerType {
4451 if (has_relevant_vertices)
4453 if (marked_cell_tree.empty())
4456 construct_locally_owned_cell_bounding_boxes_rtree(
4461 for (
unsigned int i = 0; i < request.size(); ++i)
4467 for (
const auto &box_cell :
4469 boost::geometry::index::adaptors::queried(
4470 boost::geometry::index::intersects(bb)))
4472 #ifdef DEAL_II_WITH_CGAL
4473 const auto &cell = box_cell.second;
4474 const auto &request_index = request[i].first;
4475 auto requested_intersection = request[i].second;
4476 CGALWrappers::resort_dealii_vertices_to_cgal_order(
4477 structdim, requested_intersection);
4479 const auto &try_intersection =
4480 CGALWrappers::get_vertices_in_cgal_order(
4483 const auto &found_intersections = CGALWrappers::
4484 compute_intersection_of_cells<dim, structdim, spacedim>(
4485 try_intersection, requested_intersection, tolerance);
4487 if (found_intersections.size() > 0)
4489 for (
const auto &found_intersection :
4490 found_intersections)
4492 answer.emplace_back(request_index,
4493 found_intersection);
4495 send_components.emplace_back(
4496 std::make_pair(cell->level(), cell->index()),
4499 found_intersection);
4514 const auto process_answer = [&](
const unsigned int other_rank,
4515 const AnswerType &answer) {
4516 for (
unsigned int i = 0; i < answer.size(); ++i)
4517 recv_components.emplace_back(other_rank,
4522 Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
4523 potential_owners_ranks,
4531 std::stable_sort(recv_components.begin(),
4532 recv_components.end(),
4533 [&](
const auto &a,
const auto &
b) {
4535 if (std::get<1>(a) != std::get<1>(b))
4536 return std::get<1>(a) < std::get<1>(b);
4539 return std::get<0>(a) < std::get<0>(b);
4544 std::stable_sort(send_components.begin(),
4545 send_components.end(),
4546 [&](
const auto &a,
const auto &
b) {
4548 if (std::get<1>(a) != std::get<1>(b))
4549 return std::get<1>(a) < std::get<1>(b);
4552 return std::get<2>(a) < std::get<2>(b);
4556 recv_ptrs.assign(intersection_requests.size() + 1, 0);
4557 for (
const auto &rc : recv_components)
4558 ++recv_ptrs[std::get<1>(rc) + 1];
4559 for (
unsigned int i = 0; i < intersection_requests.size(); ++i)
4560 recv_ptrs[i + 1] += recv_ptrs[i];
4569 template <
int spacedim>
4574 auto id_and_v = std::min_element(
4579 return p1.second.distance(p) < p2.second.distance(p);
4581 return id_and_v->first;
4585 template <
int dim,
int spacedim>
4586 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4593 const std::vector<bool> &marked_vertices,
4594 const double tolerance)
4599 const auto &vertex_to_cell_centers =
4607 vertex_to_cell_centers,
4610 used_vertices_rtree,
4614 template <
int spacedim>
4615 std::vector<std::vector<BoundingBox<spacedim>>>
4618 const MPI_Comm mpi_communicator)
4620 #ifndef DEAL_II_WITH_MPI
4622 (void)mpi_communicator;
4625 "GridTools::exchange_local_bounding_boxes() requires MPI."));
4629 unsigned int n_bboxes = local_bboxes.size();
4631 int n_local_data = 2 * spacedim * n_bboxes;
4634 std::vector<double> loc_data_array(n_local_data);
4635 for (
unsigned int i = 0; i < n_bboxes; ++i)
4636 for (
unsigned int d = 0;
d < spacedim; ++
d)
4639 loc_data_array[2 * i * spacedim +
d] =
4640 local_bboxes[i].get_boundary_points().first[
d];
4641 loc_data_array[2 * i * spacedim + spacedim +
d] =
4642 local_bboxes[i].get_boundary_points().second[
d];
4649 std::vector<int> size_all_data(n_procs);
4652 int ierr = MPI_Allgather(&n_local_data,
4655 size_all_data.data(),
4663 std::vector<int> rdispls(n_procs);
4665 for (
unsigned int i = 1; i < n_procs; ++i)
4666 rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
4670 std::vector<double> data_array(rdispls.back() + size_all_data.back());
4672 ierr = MPI_Allgatherv(loc_data_array.data(),
4676 size_all_data.data(),
4683 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
4684 unsigned int begin_idx = 0;
4685 for (
unsigned int i = 0; i < n_procs; ++i)
4688 unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
4689 global_bboxes[i].resize(n_bbox_i);
4690 for (
unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
4693 for (
unsigned int d = 0;
d < spacedim; ++
d)
4695 p1[
d] = data_array[begin_idx + 2 * bbox * spacedim +
d];
4697 data_array[begin_idx + 2 * bbox * spacedim + spacedim +
d];
4700 global_bboxes[i][bbox] = loc_bbox;
4703 begin_idx += size_all_data[i];
4705 return global_bboxes;
4711 template <
int spacedim>
4715 const MPI_Comm mpi_communicator)
4717 #ifndef DEAL_II_WITH_MPI
4718 (void)mpi_communicator;
4720 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> boxes_index(
4721 local_description.size());
4723 for (
unsigned int i = 0; i < local_description.size(); ++i)
4724 boxes_index[i] = std::make_pair(local_description[i], 0u);
4728 const std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes =
4732 const unsigned int n_procs =
4736 std::vector<unsigned int> bboxes_position(n_procs);
4738 unsigned int tot_bboxes = 0;
4739 for (
const auto &process_bboxes : global_bboxes)
4740 tot_bboxes += process_bboxes.size();
4743 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>>
4745 flat_global_bboxes.reserve(tot_bboxes);
4746 unsigned int process_index = 0;
4747 for (
const auto &process_bboxes : global_bboxes)
4750 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>>
4751 boxes_and_indices(process_bboxes.size());
4754 for (
unsigned int i = 0; i < process_bboxes.size(); ++i)
4755 boxes_and_indices[i] =
4756 std::make_pair(process_bboxes[i], process_index);
4758 flat_global_bboxes.insert(flat_global_bboxes.end(),
4759 boxes_and_indices.begin(),
4760 boxes_and_indices.end());
4768 flat_global_bboxes.begin(), flat_global_bboxes.end());
4774 template <
int dim,
int spacedim>
4784 static const int lookup_table_2d[2][2] =
4791 static const int lookup_table_3d[2][2][2][4] =
4813 if (pair.first.first->level() != pair.second.first.first->level())
4816 const auto face_a = pair.first.first->face(pair.first.second);
4818 pair.second.first.first->face(pair.second.first.second);
4819 const auto mask = pair.second.second;
4824 for (
unsigned int i = 0; i < face_a->n_vertices(); ++i)
4826 const bool face_orientation = mask[0];
4827 const bool face_flip = mask[1];
4828 const bool face_rotation = mask[2];
4838 j = lookup_table_2d[face_flip][i];
4841 j = lookup_table_3d[face_orientation][face_flip]
4849 const auto vertex_a = face_a->vertex_index(i);
4850 const auto vertex_b = face_b->vertex_index(j);
4851 unsigned int temp =
std::min(vertex_a, vertex_b);
4855 temp =
std::min(temp, it_a->second);
4859 temp =
std::min(temp, it_b->second);
4862 it_a->second = temp;
4867 it_b->second = temp;
4877 if (p.first == p.second)
4879 unsigned int temp = p.second;
4897 template <
int dim,
int spacedim>
4898 std::map<unsigned int, std::set<::types::subdomain_id>>
4917 for (
const auto &cell :
4919 for (
const unsigned int v : cell->vertex_indices())
4920 vertex_of_own_cell[cell->vertex_index(v)] =
true;
4924 std::map<unsigned int, std::set<types::subdomain_id>> result;
4928 if (cell->is_ghost())
4933 for (
const unsigned int v : cell->vertex_indices())
4936 if (vertex_of_own_cell[cell->vertex_index(v)])
4937 result[cell->vertex_index(v)].insert(owner);
4940 auto coinciding_vertex_group =
4942 if (coinciding_vertex_group !=
4944 for (
auto coinciding_vertex :
4946 if (vertex_of_own_cell[coinciding_vertex])
4947 result[coinciding_vertex].insert(owner);
4959 unsigned int n_vertices,
4960 unsigned int n_sub_vertices,
4961 unsigned int n_configurations,
4962 unsigned int n_lines,
4963 unsigned int n_cols,
4964 typename value_type>
4967 const std::array<unsigned int, n_configurations> &cut_line_table,
4970 const std::vector<value_type> &ls_values,
4972 const std::vector<unsigned int> &mask,
4973 const double iso_level,
4974 const double tolerance,
4977 const bool write_back_cell_data)
4981 constexpr
unsigned int X =
static_cast<unsigned int>(-1);
4984 unsigned int configuration = 0;
4985 for (
unsigned int v = 0; v < n_vertices; ++v)
4986 if (ls_values[mask[v]] < iso_level)
4987 configuration |= (1 << v);
4990 if (cut_line_table[configuration] == 0)
4995 const auto interpolate = [&](
const unsigned int i,
const unsigned int j) {
4996 if (std::abs(iso_level - ls_values[mask[i]]) < tolerance)
4997 return points[mask[i]];
4998 if (std::abs(iso_level - ls_values[mask[j]]) < tolerance)
4999 return points[mask[j]];
5000 if (std::abs(ls_values[mask[i]] - ls_values[mask[j]]) < tolerance)
5001 return points[mask[i]];
5003 const double mu = (iso_level - ls_values[mask[i]]) /
5004 (ls_values[mask[j]] - ls_values[mask[i]]);
5007 mu * (points[mask[j]] - points[mask[i]]));
5011 std::array<Point<dim>, n_lines> vertex_list_all;
5012 for (
unsigned int l = 0;
l < n_lines; ++
l)
5013 if (cut_line_table[configuration] & (1 <<
l))
5014 vertex_list_all[
l] =
5015 interpolate(line_to_vertex_table[
l][0], line_to_vertex_table[
l][1]);
5018 unsigned int local_vertex_count = 0;
5019 std::array<Point<dim>, n_lines> vertex_list_reduced;
5020 std::array<unsigned int, n_lines> local_remap;
5021 std::fill(local_remap.begin(), local_remap.end(), X);
5022 for (
int i = 0; new_line_table[configuration][i] != X; ++i)
5023 if (local_remap[new_line_table[configuration][i]] == X)
5025 vertex_list_reduced[local_vertex_count] =
5026 vertex_list_all[new_line_table[configuration][i]];
5027 local_remap[new_line_table[configuration][i]] = local_vertex_count;
5028 local_vertex_count++;
5032 const unsigned int n_vertices_old =
vertices.size();
5033 for (
unsigned int i = 0; i < local_vertex_count; ++i)
5034 vertices.push_back(vertex_list_reduced[i]);
5037 if (write_back_cell_data && dim > 1)
5039 for (
unsigned int i = 0; new_line_table[configuration][i] != X;
5040 i += n_sub_vertices)
5042 cells.resize(cells.size() + 1);
5043 cells.back().vertices.resize(n_sub_vertices);
5045 for (
unsigned int v = 0; v < n_sub_vertices; ++v)
5046 cells.back().vertices[v] =
5047 local_remap[new_line_table[configuration][i + v]] +
5056 template <
int dim,
typename VectorType>
5060 const unsigned int n_subdivisions,
5061 const double tolerance)
5062 : n_subdivisions(n_subdivisions)
5063 , tolerance(tolerance)
5064 , fe_values(mapping,
5066 create_quadrature_rule(n_subdivisions),
5072 template <
int dim,
typename VectorType>
5075 const unsigned int n_subdivisions)
5081 for (
unsigned int i = 0; i <= n_subdivisions; ++i)
5086 for (
unsigned int j = 0; j <= n_subdivisions; ++j)
5087 for (
unsigned int i = 0; i <= n_subdivisions; ++i)
5089 1.0 / n_subdivisions * j);
5093 for (
unsigned int k = 0; k <= n_subdivisions; ++k)
5094 for (
unsigned int j = 0; j <= n_subdivisions; ++j)
5095 for (
unsigned int i = 0; i <= n_subdivisions; ++i)
5097 1.0 / n_subdivisions * j,
5098 1.0 / n_subdivisions * k);
5107 template <
int dim,
typename VectorType>
5111 const VectorType &ls_vector,
5112 const double iso_level,
5119 "Not implemented for dim==1. Use the alternative process()-function "
5120 "not returning a vector of CellData objects."));
5124 process_cell(cell, ls_vector, iso_level,
vertices, cells);
5127 template <
int dim,
typename VectorType>
5131 const VectorType &ls_vector,
5132 const double iso_level,
5137 process_cell(cell, ls_vector, iso_level,
vertices);
5143 template <
int dim,
typename VectorType>
5147 const VectorType &ls_vector,
5148 const double iso_level,
5155 "Not implemented for dim==1. Use the alternative process_cell()-function "
5156 "not returning a vector of CellData objects."));
5158 std::vector<value_type> ls_values;
5160 fe_values.reinit(cell);
5161 ls_values.resize(fe_values.n_quadrature_points);
5162 fe_values.get_function_values(ls_vector, ls_values);
5164 ls_values, fe_values.get_quadrature_points(), iso_level,
vertices, cells);
5167 template <
int dim,
typename VectorType>
5171 const VectorType &ls_vector,
5172 const double iso_level,
5176 std::vector<
CellData<dim == 1 ? 1 : dim - 1>> dummy_cells;
5178 std::vector<value_type> ls_values;
5180 fe_values.reinit(cell);
5181 ls_values.resize(fe_values.n_quadrature_points);
5182 fe_values.get_function_values(ls_vector, ls_values);
5184 process_cell(ls_values,
5185 fe_values.get_quadrature_points(),
5193 template <
int dim,
typename VectorType>
5196 std::vector<value_type> &ls_values,
5198 const double iso_level,
5201 const bool write_back_cell_data)
const
5203 const unsigned p = n_subdivisions + 1;
5207 for (
unsigned int i = 0; i < n_subdivisions; ++i)
5209 std::vector<unsigned int> mask{i + 0, i + 1};
5212 if (std::abs(iso_level - ls_values[mask[0]]) < tolerance)
5213 vertices.emplace_back(points[mask[0]]);
5214 else if (std::abs(iso_level - ls_values[mask[1]]) < tolerance)
5216 if (i + 1 == n_subdivisions)
5217 vertices.emplace_back(points[mask[1]]);
5220 else if (((ls_values[mask[0]] > iso_level) &&
5221 (ls_values[mask[1]] < iso_level)) ||
5222 ((ls_values[mask[0]] < iso_level) &&
5223 (ls_values[mask[1]] > iso_level)))
5226 const double mu = (iso_level - ls_values[mask[0]]) /
5227 (ls_values[mask[1]] - ls_values[mask[0]]);
5230 vertices.emplace_back(points[mask[0]] +
5231 mu * (points[mask[1]] - points[mask[0]]));
5237 for (
unsigned int j = 0; j < n_subdivisions; ++j)
5238 for (
unsigned int i = 0; i < n_subdivisions; ++i)
5240 std::vector<unsigned int> mask{p * (j + 0) + (i + 0),
5241 p * (j + 0) + (i + 1),
5242 p * (j + 1) + (i + 1),
5243 p * (j + 1) + (i + 0)};
5251 write_back_cell_data);
5256 for (
unsigned int k = 0; k < n_subdivisions; ++k)
5257 for (
unsigned int j = 0; j < n_subdivisions; ++j)
5258 for (
unsigned int i = 0; i < n_subdivisions; ++i)
5260 std::vector<unsigned int> mask{
5261 p * p * (k + 0) + p * (j + 0) + (i + 0),
5262 p * p * (k + 0) + p * (j + 0) + (i + 1),
5263 p * p * (k + 0) + p * (j + 1) + (i + 1),
5264 p * p * (k + 0) + p * (j + 1) + (i + 0),
5265 p * p * (k + 1) + p * (j + 0) + (i + 0),
5266 p * p * (k + 1) + p * (j + 0) + (i + 1),
5267 p * p * (k + 1) + p * (j + 1) + (i + 1),
5268 p * p * (k + 1) + p * (j + 1) + (i + 0)};
5276 write_back_cell_data);
5283 template <
int dim,
typename VectorType>
5286 const std::vector<value_type> &ls_values,
5287 const std::vector<
Point<2>> &points,
5288 const std::vector<unsigned int> &mask,
5289 const double iso_level,
5292 const bool write_back_cell_data)
const
5295 constexpr
unsigned int n_vertices = 4;
5296 constexpr
unsigned int n_sub_vertices = 2;
5297 constexpr
unsigned int n_lines = 4;
5298 constexpr
unsigned int n_configurations =
Utilities::pow(2, n_vertices);
5299 constexpr
unsigned int X =
static_cast<unsigned int>(-1);
5303 constexpr std::array<unsigned int, n_configurations> cut_line_table = {
5339 {{X, X, X, X, X}}}};
5343 {{{0, 3}}, {{1, 2}}, {{0, 1}}, {{3, 2}}}};
5353 line_to_vertex_table,
5361 write_back_cell_data);
5366 template <
int dim,
typename VectorType>
5369 const std::vector<value_type> &ls_values,
5370 const std::vector<
Point<3>> &points,
5371 const std::vector<unsigned int> &mask,
5372 const double iso_level,
5375 const bool write_back_cell_data)
const
5378 constexpr
unsigned int n_vertices = 8;
5379 constexpr
unsigned int n_sub_vertices = 3;
5380 constexpr
unsigned int n_lines = 12;
5381 constexpr
unsigned int n_configurations =
Utilities::pow(2, n_vertices);
5382 constexpr
unsigned int X =
static_cast<unsigned int>(-1);
5387 constexpr std::array<unsigned int, n_configurations> cut_line_table = {{
5388 0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905,
5389 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, 0x190, 0x99, 0x393, 0x29a,
5390 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93,
5391 0xf99, 0xe90, 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, 0x53c,
5392 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, 0x3a0, 0x2a9,
5393 0x1a3, 0xaa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6,
5394 0xfaa, 0xea3, 0xda9, 0xca0, 0x460, 0x569, 0x663, 0x76a, 0x66, 0x16f,
5395 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
5396 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5,
5397 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, 0x650, 0x759, 0x453, 0x55a,
5398 0x256, 0x35f, 0x55, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53,
5399 0x859, 0x950, 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc,
5400 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, 0x8c0, 0x9c9,
5401 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6,
5402 0x4ca, 0x5c3, 0x6c9, 0x7c0, 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f,
5403 0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
5404 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5,
5405 0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, 0xb60, 0xa69, 0x963, 0x86a,
5406 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663,
5407 0x569, 0x460, 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
5408 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, 0xd30, 0xc39,
5409 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636,
5410 0x13a, 0x33, 0x339, 0x230, 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f,
5411 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190,
5412 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605,
5413 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0}};
5419 {{{X, X, X, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5420 {{0, 8, 3, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5421 {{0, 1, 9, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5422 {{1, 8, 3, 9, 8, 1, X, X, X, X, X, X, X, X, X, X}},
5423 {{1, 2, 10, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5424 {{0, 8, 3, 1, 2, 10, X, X, X, X, X, X, X, X, X, X}},
5425 {{9, 2, 10, 0, 2, 9, X, X, X, X, X, X, X, X, X, X}},
5426 {{2, 8, 3, 2, 10, 8, 10, 9, 8, X, X, X, X, X, X, X}},
5427 {{3, 11, 2, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5428 {{0, 11, 2, 8, 11, 0, X, X, X, X, X, X, X, X, X, X}},
5429 {{1, 9, 0, 2, 3, 11, X, X, X, X, X, X, X, X, X, X}},
5430 {{1, 11, 2, 1, 9, 11, 9, 8, 11, X, X, X, X, X, X, X}},
5431 {{3, 10, 1, 11, 10, 3, X, X, X, X, X, X, X, X, X, X}},
5432 {{0, 10, 1, 0, 8, 10, 8, 11, 10, X, X, X, X, X, X, X}},
5433 {{3, 9, 0, 3, 11, 9, 11, 10, 9, X, X, X, X, X, X, X}},
5434 {{9, 8, 10, 10, 8, 11, X, X, X, X, X, X, X, X, X, X}},
5435 {{4, 7, 8, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5436 {{4, 3, 0, 7, 3, 4, X, X, X, X, X, X, X, X, X, X}},
5437 {{0, 1, 9, 8, 4, 7, X, X, X, X, X, X, X, X, X, X}},
5438 {{4, 1, 9, 4, 7, 1, 7, 3, 1, X, X, X, X, X, X, X}},
5439 {{1, 2, 10, 8, 4, 7, X, X, X, X, X, X, X, X, X, X}},
5440 {{3, 4, 7, 3, 0, 4, 1, 2, 10, X, X, X, X, X, X, X}},
5441 {{9, 2, 10, 9, 0, 2, 8, 4, 7, X, X, X, X, X, X, X}},
5442 {{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, X, X, X, X}},
5443 {{8, 4, 7, 3, 11, 2, X, X, X, X, X, X, X, X, X, X}},
5444 {{11, 4, 7, 11, 2, 4, 2, 0, 4, X, X, X, X, X, X, X}},
5445 {{9, 0, 1, 8, 4, 7, 2, 3, 11, X, X, X, X, X, X, X}},
5446 {{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, X, X, X, X}},
5447 {{3, 10, 1, 3, 11, 10, 7, 8, 4, X, X, X, X, X, X, X}},
5448 {{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, X, X, X, X}},
5449 {{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, X, X, X, X}},
5450 {{4, 7, 11, 4, 11, 9, 9, 11, 10, X, X, X, X, X, X, X}},
5451 {{9, 5, 4, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5452 {{9, 5, 4, 0, 8, 3, X, X, X, X, X, X, X, X, X, X}},
5453 {{0, 5, 4, 1, 5, 0, X, X, X, X, X, X, X, X, X, X}},
5454 {{8, 5, 4, 8, 3, 5, 3, 1, 5, X, X, X, X, X, X, X}},
5455 {{1, 2, 10, 9, 5, 4, X, X, X, X, X, X, X, X, X, X}},
5456 {{3, 0, 8, 1, 2, 10, 4, 9, 5, X, X, X, X, X, X, X}},
5457 {{5, 2, 10, 5, 4, 2, 4, 0, 2, X, X, X, X, X, X, X}},
5458 {{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, X, X, X, X}},
5459 {{9, 5, 4, 2, 3, 11, X, X, X, X, X, X, X, X, X, X}},
5460 {{0, 11, 2, 0, 8, 11, 4, 9, 5, X, X, X, X, X, X, X}},
5461 {{0, 5, 4, 0, 1, 5, 2, 3, 11, X, X, X, X, X, X, X}},
5462 {{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, X, X, X, X}},
5463 {{10, 3, 11, 10, 1, 3, 9, 5, 4, X, X, X, X, X, X, X}},
5464 {{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, X, X, X, X}},
5465 {{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, X, X, X, X}},
5466 {{5, 4, 8, 5, 8, 10, 10, 8, 11, X, X, X, X, X, X, X}},
5467 {{9, 7, 8, 5, 7, 9, X, X, X, X, X, X, X, X, X, X}},
5468 {{9, 3, 0, 9, 5, 3, 5, 7, 3, X, X, X, X, X, X, X}},
5469 {{0, 7, 8, 0, 1, 7, 1, 5, 7, X, X, X, X, X, X, X}},
5470 {{1, 5, 3, 3, 5, 7, X, X, X, X, X, X, X, X, X, X}},
5471 {{9, 7, 8, 9, 5, 7, 10, 1, 2, X, X, X, X, X, X, X}},
5472 {{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, X, X, X, X}},
5473 {{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, X, X, X, X}},
5474 {{2, 10, 5, 2, 5, 3, 3, 5, 7, X, X, X, X, X, X, X}},
5475 {{7, 9, 5, 7, 8, 9, 3, 11, 2, X, X, X, X, X, X, X}},
5476 {{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, X, X, X, X}},
5477 {{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, X, X, X, X}},
5478 {{11, 2, 1, 11, 1, 7, 7, 1, 5, X, X, X, X, X, X, X}},
5479 {{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, X, X, X, X}},
5480 {{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, X}},
5481 {{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, X}},
5482 {{11, 10, 5, 7, 11, 5, X, X, X, X, X, X, X, X, X, X}},
5483 {{10, 6, 5, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5484 {{0, 8, 3, 5, 10, 6, X, X, X, X, X, X, X, X, X, X}},
5485 {{9, 0, 1, 5, 10, 6, X, X, X, X, X, X, X, X, X, X}},
5486 {{1, 8, 3, 1, 9, 8, 5, 10, 6, X, X, X, X, X, X, X}},
5487 {{1, 6, 5, 2, 6, 1, X, X, X, X, X, X, X, X, X, X}},
5488 {{1, 6, 5, 1, 2, 6, 3, 0, 8, X, X, X, X, X, X, X}},
5489 {{9, 6, 5, 9, 0, 6, 0, 2, 6, X, X, X, X, X, X, X}},
5490 {{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, X, X, X, X}},
5491 {{2, 3, 11, 10, 6, 5, X, X, X, X, X, X, X, X, X, X}},
5492 {{11, 0, 8, 11, 2, 0, 10, 6, 5, X, X, X, X, X, X, X}},
5493 {{0, 1, 9, 2, 3, 11, 5, 10, 6, X, X, X, X, X, X, X}},
5494 {{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, X, X, X, X}},
5495 {{6, 3, 11, 6, 5, 3, 5, 1, 3, X, X, X, X, X, X, X}},
5496 {{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, X, X, X, X}},
5497 {{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, X, X, X, X}},
5498 {{6, 5, 9, 6, 9, 11, 11, 9, 8, X, X, X, X, X, X, X}},
5499 {{5, 10, 6, 4, 7, 8, X, X, X, X, X, X, X, X, X, X}},
5500 {{4, 3, 0, 4, 7, 3, 6, 5, 10, X, X, X, X, X, X, X}},
5501 {{1, 9, 0, 5, 10, 6, 8, 4, 7, X, X, X, X, X, X, X}},
5502 {{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, X, X, X, X}},
5503 {{6, 1, 2, 6, 5, 1, 4, 7, 8, X, X, X, X, X, X, X}},
5504 {{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, X, X, X, X}},
5505 {{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, X, X, X, X}},
5506 {{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, X}},
5507 {{3, 11, 2, 7, 8, 4, 10, 6, 5, X, X, X, X, X, X, X}},
5508 {{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, X, X, X, X}},
5509 {{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, X, X, X, X}},
5510 {{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, X}},
5511 {{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, X, X, X, X}},
5512 {{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, X}},
5513 {{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, X}},
5514 {{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, X, X, X, X}},
5515 {{10, 4, 9, 6, 4, 10, X, X, X, X, X, X, X, X, X, X}},
5516 {{4, 10, 6, 4, 9, 10, 0, 8, 3, X, X, X, X, X, X, X}},
5517 {{10, 0, 1, 10, 6, 0, 6, 4, 0, X, X, X, X, X, X, X}},
5518 {{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, X, X, X, X}},
5519 {{1, 4, 9, 1, 2, 4, 2, 6, 4, X, X, X, X, X, X, X}},
5520 {{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, X, X, X, X}},
5521 {{0, 2, 4, 4, 2, 6, X, X, X, X, X, X, X, X, X, X}},
5522 {{8, 3, 2, 8, 2, 4, 4, 2, 6, X, X, X, X, X, X, X}},
5523 {{10, 4, 9, 10, 6, 4, 11, 2, 3, X, X, X, X, X, X, X}},
5524 {{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, X, X, X, X}},
5525 {{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, X, X, X, X}},
5526 {{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, X}},
5527 {{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, X, X, X, X}},
5528 {{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, X}},
5529 {{3, 11, 6, 3, 6, 0, 0, 6, 4, X, X, X, X, X, X, X}},
5530 {{6, 4, 8, 11, 6, 8, X, X, X, X, X, X, X, X, X, X}},
5531 {{7, 10, 6, 7, 8, 10, 8, 9, 10, X, X, X, X, X, X, X}},
5532 {{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, X, X, X, X}},
5533 {{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, X, X, X, X}},
5534 {{10, 6, 7, 10, 7, 1, 1, 7, 3, X, X, X, X, X, X, X}},
5535 {{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, X, X, X, X}},
5536 {{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, X}},
5537 {{7, 8, 0, 7, 0, 6, 6, 0, 2, X, X, X, X, X, X, X}},
5538 {{7, 3, 2, 6, 7, 2, X, X, X, X, X, X, X, X, X, X}},
5539 {{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, X, X, X, X}},
5540 {{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, X}},
5541 {{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, X}},
5542 {{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, X, X, X, X}},
5543 {{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, X}},
5544 {{0, 9, 1, 11, 6, 7, X, X, X, X, X, X, X, X, X, X}},
5545 {{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, X, X, X, X}},
5546 {{7, 11, 6, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5547 {{7, 6, 11, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5548 {{3, 0, 8, 11, 7, 6, X, X, X, X, X, X, X, X, X, X}},
5549 {{0, 1, 9, 11, 7, 6, X, X, X, X, X, X, X, X, X, X}},
5550 {{8, 1, 9, 8, 3, 1, 11, 7, 6, X, X, X, X, X, X, X}},
5551 {{10, 1, 2, 6, 11, 7, X, X, X, X, X, X, X, X, X, X}},
5552 {{1, 2, 10, 3, 0, 8, 6, 11, 7, X, X, X, X, X, X, X}},
5553 {{2, 9, 0, 2, 10, 9, 6, 11, 7, X, X, X, X, X, X, X}},
5554 {{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, X, X, X, X}},
5555 {{7, 2, 3, 6, 2, 7, X, X, X, X, X, X, X, X, X, X}},
5556 {{7, 0, 8, 7, 6, 0, 6, 2, 0, X, X, X, X, X, X, X}},
5557 {{2, 7, 6, 2, 3, 7, 0, 1, 9, X, X, X, X, X, X, X}},
5558 {{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, X, X, X, X}},
5559 {{10, 7, 6, 10, 1, 7, 1, 3, 7, X, X, X, X, X, X, X}},
5560 {{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, X, X, X, X}},
5561 {{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, X, X, X, X}},
5562 {{7, 6, 10, 7, 10, 8, 8, 10, 9, X, X, X, X, X, X, X}},
5563 {{6, 8, 4, 11, 8, 6, X, X, X, X, X, X, X, X, X, X}},
5564 {{3, 6, 11, 3, 0, 6, 0, 4, 6, X, X, X, X, X, X, X}},
5565 {{8, 6, 11, 8, 4, 6, 9, 0, 1, X, X, X, X, X, X, X}},
5566 {{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, X, X, X, X}},
5567 {{6, 8, 4, 6, 11, 8, 2, 10, 1, X, X, X, X, X, X, X}},
5568 {{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, X, X, X, X}},
5569 {{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, X, X, X, X}},
5570 {{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, X}},
5571 {{8, 2, 3, 8, 4, 2, 4, 6, 2, X, X, X, X, X, X, X}},
5572 {{0, 4, 2, 4, 6, 2, X, X, X, X, X, X, X, X, X, X}},
5573 {{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, X, X, X, X}},
5574 {{1, 9, 4, 1, 4, 2, 2, 4, 6, X, X, X, X, X, X, X}},
5575 {{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, X, X, X, X}},
5576 {{10, 1, 0, 10, 0, 6, 6, 0, 4, X, X, X, X, X, X, X}},
5577 {{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, X}},
5578 {{10, 9, 4, 6, 10, 4, X, X, X, X, X, X, X, X, X, X}},
5579 {{4, 9, 5, 7, 6, 11, X, X, X, X, X, X, X, X, X, X}},
5580 {{0, 8, 3, 4, 9, 5, 11, 7, 6, X, X, X, X, X, X, X}},
5581 {{5, 0, 1, 5, 4, 0, 7, 6, 11, X, X, X, X, X, X, X}},
5582 {{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, X, X, X, X}},
5583 {{9, 5, 4, 10, 1, 2, 7, 6, 11, X, X, X, X, X, X, X}},
5584 {{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, X, X, X, X}},
5585 {{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, X, X, X, X}},
5586 {{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, X}},
5587 {{7, 2, 3, 7, 6, 2, 5, 4, 9, X, X, X, X, X, X, X}},
5588 {{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, X, X, X, X}},
5589 {{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, X, X, X, X}},
5590 {{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, X}},
5591 {{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, X, X, X, X}},
5592 {{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, X}},
5593 {{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, X}},
5594 {{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, X, X, X, X}},
5595 {{6, 9, 5, 6, 11, 9, 11, 8, 9, X, X, X, X, X, X, X}},
5596 {{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, X, X, X, X}},
5597 {{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, X, X, X, X}},
5598 {{6, 11, 3, 6, 3, 5, 5, 3, 1, X, X, X, X, X, X, X}},
5599 {{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, X, X, X, X}},
5600 {{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, X}},
5601 {{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, X}},
5602 {{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, X, X, X, X}},
5603 {{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, X, X, X, X}},
5604 {{9, 5, 6, 9, 6, 0, 0, 6, 2, X, X, X, X, X, X, X}},
5605 {{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, X}},
5606 {{1, 5, 6, 2, 1, 6, X, X, X, X, X, X, X, X, X, X}},
5607 {{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, X}},
5608 {{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, X, X, X, X}},
5609 {{0, 3, 8, 5, 6, 10, X, X, X, X, X, X, X, X, X, X}},
5610 {{10, 5, 6, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5611 {{11, 5, 10, 7, 5, 11, X, X, X, X, X, X, X, X, X, X}},
5612 {{11, 5, 10, 11, 7, 5, 8, 3, 0, X, X, X, X, X, X, X}},
5613 {{5, 11, 7, 5, 10, 11, 1, 9, 0, X, X, X, X, X, X, X}},
5614 {{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, X, X, X, X}},
5615 {{11, 1, 2, 11, 7, 1, 7, 5, 1, X, X, X, X, X, X, X}},
5616 {{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, X, X, X, X}},
5617 {{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, X, X, X, X}},
5618 {{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, X}},
5619 {{2, 5, 10, 2, 3, 5, 3, 7, 5, X, X, X, X, X, X, X}},
5620 {{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, X, X, X, X}},
5621 {{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, X, X, X, X}},
5622 {{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, X}},
5623 {{1, 3, 5, 3, 7, 5, X, X, X, X, X, X, X, X, X, X}},
5624 {{0, 8, 7, 0, 7, 1, 1, 7, 5, X, X, X, X, X, X, X}},
5625 {{9, 0, 3, 9, 3, 5, 5, 3, 7, X, X, X, X, X, X, X}},
5626 {{9, 8, 7, 5, 9, 7, X, X, X, X, X, X, X, X, X, X}},
5627 {{5, 8, 4, 5, 10, 8, 10, 11, 8, X, X, X, X, X, X, X}},
5628 {{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, X, X, X, X}},
5629 {{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, X, X, X, X}},
5630 {{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, X}},
5631 {{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, X, X, X, X}},
5632 {{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, X}},
5633 {{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, X}},
5634 {{9, 4, 5, 2, 11, 3, X, X, X, X, X, X, X, X, X, X}},
5635 {{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, X, X, X, X}},
5636 {{5, 10, 2, 5, 2, 4, 4, 2, 0, X, X, X, X, X, X, X}},
5637 {{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, X}},
5638 {{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, X, X, X, X}},
5639 {{8, 4, 5, 8, 5, 3, 3, 5, 1, X, X, X, X, X, X, X}},
5640 {{0, 4, 5, 1, 0, 5, X, X, X, X, X, X, X, X, X, X}},
5641 {{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, X, X, X, X}},
5642 {{9, 4, 5, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5643 {{4, 11, 7, 4, 9, 11, 9, 10, 11, X, X, X, X, X, X, X}},
5644 {{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, X, X, X, X}},
5645 {{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, X, X, X, X}},
5646 {{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, X}},
5647 {{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, X, X, X, X}},
5648 {{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, X}},
5649 {{11, 7, 4, 11, 4, 2, 2, 4, 0, X, X, X, X, X, X, X}},
5650 {{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, X, X, X, X}},
5651 {{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, X, X, X, X}},
5652 {{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, X}},
5653 {{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, X}},
5654 {{1, 10, 2, 8, 7, 4, X, X, X, X, X, X, X, X, X, X}},
5655 {{4, 9, 1, 4, 1, 7, 7, 1, 3, X, X, X, X, X, X, X}},
5656 {{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, X, X, X, X}},
5657 {{4, 0, 3, 7, 4, 3, X, X, X, X, X, X, X, X, X, X}},
5658 {{4, 8, 7, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5659 {{9, 10, 8, 10, 11, 8, X, X, X, X, X, X, X, X, X, X}},
5660 {{3, 0, 9, 3, 9, 11, 11, 9, 10, X, X, X, X, X, X, X}},
5661 {{0, 1, 10, 0, 10, 8, 8, 10, 11, X, X, X, X, X, X, X}},
5662 {{3, 1, 10, 11, 3, 10, X, X, X, X, X, X, X, X, X, X}},
5663 {{1, 2, 11, 1, 11, 9, 9, 11, 8, X, X, X, X, X, X, X}},
5664 {{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, X, X, X, X}},
5665 {{0, 2, 11, 8, 0, 11, X, X, X, X, X, X, X, X, X, X}},
5666 {{3, 2, 11, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5667 {{2, 3, 8, 2, 8, 10, 10, 8, 9, X, X, X, X, X, X, X}},
5668 {{9, 10, 2, 0, 9, 2, X, X, X, X, X, X, X, X, X, X}},
5669 {{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, X, X, X, X}},
5670 {{1, 10, 2, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5671 {{1, 3, 8, 9, 1, 8, X, X, X, X, X, X, X, X, X, X}},
5672 {{0, 9, 1, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5673 {{0, 3, 8, X, X, X, X, X, X, X, X, X, X, X, X, X}},
5674 {{X, X, X, X, X, X, X, X, X, X, X, X, X, X, X, X}}}};
5699 line_to_vertex_table,
5707 write_back_cell_data);
5714 #include "grid_tools.inst"
void distribute(VectorType &vec) const
std::pair< std::vector< std::pair< int, int > >, std::vector< int > > query(const QueryType &queries)
DistributedTree(const MPI_Comm comm, const std::vector< BoundingBox< dim, Number >> &bounding_boxes)
BoundingBox< spacedim, Number > create_extended(const Number amount) const
BoundingBox< spacedim, Number > create_extended_relative(const Number relative_amount) const
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
types::global_dof_index n_dofs() const
Abstract base class for mapping classes.
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
virtual boost::container::small_vector< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr numbers::NumberTraits< Number >::real_type square() const
constexpr numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
Quadrature< spacedim > compute_affine_transformation(const std::array< Point< spacedim >, dim+1 > &vertices) const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
numbers::NumberTraits< Number >::real_type norm() const
virtual std::vector< types::boundary_id > get_boundary_ids() const
const std::vector< bool > & get_used_vertices() const
virtual MPI_Comm get_communicator() const
unsigned int n_quads() const
void load_user_indices(const std::vector< unsigned int > &v)
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
bool all_reference_cells_are_hyper_cube() const
face_iterator end_face() const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_lines() const
unsigned int n_raw_faces() const
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
cell_iterator end() const
virtual bool has_hanging_nodes() const
bool vertex_used(const unsigned int index) const
virtual void execute_coarsening_and_refinement()
const std::map< std::pair< cell_iterator, unsigned int >, std::pair< std::pair< cell_iterator, unsigned int >, std::bitset< 3 > > > & get_periodic_face_map() const
Triangulation< dim, spacedim > & get_triangulation()
unsigned int n_vertices() const
void save_user_indices(std::vector< unsigned int > &v) const
active_face_iterator begin_active_face() const
const std::vector< Point< spacedim > > & get_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) 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
IteratorRange< active_cell_iterator > active_cell_iterators() const
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsMPI()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcNeedsCGAL()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
PackagedOperation< Range > constrained_right_hand_side(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop, const Range &right_hand_side)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
LinearOperator< Range, Domain, Payload > constrained_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop)
@ update_values
Shape function values.
@ update_quadrature_points
Transformed quadrature points.
virtual std::vector< types::manifold_id > get_manifold_ids() const
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
void map_boundary_to_manifold_ids(const std::vector< types::boundary_id > &src_boundary_ids, const std::vector< types::manifold_id > &dst_manifold_ids, Triangulation< dim, spacedim > &tria, const std::vector< types::boundary_id > &reset_boundary_ids={})
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
void assign_co_dimensional_manifold_indicators(Triangulation< dim, spacedim > &tria, const std::function< types::manifold_id(const std::set< types::manifold_id > &)> &disambiguation_function=[](const std::set< types::manifold_id > &manifold_ids) { if(manifold_ids.size()==1) return *manifold_ids.begin();else return numbers::flat_manifold_id;}, bool overwrite_only_flat_manifold_ids=true)
void consistently_order_cells(std::vector< CellData< dim >> &cells)
Task< RT > new_task(const std::function< RT()> &function)
Expression fabs(const Expression &x)
Expression floor(const Expression &x)
Expression acos(const Expression &x)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
@ valid
Iterator points to a valid object.
void create_laplace_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, SparseMatrixType &matrix, const Function< spacedim, typename SparseMatrixType::value_type > *const a=nullptr, const AffineConstraints< typename SparseMatrixType::value_type > &constraints=AffineConstraints< typename SparseMatrixType::value_type >())
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)