18 #include <deal.II/base/mpi.templates.h>
23 #ifdef DEAL_II_WITH_CGAL
68 #include <boost/random/mersenne_twister.hpp>
69 #include <boost/random/uniform_real_distribution.hpp>
79 #include <unordered_map>
86 template <
int dim,
int spacedim>
96 #if defined(DEAL_II_WITH_P4EST) && defined(DEBUG)
111 std::vector<bool> boundary_vertices(
vertices.size(),
false);
117 for (; cell != endc; ++cell)
118 for (
const unsigned int face : cell->face_indices())
119 if (cell->face(face)->at_boundary())
120 for (
unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
121 boundary_vertices[cell->face(face)->vertex_index(i)] =
true;
125 double max_distance_sqr = 0;
126 std::vector<bool>::const_iterator pi = boundary_vertices.begin();
127 const unsigned int N = boundary_vertices.size();
128 for (
unsigned int i = 0; i <
N; ++i, ++pi)
130 std::vector<bool>::const_iterator pj = pi + 1;
131 for (
unsigned int j = i + 1; j <
N; ++j, ++pj)
132 if ((*pi ==
true) && (*pj ==
true) &&
137 return std::sqrt(max_distance_sqr);
142 template <
int dim,
int spacedim>
151 reference_cell.template get_default_linear_mapping<dim, spacedim>());
156 template <
int dim,
int spacedim>
162 unsigned int mapping_degree = 1;
164 mapping_degree = p->get_degree();
165 else if (
const auto *p =
167 mapping_degree = p->get_degree();
174 reference_cell.template get_gauss_type_quadrature<dim>(mapping_degree +
176 const unsigned int n_q_points = quadrature_formula.
size();
187 double local_volume = 0;
190 for (
const auto &cell :
triangulation.active_cell_iterators())
191 if (cell->is_locally_owned())
194 for (
unsigned int q = 0; q < n_q_points; ++q)
195 local_volume += fe_values.
JxW(q);
198 double global_volume = 0;
200 #ifdef DEAL_II_WITH_MPI
208 global_volume = local_volume;
210 return global_volume;
232 struct TransformR2UAffine
247 [1] = {{-1.000000}, {1.000000}};
251 {1.000000, 0.000000};
262 [2] = {{-0.500000, -0.500000},
263 {0.500000, -0.500000},
264 {-0.500000, 0.500000},
265 {0.500000, 0.500000}};
275 {0.750000, 0.250000, 0.250000, -0.250000};
281 {-0.250000, -0.250000, -0.250000},
282 {0.250000, -0.250000, -0.250000},
283 {-0.250000, 0.250000, -0.250000},
284 {0.250000, 0.250000, -0.250000},
285 {-0.250000, -0.250000, 0.250000},
286 {0.250000, -0.250000, 0.250000},
287 {-0.250000, 0.250000, 0.250000},
288 {0.250000, 0.250000, 0.250000}
307 template <
int dim,
int spacedim>
316 for (
unsigned int d = 0;
d < spacedim; ++
d)
317 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
318 for (
unsigned int e = 0;
e < dim; ++
e)
323 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
326 return std::make_pair(
A,
b);
343 for (
const auto &cell :
triangulation.active_cell_iterators())
345 if (cell->is_locally_owned())
347 double aspect_ratio_cell = 0.0;
352 for (
unsigned int q = 0; q < quadrature.
size(); ++q)
363 aspect_ratio_cell = std::numeric_limits<double>::infinity();
368 for (
unsigned int i = 0; i < dim; ++i)
369 for (
unsigned int j = 0; j < dim; ++j)
370 J(i, j) = jacobian[i][j];
376 double const ar = max_sv / min_sv;
383 aspect_ratio_cell =
std::max(aspect_ratio_cell, ar);
388 aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
392 return aspect_ratio_vector;
413 template <
int dim,
int spacedim>
419 const auto predicate = [](
const iterator &) {
return true; };
422 tria, std::function<
bool(
const iterator &)>(predicate));
448 template <
int structdim>
472 "Two CellData objects with equal vertices must "
473 "have the same material/boundary ids and manifold "
498 template <
class FaceIteratorType>
502 CellData<dim - 1> face_cell_data(face->n_vertices());
503 for (
unsigned int vertex_n = 0; vertex_n < face->n_vertices();
505 face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
506 face_cell_data.boundary_id = face->boundary_id();
507 face_cell_data.manifold_id = face->manifold_id();
509 face_data.insert(std::move(face_cell_data));
537 template <
class FaceIteratorType>
552 template <
int dim,
int spacedim>
554 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>,
SubCellData>
558 ExcMessage(
"The input triangulation must be non-empty."));
560 std::vector<Point<spacedim>>
vertices;
561 std::vector<CellData<dim>> cells;
563 unsigned int max_level_0_vertex_n = 0;
565 for (
const unsigned int cell_vertex_n : cell->vertex_indices())
566 max_level_0_vertex_n =
567 std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
568 vertices.resize(max_level_0_vertex_n + 1);
578 for (
const unsigned int cell_vertex_n : cell->vertex_indices())
582 vertices[cell->vertex_index(cell_vertex_n)] =
583 cell->vertex(cell_vertex_n);
585 cell->vertex_index(cell_vertex_n);
589 cells.push_back(cell_data);
594 for (
const unsigned int face_n : cell->face_indices())
597 const auto face = cell->face(face_n);
606 for (
unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
608 const auto line = cell->line(line_n);
614 for (
unsigned int vertex_n : line->vertex_indices())
616 line->vertex_index(vertex_n);
619 line_data.insert(std::move(line_cell_data));
628 std::vector<bool> used_vertices(
vertices.size());
630 for (
const auto v : cell_data.vertices)
631 used_vertices[v] =
true;
632 Assert(std::find(used_vertices.begin(), used_vertices.end(),
false) ==
634 ExcMessage(
"The level zero vertices should form a contiguous "
642 for (
const CellData<1> &face_line_data : line_data)
645 return std::tuple<std::vector<Point<spacedim>>,
646 std::vector<CellData<dim>>,
649 std::move(subcell_data));
654 template <
int dim,
int spacedim>
663 "Invalid SubCellData supplied according to ::check_consistency(). "
664 "This is caused by data containing objects for the wrong dimension."));
667 std::vector<bool> vertex_used(
vertices.size(),
false);
668 for (
unsigned int c = 0; c < cells.size(); ++c)
669 for (
unsigned int v = 0; v < cells[c].vertices.size(); ++v)
672 ExcMessage(
"Invalid vertex index encountered! cells[" +
676 " is invalid, because only " +
678 " vertices were supplied."));
679 vertex_used[cells[c].vertices[v]] =
true;
686 std::vector<unsigned int> new_vertex_numbers(
vertices.size(),
688 unsigned int next_free_number = 0;
689 for (
unsigned int i = 0; i <
vertices.size(); ++i)
690 if (vertex_used[i] ==
true)
692 new_vertex_numbers[i] = next_free_number;
697 for (
unsigned int c = 0; c < cells.size(); ++c)
699 v = new_vertex_numbers[v];
704 for (
unsigned int v = 0;
709 new_vertex_numbers.size(),
711 "Invalid vertex index in subcelldata.boundary_lines. "
712 "subcelldata.boundary_lines[" +
717 " is invalid, because only " +
719 " vertices were supplied."));
726 for (
unsigned int v = 0;
731 new_vertex_numbers.size(),
733 "Invalid vertex index in subcelldata.boundary_quads. "
734 "subcelldata.boundary_quads[" +
739 " is invalid, because only " +
741 " vertices were supplied."));
749 std::vector<Point<spacedim>> tmp;
750 tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(),
true));
751 for (
unsigned int v = 0; v <
vertices.size(); ++v)
752 if (vertex_used[v] ==
true)
759 template <
int dim,
int spacedim>
764 std::vector<unsigned int> & considered_vertices,
771 std::vector<unsigned int> new_vertex_numbers(
vertices.size());
772 std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
775 if (considered_vertices.size() == 0)
776 considered_vertices = new_vertex_numbers;
791 unsigned int longest_coordinate_direction = 0;
792 double longest_coordinate_length = bbox.
side_length(0);
793 for (
unsigned int d = 1;
d < spacedim; ++
d)
796 if (longest_coordinate_length < coordinate_length)
798 longest_coordinate_length = coordinate_length;
799 longest_coordinate_direction =
d;
805 std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
806 sorted_vertices.reserve(
vertices.size());
807 for (
const unsigned int vertex_n : considered_vertices)
810 sorted_vertices.emplace_back(vertex_n,
vertices[vertex_n]);
812 std::sort(sorted_vertices.begin(),
813 sorted_vertices.end(),
816 return a.second[longest_coordinate_direction] <
817 b.second[longest_coordinate_direction];
822 for (
unsigned int d = 0;
d < spacedim; ++
d)
823 if (std::abs(a[
d] -
b[
d]) > tol)
830 auto range_start = sorted_vertices.begin();
831 while (range_start != sorted_vertices.end())
833 auto range_end = range_start + 1;
834 while (range_end != sorted_vertices.end() &&
835 std::abs(range_end->second[longest_coordinate_direction] -
836 range_start->second[longest_coordinate_direction]) <
842 std::sort(range_start,
846 return a.first <
b.first;
855 for (
auto reference = range_start; reference != range_end; ++reference)
858 for (
auto it = reference + 1; it != range_end; ++it)
860 if (within_tolerance(reference->second, it->second))
862 new_vertex_numbers[it->first] = reference->first;
868 range_start = range_end;
875 for (
auto &cell : cells)
876 for (
auto &vertex_index : cell.vertices)
877 vertex_index = new_vertex_numbers[vertex_index];
879 for (
auto &vertex_index : quad.vertices)
880 vertex_index = new_vertex_numbers[vertex_index];
882 for (
auto &vertex_index : line.vertices)
883 vertex_index = new_vertex_numbers[vertex_index];
903 for (
unsigned int i = 0; i <
vertices.size(); ++i)
904 map_point_to_local_vertex_index[
vertices[i]] = i;
907 if (map_point_to_local_vertex_index.size() ==
vertices.size())
911 vertices.resize(map_point_to_local_vertex_index.size());
914 for (
const auto &p : map_point_to_local_vertex_index)
921 template <
int dim,
int spacedim>
932 if (dim == 2 && spacedim == 3)
935 std::size_t n_negative_cells = 0;
936 std::size_t cell_no = 0;
937 for (
auto &cell : cells)
953 std::swap(cell.vertices[1], cell.vertices[2]);
958 std::swap(cell.vertices[0], cell.vertices[2]);
959 std::swap(cell.vertices[1], cell.vertices[3]);
960 std::swap(cell.vertices[4], cell.vertices[6]);
961 std::swap(cell.vertices[5], cell.vertices[7]);
969 std::swap(cell.vertices[cell.vertices.size() - 2],
970 cell.vertices[cell.vertices.size() - 1]);
975 std::swap(cell.vertices[0], cell.vertices[3]);
976 std::swap(cell.vertices[1], cell.vertices[4]);
977 std::swap(cell.vertices[2], cell.vertices[5]);
983 std::swap(cell.vertices[2], cell.vertices[3]);
997 return n_negative_cells;
1001 template <
int dim,
int spacedim>
1007 const std::size_t n_negative_cells =
1016 AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
1019 "This function assumes that either all cells have positive "
1020 "volume, or that all cells have been specified in an "
1021 "inverted vertex order so that their volume is negative. "
1022 "(In the latter case, this class automatically inverts "
1023 "every cell.) However, the mesh you have specified "
1024 "appears to have both cells with positive and cells with "
1025 "negative volume. You need to check your mesh which "
1026 "cells these are and how they got there.\n"
1027 "As a hint, of the total ") +
1030 " appear to have a negative volume."));
1048 CheapEdge(
const unsigned int v0,
const unsigned int v1)
1060 return ((
v0 <
e.v0) || ((
v0 ==
e.v0) && (
v1 <
e.v1)));
1083 std::set<CheapEdge> edges;
1085 for (
typename std::vector<
CellData<dim>>::const_iterator c =
1095 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1097 const CheapEdge reverse_edge(
1100 if (edges.find(reverse_edge) != edges.end())
1105 const CheapEdge correct_edge(
1108 edges.insert(correct_edge);
1124 struct ParallelEdges
1138 static const unsigned int
1203 class AdjacentCells;
1211 class AdjacentCells<2>
1218 using const_iterator =
const AdjacentCell *;
1229 push_back(
const AdjacentCell &adjacent_cell)
1293 class AdjacentCells<3> :
public std::vector<AdjacentCell>
1315 Edge(
const CellData<dim> &cell,
const unsigned int edge_number)
1366 enum OrientationStatus
1396 Cell(
const CellData<dim> &c,
const std::vector<Edge<dim>> &edge_list)
1403 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1405 const Edge<dim>
e(c,
l);
1441 class EdgeDeltaSet<2>
1447 using const_iterator =
const unsigned int *;
1473 insert(
const unsigned int edge_index)
1534 class EdgeDeltaSet<3> :
public std::set<unsigned int>
1544 std::vector<Edge<dim>>
1551 std::vector<Edge<dim>> edge_list;
1553 for (
unsigned int i = 0; i < cells.size(); ++i)
1554 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1555 edge_list.emplace_back(cells[i],
l);
1558 std::sort(edge_list.begin(), edge_list.end());
1559 edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
1572 std::vector<Cell<dim>>
1573 build_cells_and_connect_edges(
const std::vector<
CellData<dim>> &cells,
1574 std::vector<Edge<dim>> & edges)
1576 std::vector<Cell<dim>> cell_list;
1577 cell_list.reserve(cells.size());
1578 for (
unsigned int i = 0; i < cells.size(); ++i)
1582 cell_list.emplace_back(cells[i], edges);
1586 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1587 edges[cell_list.back().edge_indices[
l]].adjacent_cells.push_back(
1588 AdjacentCell(i,
l));
1603 get_next_unoriented_cell(
const std::vector<Cell<dim>> &cells,
1604 const std::vector<Edge<dim>> &edges,
1605 const unsigned int current_cell)
1607 for (
unsigned int c = current_cell; c < cells.size(); ++c)
1608 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1610 Edge<dim>::not_oriented)
1625 orient_one_set_of_parallel_edges(
const std::vector<Cell<dim>> &cells,
1626 std::vector<Edge<dim>> & edges,
1627 const unsigned int cell,
1628 const unsigned int local_edge)
1655 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1656 (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
1672 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1673 (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
1682 EdgeDeltaSet<dim> Delta_k;
1683 EdgeDeltaSet<dim> Delta_k_minus_1;
1684 Delta_k_minus_1.insert(cells[cell].
edge_indices[local_edge]);
1686 while (Delta_k_minus_1.begin() !=
1687 Delta_k_minus_1.end())
1691 for (
typename EdgeDeltaSet<dim>::const_iterator delta =
1692 Delta_k_minus_1.begin();
1693 delta != Delta_k_minus_1.end();
1697 Edge<dim>::not_oriented,
1701 for (
typename AdjacentCells<dim>::const_iterator adjacent_cell =
1703 adjacent_cell != edges[*delta].adjacent_cells.end();
1706 const unsigned int K = adjacent_cell->cell_index;
1707 const unsigned int delta_is_edge_in_K =
1708 adjacent_cell->edge_within_cell;
1712 const unsigned int first_edge_vertex =
1713 (edges[*delta].orientation_status == Edge<dim>::forward ?
1714 edges[*delta].vertex_indices[0] :
1715 edges[*delta].vertex_indices[1]);
1716 const unsigned int first_edge_vertex_in_K =
1719 delta_is_edge_in_K, 0)];
1721 first_edge_vertex == first_edge_vertex_in_K ||
1722 first_edge_vertex ==
1724 dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
1729 for (
unsigned int o_e = 0;
1736 const unsigned int opposite_edge =
1737 cells[
K].edge_indices[ParallelEdges<
1739 const unsigned int first_opposite_edge_vertex =
1740 cells[
K].vertex_indices
1744 (first_edge_vertex == first_edge_vertex_in_K ? 0 :
1751 const typename Edge<dim>::OrientationStatus
1752 opposite_edge_orientation =
1753 (edges[opposite_edge].vertex_indices[0] ==
1754 first_opposite_edge_vertex ?
1755 Edge<dim>::forward :
1756 Edge<dim>::backward);
1761 Edge<dim>::not_oriented)
1765 edges[opposite_edge].orientation_status =
1766 opposite_edge_orientation;
1767 Delta_k.insert(opposite_edge);
1783 opposite_edge_orientation,
1789 opposite_edge_orientation)
1802 Delta_k_minus_1 = Delta_k;
1816 rotate_cell(
const std::vector<Cell<dim>> &cell_list,
1817 const std::vector<Edge<dim>> &edge_list,
1824 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++
e)
1831 starting_vertex_of_edge[
e] =
1835 starting_vertex_of_edge[
e] =
1851 if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
1852 (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
1853 origin_vertex_of_cell = starting_vertex_of_edge[0];
1854 else if ((starting_vertex_of_edge[1] ==
1855 starting_vertex_of_edge[2]) ||
1856 (starting_vertex_of_edge[1] ==
1857 starting_vertex_of_edge[3]))
1858 origin_vertex_of_cell = starting_vertex_of_edge[1];
1870 for (origin_vertex_of_cell = 0;
1871 origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
1872 ++origin_vertex_of_cell)
1873 if (std::count(starting_vertex_of_edge,
1874 starting_vertex_of_edge +
1879 Assert(origin_vertex_of_cell <
1904 const unsigned int tmp = raw_cells[
cell_index].vertices[0];
1927 static const unsigned int cube_permutations[8][8] = {
1928 {0, 1, 2, 3, 4, 5, 6, 7},
1929 {1, 5, 3, 7, 0, 4, 2, 6},
1930 {2, 6, 0, 4, 3, 7, 1, 5},
1931 {3, 2, 1, 0, 7, 6, 5, 4},
1932 {4, 0, 6, 2, 5, 1, 7, 3},
1933 {5, 4, 7, 6, 1, 0, 3, 2},
1934 {6, 7, 4, 5, 2, 3, 0, 1},
1935 {7, 3, 5, 1, 6, 2, 4, 0}};
1940 temp_vertex_indices[v] =
1942 .vertices[cube_permutations[origin_vertex_of_cell][v]];
1944 raw_cells[
cell_index].vertices[v] = temp_vertex_indices[v];
1968 std::vector<Edge<dim>> edge_list = build_edges(cells);
1969 std::vector<Cell<dim>> cell_list =
1970 build_cells_and_connect_edges(cells, edge_list);
1974 unsigned int next_cell_with_unoriented_edge = 0;
1975 while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
1976 cell_list, edge_list, next_cell_with_unoriented_edge)) !=
1990 for (
unsigned int l = 0;
l < dim; ++
l)
1991 if (edge_list[cell_list[next_cell_with_unoriented_edge]
1994 orient_one_set_of_parallel_edges(
1997 next_cell_with_unoriented_edge,
2002 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
2003 Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
2012 for (
unsigned int c = 0; c < cells.size(); ++c)
2013 rotate_cell(cell_list, edge_list, c, cells);
2028 Assert(cells.size() != 0,
2030 "List of elements to orient must have at least one cell"));
2038 if (!is_consistent(cells))
2057 template <
int spacedim>
2115 template <
int spacedim>
2134 template <
int dim,
int spacedim>
2144 template <
int dim,
int spacedim>
2153 "GridTools::rotate() is only available for spacedim = 2."));
2188 const unsigned int axis,
2200 template <
int dim,
int spacedim>
2222 const unsigned int n_dofs = S.
n();
2234 const auto constrained_rhs =
2236 solver.
solve(SF, u, constrained_rhs, prec);
2249 const bool solve_for_absolute_positions)
2275 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
2284 std::array<AffineConstraints<double>, dim> constraints;
2285 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
2293 for (
const unsigned int vertex_no : cell->vertex_indices())
2295 const unsigned int vertex_index = cell->vertex_index(vertex_no);
2296 const Point<dim> & vertex_point = cell->vertex(vertex_no);
2298 const typename std::map<unsigned int, Point<dim>>::const_iterator
2299 map_iter = new_points.find(vertex_index);
2301 if (map_iter != map_end)
2302 for (
unsigned int i = 0; i < dim; ++i)
2304 constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
2305 constraints[i].set_inhomogeneity(
2306 cell->vertex_dof_index(vertex_no, 0),
2307 (solve_for_absolute_positions ?
2308 map_iter->second(i) :
2309 map_iter->second(i) - vertex_point[i]));
2314 for (
unsigned int i = 0; i < dim; ++i)
2315 constraints[i].close();
2319 for (
unsigned int i = 0; i < dim; ++i)
2324 for (
unsigned int i = 0; i < dim; ++i)
2331 std::vector<bool> vertex_touched(
triangulation.n_vertices(),
false);
2333 for (
const unsigned int vertex_no : cell->vertex_indices())
2334 if (vertex_touched[cell->vertex_index(vertex_no)] ==
false)
2339 cell->vertex_dof_index(vertex_no, 0);
2340 for (
unsigned int i = 0; i < dim; ++i)
2341 if (solve_for_absolute_positions)
2342 v(i) = us[i](dof_index);
2344 v(i) += us[i](dof_index);
2346 vertex_touched[cell->vertex_index(vertex_no)] =
true;
2350 template <
int dim,
int spacedim>
2351 std::map<unsigned int, Point<spacedim>>
2354 std::map<unsigned int, Point<spacedim>> vertex_map;
2358 for (; cell != endc; ++cell)
2360 for (
unsigned int i : cell->face_indices())
2364 if (face->at_boundary())
2366 for (
unsigned j = 0; j < face->
n_vertices(); ++j)
2369 const unsigned int vertex_index = face->vertex_index(j);
2370 vertex_map[vertex_index] = vertex;
2382 template <
int dim,
int spacedim>
2386 const bool keep_boundary,
2387 const unsigned int seed)
2404 double almost_infinite_length = 0;
2409 almost_infinite_length += cell->diameter();
2411 std::vector<double> minimal_length(
triangulation.n_vertices(),
2412 almost_infinite_length);
2415 std::vector<bool> at_boundary(keep_boundary ?
triangulation.n_vertices() :
2420 const bool is_parallel_shared =
2423 for (
const auto &cell :
triangulation.active_cell_iterators())
2424 if (is_parallel_shared || cell->is_locally_owned())
2428 for (
unsigned int i = 0; i < cell->n_lines(); ++i)
2431 line = cell->line(i);
2433 if (keep_boundary && line->at_boundary())
2435 at_boundary[line->vertex_index(0)] =
true;
2436 at_boundary[line->vertex_index(1)] =
true;
2439 minimal_length[line->vertex_index(0)] =
2441 minimal_length[line->vertex_index(0)]);
2442 minimal_length[line->vertex_index(1)] =
2444 minimal_length[line->vertex_index(1)]);
2450 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
2451 if (cell->at_boundary(vertex) ==
true)
2452 at_boundary[cell->vertex_index(vertex)] =
true;
2454 minimal_length[cell->vertex_index(0)] =
2456 minimal_length[cell->vertex_index(0)]);
2457 minimal_length[cell->vertex_index(1)] =
2459 minimal_length[cell->vertex_index(1)]);
2464 boost::random::mt19937 rng(seed);
2465 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
2469 if (
auto distributed_triangulation =
2473 const std::vector<bool> locally_owned_vertices =
2475 std::vector<bool> vertex_moved(
triangulation.n_vertices(),
false);
2478 for (
const auto &cell :
triangulation.active_cell_iterators())
2479 if (cell->is_locally_owned())
2481 for (
const unsigned int vertex_no : cell->vertex_indices())
2483 const unsigned global_vertex_no =
2484 cell->vertex_index(vertex_no);
2489 if ((keep_boundary && at_boundary[global_vertex_no]) ||
2490 vertex_moved[global_vertex_no] ||
2491 !locally_owned_vertices[global_vertex_no])
2496 for (
unsigned int d = 0;
d < spacedim; ++
d)
2497 shift_vector(
d) = uniform_distribution(rng);
2499 shift_vector *= factor * minimal_length[global_vertex_no] /
2500 std::sqrt(shift_vector.
square());
2503 cell->vertex(vertex_no) += shift_vector;
2504 vertex_moved[global_vertex_no] =
true;
2508 distributed_triangulation->communicate_locally_moved_vertices(
2509 locally_owned_vertices);
2519 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
2520 const std::vector<Point<spacedim>> &old_vertex_locations =
2523 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2527 if (keep_boundary && at_boundary[vertex])
2528 new_vertex_locations[vertex] = old_vertex_locations[vertex];
2533 for (
unsigned int d = 0;
d < spacedim; ++
d)
2534 shift_vector(
d) = uniform_distribution(rng);
2536 shift_vector *= factor * minimal_length[vertex] /
2537 std::sqrt(shift_vector.
square());
2540 new_vertex_locations[vertex] =
2541 old_vertex_locations[vertex] + shift_vector;
2546 for (
const auto &cell :
triangulation.active_cell_iterators())
2547 for (
const unsigned int vertex_no : cell->vertex_indices())
2548 cell->vertex(vertex_no) =
2549 new_vertex_locations[cell->vertex_index(vertex_no)];
2565 for (; cell != endc; ++cell)
2566 if (!cell->is_artificial())
2567 for (
const unsigned int face : cell->face_indices())
2568 if (cell->face(face)->has_children() &&
2569 !cell->face(face)->at_boundary())
2573 cell->face(face)->child(0)->vertex(1) =
2574 (cell->face(face)->vertex(0) +
2575 cell->face(face)->vertex(1)) /
2579 cell->face(face)->child(0)->vertex(1) =
2580 .5 * (cell->face(face)->vertex(0) +
2581 cell->face(face)->vertex(1));
2582 cell->face(face)->child(0)->vertex(2) =
2583 .5 * (cell->face(face)->vertex(0) +
2584 cell->face(face)->vertex(2));
2585 cell->face(face)->child(1)->vertex(3) =
2586 .5 * (cell->face(face)->vertex(1) +
2587 cell->face(face)->vertex(3));
2588 cell->face(face)->child(2)->vertex(3) =
2589 .5 * (cell->face(face)->vertex(2) +
2590 cell->face(face)->vertex(3));
2593 cell->face(face)->child(0)->vertex(3) =
2594 .25 * (cell->face(face)->vertex(0) +
2595 cell->face(face)->vertex(1) +
2596 cell->face(face)->vertex(2) +
2597 cell->face(face)->vertex(3));
2605 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2609 const
Point<spacedim> & p,
2610 const std::vector<
bool> &marked_vertices)
2619 marked_vertices.size() == 0,
2621 marked_vertices.size()));
2628 marked_vertices.size() == 0 ||
2629 std::equal(marked_vertices.begin(),
2630 marked_vertices.end(),
2632 [](
bool p,
bool q) { return !p || q; }),
2634 "marked_vertices should be a subset of used vertices in the triangulation "
2635 "but marked_vertices contains one or more vertices that are not used vertices!"));
2639 const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
2645 std::vector<bool>::const_iterator
first =
2646 std::find(vertices_to_use.begin(), vertices_to_use.end(),
true);
2651 unsigned int best_vertex = std::distance(vertices_to_use.begin(),
first);
2652 double best_dist = (p -
vertices[best_vertex]).norm_square();
2656 for (
unsigned int j = best_vertex + 1; j <
vertices.size(); ++j)
2657 if (vertices_to_use[j])
2659 const double dist = (p -
vertices[j]).norm_square();
2660 if (dist < best_dist)
2672 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2676 const MeshType<dim, spacedim> &mesh,
2677 const
Point<spacedim> & p,
2678 const std::vector<
bool> &marked_vertices)
2681 if (mapping.preserves_vertex_locations() ==
true)
2691 marked_vertices.size() == 0,
2693 marked_vertices.size()));
2701 marked_vertices.size() == 0 ||
2702 std::equal(marked_vertices.begin(),
2703 marked_vertices.end(),
2705 [](
bool p,
bool q) { return !p || q; }),
2707 "marked_vertices should be a subset of used vertices in the triangulation "
2708 "but marked_vertices contains one or more vertices that are not used vertices!"));
2711 if (marked_vertices.size() != 0)
2714 if (marked_vertices[it->first] ==
false)
2729 template <
int dim,
int spacedim>
2730 std::vector<std::vector<Tensor<1, spacedim>>>
2738 const unsigned int n_vertices = vertex_to_cells.size();
2743 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
2745 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2748 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
2749 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
2751 typename std::set<typename Triangulation<dim, spacedim>::
2752 active_cell_iterator>::iterator it =
2753 vertex_to_cells[vertex].begin();
2754 for (
unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
2756 vertex_to_cell_centers[vertex][cell] =
2757 (*it)->center() -
vertices[vertex];
2758 vertex_to_cell_centers[vertex][cell] /=
2759 vertex_to_cell_centers[vertex][cell].
norm();
2762 return vertex_to_cell_centers;
2768 template <
int spacedim>
2771 const unsigned int a,
2772 const unsigned int b,
2776 const double scalar_product_a = center_directions[a] * point_direction;
2777 const double scalar_product_b = center_directions[
b] * point_direction;
2782 return (scalar_product_a > scalar_product_b);
2786 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2790 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
2792 std::pair<typename ::internal::
2793 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
2798 const MeshType<dim, spacedim> &mesh,
2801 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
2804 &vertex_to_cell_centers,
2805 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
2806 const std::vector<bool> &marked_vertices,
2808 & used_vertices_rtree,
2809 const double tolerance,
2813 *relevant_cell_bounding_boxes_rtree)
2815 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2818 cell_and_position.first = mesh.end();
2822 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2824 cell_and_position_approx;
2826 if (relevant_cell_bounding_boxes_rtree !=
nullptr &&
2827 !relevant_cell_bounding_boxes_rtree->empty())
2833 for (
unsigned int d = 0;
d < spacedim; ++
d)
2835 p1[
d] = p1[
d] - tolerance;
2836 p2[
d] = p2[
d] + tolerance;
2841 if (relevant_cell_bounding_boxes_rtree->qbegin(
2842 boost::geometry::index::intersects(bb)) ==
2843 relevant_cell_bounding_boxes_rtree->qend())
2844 return cell_and_position;
2847 bool found_cell =
false;
2848 bool approx_cell =
false;
2850 unsigned int closest_vertex_index = 0;
2855 if (marked_vertices.size() && !used_vertices_rtree.empty())
2858 std::find(marked_vertices.begin(), marked_vertices.end(),
true);
2859 Assert(itr != marked_vertices.end(),
2861 closest_vertex_index = std::distance(marked_vertices.begin(), itr);
2865 auto current_cell = cell_hint;
2868 const auto cell_marked = [&mesh, &marked_vertices](
const auto &cell) {
2869 if (marked_vertices.size() == 0)
2872 if (cell != mesh.active_cell_iterators().end())
2873 for (
unsigned int i = 0; i < cell->n_vertices(); ++i)
2874 if (marked_vertices[cell->vertex_index(i)])
2881 const auto any_cell_marked = [&cell_marked](
const auto &cells) {
2882 return std::any_of(cells.begin(),
2884 [&cell_marked](
const auto &cell) {
2885 return cell_marked(cell);
2889 while (found_cell ==
false)
2894 cell_marked(cell_hint))
2896 const auto cell_vertices = mapping.
get_vertices(current_cell);
2897 const unsigned int closest_vertex =
2898 find_closest_vertex_of_cell<dim, spacedim>(current_cell,
2901 vertex_to_point = p - cell_vertices[closest_vertex];
2902 closest_vertex_index = current_cell->vertex_index(closest_vertex);
2911 #if defined(DEAL_II_WITH_BOOST_BUNDLED) || \
2912 !(defined(__clang_major__) && __clang_major__ >= 16) || \
2913 BOOST_VERSION >= 108100
2914 if (!used_vertices_rtree.empty())
2917 using ValueType = std::pair<Point<spacedim>,
unsigned int>;
2918 std::function<
bool(
const ValueType &)> marked;
2919 if (marked_vertices.size() == mesh.n_vertices())
2920 marked = [&marked_vertices](
const ValueType &value) ->
bool {
2921 return marked_vertices[value.second];
2924 marked = [](
const ValueType &) ->
bool {
return true; };
2926 std::vector<std::pair<Point<spacedim>,
unsigned int>> res;
2927 used_vertices_rtree.query(
2928 boost::geometry::index::nearest(p, 1) &&
2929 boost::geometry::index::satisfies(marked),
2930 std::back_inserter(res));
2935 ::
ExcMessage(
"There can not be multiple results"));
2938 if (any_cell_marked(vertex_to_cells[res[0].
second]))
2939 closest_vertex_index = res[0].second;
2945 mapping, mesh, p, marked_vertices);
2947 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
2953 Assert(any_cell_marked(vertex_to_cells[closest_vertex_index]),
2958 const double vertex_point_norm = vertex_to_point.
norm();
2959 if (vertex_point_norm > 0)
2960 vertex_to_point /= vertex_point_norm;
2962 const unsigned int n_neighbor_cells =
2963 vertex_to_cells[closest_vertex_index].size();
2966 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
2968 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
2969 neighbor_permutation[i] = i;
2971 auto comp = [&](
const unsigned int a,
const unsigned int b) ->
bool {
2972 return internal::compare_point_association<spacedim>(
2976 vertex_to_cell_centers[closest_vertex_index]);
2979 std::sort(neighbor_permutation.begin(),
2980 neighbor_permutation.end(),
2985 double best_distance = tolerance;
2989 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
2993 auto cell = vertex_to_cells[closest_vertex_index].begin();
2996 if (!(*cell)->is_artificial())
3003 cell_and_position.first = *cell;
3004 cell_and_position.second = p_unit;
3006 approx_cell =
false;
3015 if (dist < best_distance)
3017 best_distance = dist;
3018 cell_and_position_approx.first = *cell;
3019 cell_and_position_approx.second = p_unit;
3028 if (found_cell ==
true)
3029 return cell_and_position;
3030 else if (approx_cell ==
true)
3031 return cell_and_position_approx;
3047 mapping, mesh, p, marked_vertices, tolerance);
3049 current_cell =
typename MeshType<dim, spacedim>::active_cell_iterator();
3051 return cell_and_position;
3056 template <
int dim,
int spacedim>
3065 unsigned int closest_vertex = 0;
3067 for (
unsigned int v = 1; v < cell->
n_vertices(); ++v)
3070 if (vertex_distance < minimum_distance)
3073 minimum_distance = vertex_distance;
3076 return closest_vertex;
3083 namespace BoundingBoxPredicate
3085 template <
class MeshType>
3087 concepts::is_triangulation_or_dof_handler<MeshType>)
3091 cell_iterator &parent_cell,
3092 const std::function<
bool(
3093 const typename MeshType::
3094 active_cell_iterator &)>
3097 bool has_predicate =
3099 std::vector<typename MeshType::active_cell_iterator> active_cells;
3100 if (parent_cell->is_active())
3101 active_cells = {parent_cell};
3105 active_cells = get_active_child_cells<MeshType>(parent_cell);
3107 const unsigned int spacedim = MeshType::space_dimension;
3111 while (i < active_cells.size() && !predicate(active_cells[i]))
3115 if (active_cells.size() == 0 || i == active_cells.size())
3118 return std::make_tuple(bbox, has_predicate);
3125 for (; i < active_cells.size(); ++i)
3126 if (predicate(active_cells[i]))
3128 for (
unsigned int d = 0;
d < spacedim; ++
d)
3130 minp[
d] =
std::min(minp[
d], active_cells[i]->vertex(v)[
d]);
3131 maxp[
d] =
std::max(maxp[
d], active_cells[i]->vertex(v)[
d]);
3134 has_predicate =
true;
3136 return std::make_tuple(bbox, has_predicate);
3143 template <
class MeshType>
3147 const MeshType &mesh,
3148 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
3150 const unsigned int refinement_level,
3151 const bool allow_merge,
3152 const unsigned int max_boxes)
3159 refinement_level <= mesh.n_levels(),
3161 "Error: refinement level is higher then total levels in the triangulation!"));
3163 const unsigned int spacedim = MeshType::space_dimension;
3164 std::vector<BoundingBox<spacedim>> bounding_boxes;
3168 for (
unsigned int i = 0; i < refinement_level; ++i)
3169 for (
const typename MeshType::cell_iterator &cell :
3170 mesh.active_cell_iterators_on_level(i))
3172 bool has_predicate =
false;
3174 std::tie(bbox, has_predicate) =
3176 MeshType>(cell, predicate);
3178 bounding_boxes.push_back(bbox);
3182 for (
const typename MeshType::cell_iterator &cell :
3183 mesh.cell_iterators_on_level(refinement_level))
3185 bool has_predicate =
false;
3187 std::tie(bbox, has_predicate) =
3189 MeshType>(cell, predicate);
3191 bounding_boxes.push_back(bbox);
3196 return bounding_boxes;
3202 std::vector<unsigned int> merged_boxes_idx;
3203 bool found_neighbors =
true;
3208 while (found_neighbors)
3210 found_neighbors =
false;
3211 for (
unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
3213 if (std::find(merged_boxes_idx.begin(),
3214 merged_boxes_idx.end(),
3215 i) == merged_boxes_idx.end())
3216 for (
unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
3217 if (std::find(merged_boxes_idx.begin(),
3218 merged_boxes_idx.end(),
3219 j) == merged_boxes_idx.end() &&
3220 bounding_boxes[i].get_neighbor_type(
3221 bounding_boxes[j]) ==
3224 bounding_boxes[i].merge_with(bounding_boxes[j]);
3225 merged_boxes_idx.push_back(j);
3226 found_neighbors =
true;
3232 std::vector<BoundingBox<spacedim>> merged_b_boxes;
3233 for (
unsigned int i = 0; i < bounding_boxes.size(); ++i)
3234 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
3235 merged_boxes_idx.end())
3236 merged_b_boxes.push_back(bounding_boxes[i]);
3241 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
3243 std::vector<double> volumes;
3244 for (
unsigned int i = 0; i < merged_b_boxes.size(); ++i)
3245 volumes.push_back(merged_b_boxes[i].volume());
3247 while (merged_b_boxes.size() > max_boxes)
3249 unsigned int min_idx =
3250 std::min_element(volumes.begin(), volumes.end()) -
3252 volumes.erase(volumes.begin() + min_idx);
3254 bool not_removed =
true;
3255 for (
unsigned int i = 0;
3256 i < merged_b_boxes.size() && not_removed;
3261 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
3262 merged_b_boxes[min_idx]) ==
3264 merged_b_boxes[i].get_neighbor_type(
3265 merged_b_boxes[min_idx]) ==
3268 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
3269 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
3270 not_removed =
false;
3273 ExcMessage(
"Error: couldn't merge bounding boxes!"));
3276 Assert(merged_b_boxes.size() <= max_boxes,
3278 "Error: couldn't reach target number of bounding boxes!"));
3279 return merged_b_boxes;
3285 template <
int spacedim>
3287 std::tuple<std::vector<std::vector<unsigned int>>,
3288 std::map<unsigned int, unsigned int>,
3289 std::map<unsigned int, std::vector<unsigned int>>>
3297 unsigned int n_procs = global_bboxes.size();
3298 std::vector<std::vector<unsigned int>> point_owners(n_procs);
3299 std::map<unsigned int, unsigned int> map_owners_found;
3300 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3302 unsigned int n_points = points.size();
3303 for (
unsigned int pt = 0; pt < n_points; ++pt)
3306 std::vector<unsigned int> owners_found;
3308 for (
unsigned int rk = 0; rk < n_procs; ++rk)
3311 if (bbox.point_inside(points[pt]))
3313 point_owners[rk].emplace_back(pt);
3314 owners_found.emplace_back(rk);
3318 Assert(owners_found.size() > 0,
3319 ExcMessage(
"No owners found for the point " +
3321 if (owners_found.size() == 1)
3322 map_owners_found[pt] = owners_found[0];
3325 map_owners_guessed[pt] = owners_found;
3328 return std::make_tuple(std::move(point_owners),
3329 std::move(map_owners_found),
3330 std::move(map_owners_guessed));
3333 template <
int spacedim>
3335 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
3336 std::map<unsigned int, unsigned int>,
3337 std::map<unsigned int, std::vector<unsigned int>>>
3345 std::map<unsigned int, std::vector<unsigned int>> point_owners;
3346 std::map<unsigned int, unsigned int> map_owners_found;
3347 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3348 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> search_result;
3350 unsigned int n_points = points.size();
3351 for (
unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
3353 search_result.clear();
3356 covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
3357 std::back_inserter(search_result));
3360 std::set<unsigned int> owners_found;
3362 for (
const auto &rank_bbox : search_result)
3366 const bool pt_inserted = owners_found.insert(pt_n).second;
3368 point_owners[rank_bbox.second].emplace_back(pt_n);
3370 Assert(owners_found.size() > 0,
3371 ExcMessage(
"No owners found for the point " +
3373 if (owners_found.size() == 1)
3374 map_owners_found[pt_n] = *owners_found.begin();
3379 std::back_inserter(map_owners_guessed[pt_n]));
3382 return std::make_tuple(std::move(point_owners),
3383 std::move(map_owners_found),
3384 std::move(map_owners_guessed));
3388 template <
int dim,
int spacedim>
3390 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3394 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3399 for (; cell != endc; ++cell)
3400 for (
const unsigned int i : cell->vertex_indices())
3405 for (; cell != endc; ++cell)
3407 for (
unsigned int i : cell->face_indices())
3409 if ((cell->at_boundary(i) ==
false) &&
3410 (cell->neighbor(i)->is_active()))
3413 adjacent_cell = cell->neighbor(i);
3414 for (
unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
3423 for (
unsigned int i = 0; i < cell->
n_lines(); ++i)
3424 if (cell->line(i)->has_children())
3438 template <
int dim,
int spacedim>
3439 std::map<unsigned int, types::global_vertex_index>
3443 std::map<unsigned int, types::global_vertex_index>
3444 local_to_global_vertex_index;
3446 #ifndef DEAL_II_WITH_MPI
3453 ExcMessage(
"This function does not make any sense "
3454 "for parallel::distributed::Triangulation "
3455 "objects if you do not have MPI enabled."));
3459 using active_cell_iterator =
3461 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
3466 unsigned int max_cellid_size = 0;
3467 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
3469 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
3475 std::set<active_cell_iterator> missing_vert_cells;
3476 std::set<unsigned int> used_vertex_index;
3477 for (
const auto &cell :
triangulation.active_cell_iterators())
3479 if (cell->is_locally_owned())
3481 for (
const unsigned int i : cell->vertex_indices())
3484 for (
const auto &adjacent_cell :
3485 vertex_to_cell[cell->vertex_index(i)])
3486 lowest_subdomain_id =
std::min(lowest_subdomain_id,
3487 adjacent_cell->subdomain_id());
3490 if (lowest_subdomain_id == cell->subdomain_id())
3494 if (used_vertex_index.find(cell->vertex_index(i)) ==
3495 used_vertex_index.end())
3498 local_to_global_vertex_index[cell->vertex_index(i)] =
3503 for (
const auto &adjacent_cell :
3504 vertex_to_cell[cell->vertex_index(i)])
3505 if (adjacent_cell->subdomain_id() !=
3506 cell->subdomain_id())
3510 tmp(adjacent_cell->subdomain_id(),
3511 cell->vertex_index(i));
3512 if (vertices_added.find(tmp) ==
3513 vertices_added.end())
3515 vertices_to_send[adjacent_cell
3518 cell->vertex_index(i),
3519 cell->id().to_string());
3520 if (cell->id().to_string().size() >
3523 cell->id().to_string().size();
3524 vertices_added.insert(tmp);
3527 used_vertex_index.insert(cell->vertex_index(i));
3534 vertices_to_recv[lowest_subdomain_id].insert(
3535 cell->vertex_index(i));
3536 missing_vert_cells.insert(cell);
3543 if (cell->is_ghost())
3545 for (
unsigned int i : cell->face_indices())
3547 if (cell->at_boundary(i) ==
false)
3549 if (cell->neighbor(i)->is_active())
3552 spacedim>::active_cell_iterator
3553 adjacent_cell = cell->neighbor(i);
3554 if ((adjacent_cell->is_locally_owned()))
3557 adjacent_cell->subdomain_id();
3558 if (cell->subdomain_id() < adj_subdomain_id)
3559 for (
unsigned int j = 0;
3560 j < cell->face(i)->n_vertices();
3563 vertices_to_recv[cell->subdomain_id()].insert(
3564 cell->face(i)->vertex_index(j));
3565 missing_vert_cells.insert(cell);
3581 int ierr = MPI_Exscan(&next_index,
3589 for (
auto &global_index_it : local_to_global_vertex_index)
3590 global_index_it.second +=
shift;
3604 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
3605 vertices_to_send.size());
3606 std::vector<MPI_Request> first_requests(vertices_to_send.size());
3610 std::string>>>::iterator
3611 vert_to_send_it = vertices_to_send.begin(),
3612 vert_to_send_end = vertices_to_send.end();
3613 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3614 ++vert_to_send_it, ++i)
3616 int destination = vert_to_send_it->first;
3617 const unsigned int n_vertices = vert_to_send_it->second.size();
3618 const int buffer_size = 2 * n_vertices;
3619 vertices_send_buffers[i].resize(buffer_size);
3622 for (
unsigned int j = 0; j < n_vertices; ++j)
3624 vertices_send_buffers[i][2 * j] =
3625 std::get<0>(vert_to_send_it->second[j]);
3626 vertices_send_buffers[i][2 * j + 1] =
3627 local_to_global_vertex_index[std::get<1>(
3628 vert_to_send_it->second[j])];
3632 ierr = MPI_Isend(vertices_send_buffers[i].data(),
3638 &first_requests[i]);
3643 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
3644 vertices_to_recv.size());
3645 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
3646 vert_to_recv_it = vertices_to_recv.begin(),
3647 vert_to_recv_end = vertices_to_recv.end();
3648 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3649 ++vert_to_recv_it, ++i)
3651 int source = vert_to_recv_it->first;
3652 const unsigned int n_vertices = vert_to_recv_it->second.size();
3653 const int buffer_size = 2 * n_vertices;
3654 vertices_recv_buffers[i].resize(buffer_size);
3657 ierr = MPI_Recv(vertices_recv_buffers[i].data(),
3669 std::vector<std::vector<char>> cellids_send_buffers(
3670 vertices_to_send.size());
3671 std::vector<MPI_Request> second_requests(vertices_to_send.size());
3672 vert_to_send_it = vertices_to_send.begin();
3673 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3674 ++vert_to_send_it, ++i)
3676 int destination = vert_to_send_it->first;
3677 const unsigned int n_vertices = vert_to_send_it->second.size();
3678 const int buffer_size = max_cellid_size * n_vertices;
3679 cellids_send_buffers[i].resize(buffer_size);
3682 unsigned int pos = 0;
3683 for (
unsigned int j = 0; j < n_vertices; ++j)
3685 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
3686 for (
unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
3688 if (k < cell_id.size())
3689 cellids_send_buffers[i][pos] = cell_id[k];
3693 cellids_send_buffers[i][pos] =
'-';
3698 ierr = MPI_Isend(cellids_send_buffers[i].data(),
3704 &second_requests[i]);
3709 std::vector<std::vector<char>> cellids_recv_buffers(
3710 vertices_to_recv.size());
3711 vert_to_recv_it = vertices_to_recv.begin();
3712 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3713 ++vert_to_recv_it, ++i)
3715 int source = vert_to_recv_it->first;
3716 const unsigned int n_vertices = vert_to_recv_it->second.size();
3717 const int buffer_size = max_cellid_size * n_vertices;
3718 cellids_recv_buffers[i].resize(buffer_size);
3721 ierr = MPI_Recv(cellids_recv_buffers[i].data(),
3733 vert_to_recv_it = vertices_to_recv.begin();
3734 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3735 ++i, ++vert_to_recv_it)
3737 for (
unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
3739 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
3741 vertices_recv_buffers[i][2 * j + 1];
3742 const std::string cellid_recv(
3743 &cellids_recv_buffers[i][max_cellid_size * j],
3744 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
3746 typename std::set<active_cell_iterator>::iterator
3747 cell_set_it = missing_vert_cells.begin(),
3748 end_cell_set = missing_vert_cells.end();
3749 for (; (found ==
false) && (cell_set_it != end_cell_set);
3752 typename std::set<active_cell_iterator>::iterator
3754 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
3756 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
3757 for (; candidate_cell != end_cell; ++candidate_cell)
3759 std::string current_cellid =
3760 (*candidate_cell)->id().to_string();
3761 current_cellid.resize(max_cellid_size,
'-');
3762 if (current_cellid.compare(cellid_recv) == 0)
3764 local_to_global_vertex_index
3765 [(*candidate_cell)->vertex_index(local_pos_recv)] =
3777 return local_to_global_vertex_index;
3782 template <
int dim,
int spacedim>
3798 for (
const auto &cell :
triangulation.active_cell_iterators())
3800 const unsigned int index = cell->active_cell_index();
3801 cell_connectivity.
add(index, index);
3802 for (
auto f : cell->face_indices())
3803 if ((cell->at_boundary(f) ==
false) &&
3804 (cell->neighbor(f)->has_children() ==
false))
3806 const unsigned int other_index =
3807 cell->neighbor(f)->active_cell_index();
3808 cell_connectivity.
add(index, other_index);
3809 cell_connectivity.
add(other_index, index);
3816 template <
int dim,
int spacedim>
3822 std::vector<std::vector<unsigned int>> vertex_to_cell(
3824 for (
const auto &cell :
triangulation.active_cell_iterators())
3826 for (
const unsigned int v : cell->vertex_indices())
3827 vertex_to_cell[cell->vertex_index(v)].push_back(
3828 cell->active_cell_index());
3833 for (
const auto &cell :
triangulation.active_cell_iterators())
3835 for (
const unsigned int v : cell->vertex_indices())
3836 for (
unsigned int n = 0;
3837 n < vertex_to_cell[cell->vertex_index(v)].size();
3839 cell_connectivity.
add(cell->active_cell_index(),
3840 vertex_to_cell[cell->vertex_index(v)][n]);
3845 template <
int dim,
int spacedim>
3849 const unsigned int level,
3852 std::vector<std::vector<unsigned int>> vertex_to_cell(
3859 for (
const unsigned int v : cell->vertex_indices())
3860 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
3870 for (
const unsigned int v : cell->vertex_indices())
3871 for (
unsigned int n = 0;
3872 n < vertex_to_cell[cell->vertex_index(v)].size();
3874 cell_connectivity.
add(cell->index(),
3875 vertex_to_cell[cell->vertex_index(v)][n]);
3881 template <
int dim,
int spacedim>
3889 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3890 "are already partitioned implicitly and can not be "
3891 "partitioned again explicitly."));
3893 std::vector<unsigned int> cell_weights;
3903 for (
const auto &cell :
triangulation.active_cell_iterators())
3904 if (cell->is_locally_owned())
3905 cell_weights[cell->active_cell_index()] =
3915 if (
const auto shared_tria =
3919 shared_tria->get_communicator(),
3923 Assert(std::accumulate(cell_weights.begin(),
3925 std::uint64_t(0)) > 0,
3926 ExcMessage(
"The global sum of weights over all active cells "
3927 "is zero. Please verify how you generate weights."));
3939 template <
int dim,
int spacedim>
3942 const std::vector<unsigned int> &cell_weights,
3948 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3949 "are already partitioned implicitly and can not be "
3950 "partitioned again explicitly."));
3954 if (n_partitions == 1)
3956 for (
const auto &cell :
triangulation.active_cell_iterators())
3957 cell->set_subdomain_id(0);
3971 sp_cell_connectivity.
copy_from(cell_connectivity);
3974 sp_cell_connectivity,
3981 template <
int dim,
int spacedim>
3990 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3991 "are already partitioned implicitly and can not be "
3992 "partitioned again explicitly."));
3994 std::vector<unsigned int> cell_weights;
4004 for (
const auto &cell :
triangulation.active_cell_iterators() |
4006 cell_weights[cell->active_cell_index()] =
4016 if (
const auto shared_tria =
4020 shared_tria->get_communicator(),
4024 Assert(std::accumulate(cell_weights.begin(),
4026 std::uint64_t(0)) > 0,
4027 ExcMessage(
"The global sum of weights over all active cells "
4028 "is zero. Please verify how you generate weights."));
4034 cell_connection_graph,
4041 template <
int dim,
int spacedim>
4044 const std::vector<unsigned int> &cell_weights,
4051 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4052 "are already partitioned implicitly and can not be "
4053 "partitioned again explicitly."));
4056 ExcMessage(
"Connectivity graph has wrong size"));
4058 ExcMessage(
"Connectivity graph has wrong size"));
4064 if (n_partitions == 1)
4066 for (
const auto &cell :
triangulation.active_cell_iterators())
4067 cell->set_subdomain_id(0);
4075 std::vector<unsigned int> partition_indices(
triangulation.n_active_cells());
4083 for (
const auto &cell :
triangulation.active_cell_iterators())
4084 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
4096 unsigned int & current_proc_idx,
4097 unsigned int & current_cell_idx,
4099 const unsigned int n_partitions)
4101 if (cell->is_active())
4103 while (current_cell_idx >=
4105 (current_proc_idx + 1) / n_partitions))
4107 cell->set_subdomain_id(current_proc_idx);
4112 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4122 template <
int dim,
int spacedim>
4126 const bool group_siblings)
4130 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4131 "are already partitioned implicitly and can not be "
4132 "partitioned again explicitly."));
4140 if (n_partitions == 1)
4142 for (
const auto &cell :
triangulation.active_cell_iterators())
4143 cell->set_subdomain_id(0);
4149 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
4150 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
4156 coarse_cell_to_p4est_tree_permutation.resize(
triangulation.n_cells(0));
4158 coarse_cell_to_p4est_tree_permutation);
4160 p4est_tree_to_coarse_cell_permutation =
4163 unsigned int current_proc_idx = 0;
4164 unsigned int current_cell_idx = 0;
4169 for (
unsigned int idx = 0; idx <
triangulation.n_cells(0); ++idx)
4171 const unsigned int coarse_cell_idx =
4172 p4est_tree_to_coarse_cell_permutation[idx];
4195 for (; cell != endc; ++cell)
4197 if (cell->is_active())
4199 bool all_children_active =
true;
4200 std::map<unsigned int, unsigned int> map_cpu_n_cells;
4201 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4202 if (!cell->child(n)->is_active())
4204 all_children_active =
false;
4208 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
4210 if (!all_children_active)
4213 unsigned int new_owner = cell->child(0)->subdomain_id();
4214 for (std::map<unsigned int, unsigned int>::iterator it =
4215 map_cpu_n_cells.begin();
4216 it != map_cpu_n_cells.end();
4218 if (it->second > map_cpu_n_cells[new_owner])
4219 new_owner = it->first;
4221 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4222 cell->child(n)->set_subdomain_id(new_owner);
4228 template <
int dim,
int spacedim>
4233 for (
int lvl = n_levels - 1; lvl >= 0; --lvl)
4235 for (
const auto &cell :
triangulation.cell_iterators_on_level(lvl))
4237 if (cell->is_active())
4238 cell->set_level_subdomain_id(cell->subdomain_id());
4241 Assert(cell->child(0)->level_subdomain_id() !=
4244 cell->set_level_subdomain_id(
4245 cell->child(0)->level_subdomain_id());
4257 template <
int dim,
int spacedim>
4262 const std::vector<CellId> & cell_ids,
4263 std::vector<types::subdomain_id> &subdomain_ids)
4265 #ifndef DEAL_II_WITH_P4EST
4268 (void)subdomain_ids;
4272 "You are attempting to use a functionality that is only available "
4273 "if deal.II was configured to use p4est, but cmake did not find a "
4274 "valid p4est library."));
4279 for (
const auto &cell_id : cell_ids)
4282 typename ::internal::p4est::types<dim>::quadrant p4est_cell,
4285 ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
4286 for (
const auto &child_index : cell_id.get_child_indices())
4288 ::internal::p4est::init_quadrant_children<dim>(
4289 p4est_cell, p4est_children);
4291 p4est_children[
static_cast<unsigned int>(child_index)];
4297 const_cast<typename ::internal::p4est::types<dim>::forest
4299 cell_id.get_coarse_cell_id(),
4306 subdomain_ids.push_back(owner);
4313 template <
int spacedim>
4317 const std::vector<CellId> &,
4318 std::vector<types::subdomain_id> &)
4327 template <
int dim,
int spacedim>
4328 std::vector<types::subdomain_id>
4330 const std::vector<CellId> & cell_ids)
4332 std::vector<types::subdomain_id> subdomain_ids;
4333 subdomain_ids.reserve(cell_ids.size());
4342 *parallel_tria =
dynamic_cast<
4356 const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
4357 shared_tria->get_true_subdomain_ids_of_cells();
4359 for (
const auto &cell_id : cell_ids)
4361 const unsigned int active_cell_index =
4362 shared_tria->create_cell_iterator(cell_id)->active_cell_index();
4363 subdomain_ids.push_back(
4364 true_subdomain_ids_of_cells[active_cell_index]);
4371 for (
const auto &cell_id : cell_ids)
4373 subdomain_ids.push_back(
4374 triangulation.create_cell_iterator(cell_id)->subdomain_id());
4378 return subdomain_ids;
4383 template <
int dim,
int spacedim>
4386 std::vector<types::subdomain_id> & subdomain)
4391 for (
const auto &cell :
triangulation.active_cell_iterators())
4392 subdomain[cell->active_cell_index()] = cell->subdomain_id();
4397 template <
int dim,
int spacedim>
4403 unsigned int count = 0;
4404 for (
const auto &cell :
triangulation.active_cell_iterators())
4405 if (cell->subdomain_id() == subdomain)
4413 template <
int dim,
int spacedim>
4418 std::vector<bool> locally_owned_vertices =
4425 if (
const auto *tr =
dynamic_cast<
4428 for (
const auto &cell :
triangulation.active_cell_iterators())
4429 if (cell->is_artificial() ||
4430 (cell->is_ghost() &&
4431 (cell->subdomain_id() < tr->locally_owned_subdomain())))
4432 for (
const unsigned int v : cell->vertex_indices())
4433 locally_owned_vertices[cell->vertex_index(v)] =
false;
4435 return locally_owned_vertices;
4440 template <
int dim,
int spacedim>
4446 for (
const auto &cell :
triangulation.active_cell_iterators())
4447 if (!cell->is_artificial())
4448 min_diameter =
std::min(min_diameter, cell->diameter(mapping));
4450 double global_min_diameter = 0;
4452 #ifdef DEAL_II_WITH_MPI
4456 global_min_diameter =
4460 global_min_diameter = min_diameter;
4462 return global_min_diameter;
4467 template <
int dim,
int spacedim>
4472 double max_diameter = 0.;
4473 for (
const auto &cell :
triangulation.active_cell_iterators())
4474 if (!cell->is_artificial())
4475 max_diameter =
std::max(max_diameter, cell->diameter(mapping));
4477 double global_max_diameter = 0;
4479 #ifdef DEAL_II_WITH_MPI
4483 global_max_diameter =
4487 global_max_diameter = max_diameter;
4489 return global_max_diameter;
4496 namespace FixUpDistortedChildCells
4519 template <
typename Iterator,
int spacedim>
4524 const unsigned int structdim =
4525 Iterator::AccessorType::structure_dimension;
4526 Assert(spacedim == Iterator::AccessorType::dimension,
4532 Assert(object->refinement_case() ==
4540 Tensor<spacedim - structdim, spacedim>
4543 for (
const unsigned int i : object->vertex_indices())
4544 parent_vertices[i] =
object->vertex(i);
4547 parent_vertices, parent_alternating_forms);
4549 const Tensor<spacedim - structdim, spacedim>
4550 average_parent_alternating_form =
4551 std::accumulate(parent_alternating_forms,
4552 parent_alternating_forms +
4566 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4570 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4571 for (
const unsigned int i : object->child(c)->vertex_indices())
4572 child_vertices[c][i] =
object->child(c)->vertex(i);
4580 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4582 1] = object_mid_point;
4584 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4586 child_vertices[c], child_alternating_forms[c]);
4598 double objective = 0;
4599 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4600 for (
const unsigned int i : object->child(c)->vertex_indices())
4602 (child_alternating_forms[c][i] -
4603 average_parent_alternating_form / std::pow(2., 1. * structdim))
4615 template <
typename Iterator>
4618 const unsigned int f,
4619 std::integral_constant<int, 1>)
4621 return object->vertex(f);
4631 template <
typename Iterator>
4634 const unsigned int f,
4635 std::integral_constant<int, 2>)
4637 return object->line(f)->center();
4647 template <
typename Iterator>
4650 const unsigned int f,
4651 std::integral_constant<int, 3>)
4653 return object->face(f)->center();
4680 template <
typename Iterator>
4684 const unsigned int structdim =
4685 Iterator::AccessorType::structure_dimension;
4687 double diameter =
object->diameter();
4688 for (
const unsigned int f : object->face_indices())
4689 for (
unsigned int e = f + 1;
e <
object->n_faces(); ++
e)
4694 std::integral_constant<int, structdim>())
4696 object,
e, std::integral_constant<int, structdim>())));
4707 template <
typename Iterator>
4711 const unsigned int structdim =
4712 Iterator::AccessorType::structure_dimension;
4713 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
4720 Assert(object->refinement_case() ==
4730 unsigned int iteration = 0;
4735 double initial_delta = 0;
4742 const double step_length =
diameter / 4 / (iteration + 1);
4747 for (
unsigned int d = 0;
d < spacedim; ++
d)
4749 const double eps = step_length / 10;
4763 if (gradient.norm() == 0)
4772 std::min(2 * current_value / (gradient * gradient),
4773 step_length / gradient.norm()) *
4778 const double previous_value = current_value;
4782 initial_delta = (previous_value - current_value);
4785 if ((iteration >= 1) &&
4786 ((previous_value - current_value < 0) ||
4787 (
std::fabs(previous_value - current_value) <
4788 0.001 * initial_delta)))
4793 while (iteration < 20);
4809 double old_min_product, new_min_product;
4814 parent_vertices[i] = object->vertex(i);
4816 Tensor<spacedim - structdim, spacedim>
4819 parent_vertices, parent_alternating_forms);
4825 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4826 for (
const unsigned int i : object->child(c)->vertex_indices())
4827 child_vertices[c][i] =
object->child(c)->vertex(i);
4829 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4833 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4835 child_vertices[c], child_alternating_forms[c]);
4838 child_alternating_forms[0][0] * parent_alternating_forms[0];
4839 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4840 for (
const unsigned int i : object->child(c)->vertex_indices())
4841 for (
const unsigned int j : object->vertex_indices())
4842 old_min_product = std::min<double>(old_min_product,
4843 child_alternating_forms[c][i] *
4844 parent_alternating_forms[j]);
4852 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4854 1] = object_mid_point;
4856 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4858 child_vertices[c], child_alternating_forms[c]);
4861 child_alternating_forms[0][0] * parent_alternating_forms[0];
4862 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4863 for (
const unsigned int i : object->child(c)->vertex_indices())
4864 for (
const unsigned int j : object->vertex_indices())
4865 new_min_product = std::min<double>(new_min_product,
4866 child_alternating_forms[c][i] *
4867 parent_alternating_forms[j]);
4875 if (new_min_product >= old_min_product)
4876 object->child(0)->vertex(
4883 return (
std::max(new_min_product, old_min_product) > 0);
4889 template <
int dim,
int spacedim>
4892 const typename ::Triangulation<dim, spacedim>::cell_iterator
4894 std::integral_constant<int, dim>,
4895 std::integral_constant<int, spacedim>)
4905 for (
auto f : cell->face_indices())
4908 Assert(cell->face(f)->refinement_case() ==
4912 bool subface_is_more_refined =
false;
4913 for (
unsigned int g = 0;
4914 g < GeometryInfo<dim>::max_children_per_face;
4916 if (cell->face(f)->child(g)->has_children())
4918 subface_is_more_refined =
true;
4922 if (subface_is_more_refined ==
true)
4933 template <
int dim,
int spacedim>
4941 dim != 1 && spacedim != 1,
4942 "This function is only valid when dim != 1 or spacedim != 1.");
4946 for (
typename std::list<
4948 cell_ptr = distorted_cells.distorted_cells.
begin();
4949 cell_ptr != distorted_cells.distorted_cells.
end();
4955 Assert(!cell->is_active(),
4957 "This function is only valid for a list of cells that "
4958 "have children (i.e., no cell in the list may be active)."));
4962 std::integral_constant<int, dim>(),
4963 std::integral_constant<int, spacedim>());
4967 unfixable_subset.distorted_cells.push_back(cell);
4970 return unfixable_subset;
4975 template <
int dim,
int spacedim>
4978 const bool reset_boundary_ids)
4981 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
4982 auto m_it = dst_manifold_ids.begin();
4983 for (
const auto b : src_boundary_ids)
4988 const std::vector<types::boundary_id> reset_boundary_id =
4989 reset_boundary_ids ?
4990 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
5000 template <
int dim,
int spacedim>
5003 const std::vector<types::boundary_id> &src_boundary_ids,
5004 const std::vector<types::manifold_id> &dst_manifold_ids,
5006 const std::vector<types::boundary_id> &reset_boundary_ids_)
5009 const auto reset_boundary_ids =
5010 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
5019 for (
auto f : cell->face_indices())
5020 if (cell->face(f)->at_boundary())
5021 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
5023 const auto bid = cell->face(f)->line(
e)->boundary_id();
5024 const unsigned int ind = std::find(src_boundary_ids.begin(),
5025 src_boundary_ids.end(),
5027 src_boundary_ids.begin();
5028 if (ind < src_boundary_ids.size())
5029 cell->face(f)->line(
e)->set_manifold_id(
5030 dst_manifold_ids[ind]);
5035 for (
auto f : cell->face_indices())
5036 if (cell->face(f)->at_boundary())
5038 const auto bid = cell->face(f)->boundary_id();
5039 const unsigned int ind =
5040 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
5041 src_boundary_ids.begin();
5043 if (ind < src_boundary_ids.size())
5046 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
5048 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
5052 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
5054 const auto bid = cell->face(f)->line(
e)->boundary_id();
5055 const unsigned int ind = std::find(src_boundary_ids.begin(),
5056 src_boundary_ids.end(),
5058 src_boundary_ids.begin();
5059 if (ind < src_boundary_ids.size())
5060 cell->face(f)->line(
e)->set_boundary_id(
5061 reset_boundary_ids[ind]);
5067 template <
int dim,
int spacedim>
5070 const bool compute_face_ids)
5076 for (; cell != endc; ++cell)
5078 cell->set_manifold_id(cell->material_id());
5079 if (compute_face_ids ==
true)
5081 for (
auto f : cell->face_indices())
5083 if (cell->at_boundary(f) ==
false)
5084 cell->face(f)->set_manifold_id(
5086 cell->neighbor(f)->material_id()));
5088 cell->face(f)->set_manifold_id(cell->material_id());
5095 template <
int dim,
int spacedim>
5100 const std::set<types::manifold_id> &)> &disambiguation_function,
5101 bool overwrite_only_flat_manifold_ids)
5106 const unsigned int n_subobjects =
5110 std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
5111 std::vector<unsigned int> backup;
5115 unsigned next_index = 1;
5119 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
5121 if (cell->line(
l)->user_index() == 0)
5124 manifold_ids[next_index].insert(cell->manifold_id());
5125 cell->line(
l)->set_user_index(next_index++);
5128 manifold_ids[cell->line(
l)->user_index()].insert(
5129 cell->manifold_id());
5132 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
5134 if (cell->quad(
l)->user_index() == 0)
5137 manifold_ids[next_index].insert(cell->manifold_id());
5138 cell->quad(
l)->set_user_index(next_index++);
5141 manifold_ids[cell->quad(
l)->user_index()].insert(
5142 cell->manifold_id());
5148 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
5150 const auto id = cell->line(
l)->user_index();
5154 if (cell->line(
l)->manifold_id() ==
5156 overwrite_only_flat_manifold_ids ==
false)
5157 cell->line(
l)->set_manifold_id(
5158 disambiguation_function(manifold_ids[
id]));
5159 cell->line(
l)->set_user_index(0);
5163 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
5165 const auto id = cell->quad(
l)->user_index();
5169 if (cell->quad(
l)->manifold_id() ==
5171 overwrite_only_flat_manifold_ids ==
false)
5172 cell->quad(
l)->set_manifold_id(
5173 disambiguation_function(manifold_ids[
id]));
5174 cell->quad(
l)->set_user_index(0);
5183 template <
int dim,
int spacedim>
5184 std::pair<unsigned int, double>
5188 double max_ratio = 1;
5189 unsigned int index = 0;
5191 for (
unsigned int i = 0; i < dim; ++i)
5192 for (
unsigned int j = i + 1; j < dim; ++j)
5194 unsigned int ax = i % dim;
5195 unsigned int next_ax = j % dim;
5198 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
5200 if (ratio > max_ratio)
5205 else if (1.0 / ratio > max_ratio)
5207 max_ratio = 1.0 / ratio;
5211 return std::make_pair(index, max_ratio);
5215 template <
int dim,
int spacedim>
5218 const bool isotropic,
5219 const unsigned int max_iterations)
5221 unsigned int iter = 0;
5222 bool continue_refinement =
true;
5224 while (continue_refinement && (iter < max_iterations))
5228 continue_refinement =
false;
5231 for (
const unsigned int j : cell->face_indices())
5232 if (cell->at_boundary(j) ==
false &&
5233 cell->neighbor(j)->has_children())
5237 cell->set_refine_flag();
5238 continue_refinement =
true;
5241 continue_refinement |= cell->flag_for_face_refinement(j);
5248 template <
int dim,
int spacedim>
5251 const double max_ratio,
5252 const unsigned int max_iterations)
5254 unsigned int iter = 0;
5255 bool continue_refinement =
true;
5257 while (continue_refinement && (iter < max_iterations))
5260 continue_refinement =
false;
5263 std::pair<unsigned int, double> info =
5264 GridTools::get_longest_direction<dim, spacedim>(cell);
5265 if (info.second > max_ratio)
5267 cell->set_refine_flag(
5269 continue_refinement =
true;
5277 template <
int dim,
int spacedim>
5280 const double limit_angle_fraction)
5288 "have hanging nodes."));
5292 bool has_cells_with_more_than_dim_faces_on_boundary =
true;
5293 bool has_cells_with_dim_faces_on_boundary =
false;
5295 unsigned int refinement_cycles = 0;
5297 while (has_cells_with_more_than_dim_faces_on_boundary)
5299 has_cells_with_more_than_dim_faces_on_boundary =
false;
5303 unsigned int boundary_face_counter = 0;
5304 for (
auto f : cell->face_indices())
5305 if (cell->face(f)->at_boundary())
5306 boundary_face_counter++;
5307 if (boundary_face_counter > dim)
5309 has_cells_with_more_than_dim_faces_on_boundary =
true;
5312 else if (boundary_face_counter == dim)
5313 has_cells_with_dim_faces_on_boundary =
true;
5315 if (has_cells_with_more_than_dim_faces_on_boundary)
5318 refinement_cycles++;
5322 if (has_cells_with_dim_faces_on_boundary)
5325 refinement_cycles++;
5329 while (refinement_cycles > 0)
5332 cell->set_coarsen_flag();
5334 refinement_cycles--;
5344 std::vector<CellData<dim>> cells_to_add;
5348 const unsigned int v0 = 0,
v1 = 1, v2 = (dim > 1 ? 2 : 0),
5349 v3 = (dim > 1 ? 3 : 0);
5353 double angle_fraction = 0;
5359 p0[spacedim > 1 ? 1 : 0] = 1;
5363 if (cell->face(
v0)->at_boundary() && cell->face(v3)->at_boundary())
5365 p0 = cell->vertex(
v0) - cell->vertex(v2);
5366 p1 = cell->vertex(v3) - cell->vertex(v2);
5367 vertex_at_corner = v2;
5369 else if (cell->face(v3)->at_boundary() &&
5370 cell->face(
v1)->at_boundary())
5372 p0 = cell->vertex(v2) - cell->vertex(v3);
5373 p1 = cell->vertex(
v1) - cell->vertex(v3);
5374 vertex_at_corner = v3;
5376 else if (cell->face(1)->at_boundary() &&
5377 cell->face(2)->at_boundary())
5379 p0 = cell->vertex(
v0) - cell->vertex(
v1);
5380 p1 = cell->vertex(v3) - cell->vertex(
v1);
5381 vertex_at_corner =
v1;
5383 else if (cell->face(2)->at_boundary() &&
5384 cell->face(0)->at_boundary())
5386 p0 = cell->vertex(v2) - cell->vertex(
v0);
5387 p1 = cell->vertex(
v1) - cell->vertex(
v0);
5388 vertex_at_corner =
v0;
5399 if (angle_fraction > limit_angle_fraction)
5401 auto flags_removal = [&](
unsigned int f1,
5404 unsigned int n2) ->
void {
5405 cells_to_remove[cell->active_cell_index()] =
true;
5406 cells_to_remove[cell->neighbor(n1)->active_cell_index()] =
true;
5407 cells_to_remove[cell->neighbor(n2)->active_cell_index()] =
true;
5409 faces_to_remove[cell->face(f1)->index()] =
true;
5410 faces_to_remove[cell->face(f2)->index()] =
true;
5412 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] =
true;
5413 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] =
true;
5416 auto cell_creation = [&](
const unsigned int vv0,
5417 const unsigned int vv1,
5418 const unsigned int f0,
5419 const unsigned int f1,
5421 const unsigned int n0,
5422 const unsigned int v0n0,
5423 const unsigned int v1n0,
5425 const unsigned int n1,
5426 const unsigned int v0n1,
5427 const unsigned int v1n1) {
5433 c1.
vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
5434 c1.
vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
5440 c2.
vertices[
v1] = cell->neighbor(n1)->vertex_index(v0n1);
5441 c2.
vertices[v2] = cell->vertex_index(vv1);
5442 c2.
vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
5447 l1.
vertices[0] = cell->vertex_index(vv0);
5448 l1.
vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
5454 l2.
vertices[0] = cell->vertex_index(vv0);
5455 l2.
vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
5461 cells_to_add.push_back(c1);
5462 cells_to_add.push_back(c2);
5467 switch (vertex_at_corner)
5470 flags_removal(0, 2, 3, 1);
5471 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
5474 flags_removal(1, 2, 3, 0);
5475 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
5478 flags_removal(3, 0, 1, 2);
5479 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
5482 flags_removal(3, 1, 0, 2);
5483 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
5496 if (cells_to_add.size() == 0)
5498 while (refinement_cycles > 0)
5501 cell->set_coarsen_flag();
5503 refinement_cycles--;
5511 if (cells_to_remove[cell->active_cell_index()] ==
false)
5514 for (
const unsigned int v : cell->vertex_indices())
5515 c.
vertices[v] = cell->vertex_index(v);
5518 cells_to_add.push_back(c);
5526 for (; face != endf; ++face)
5527 if ((face->at_boundary() ||
5529 faces_to_remove[face->index()] ==
false)
5531 for (
unsigned int l = 0;
l < face->
n_lines(); ++
l)
5536 for (
const unsigned int v : face->vertex_indices())
5537 line.
vertices[v] = face->vertex_index(v);
5543 for (
const unsigned int v : face->line(
l)->vertex_indices())
5544 line.
vertices[v] = face->line(
l)->vertex_index(v);
5553 for (
const unsigned int v : face->vertex_indices())
5554 quad.
vertices[v] = face->vertex_index(v);
5562 subcelldata_to_add);
5567 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
5586 template <
int dim,
int spacedim>
5589 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5590 std::vector<std::vector<Point<dim>>>,
5591 std::vector<std::vector<unsigned int>>>
5603 auto &cells = std::get<0>(cqmp);
5604 auto &qpoints = std::get<1>(cqmp);
5605 auto &maps = std::get<2>(cqmp);
5607 return std::make_tuple(std::move(cells),
5614 template <
int dim,
int spacedim>
5617 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5618 std::vector<std::vector<Point<dim>>>,
5619 std::vector<std::vector<unsigned int>>,
5620 std::vector<unsigned int>>
5630 Assert((dim == spacedim),
5631 ExcMessage(
"Only implemented for dim==spacedim."));
5640 const unsigned int np = points.size();
5642 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5644 std::vector<std::vector<Point<dim>>> qpoints_out;
5645 std::vector<std::vector<unsigned int>> maps_out;
5646 std::vector<unsigned int> missing_points_out;
5650 return std::make_tuple(std::move(cells_out),
5651 std::move(qpoints_out),
5652 std::move(maps_out),
5653 std::move(missing_points_out));
5661 std::vector<std::pair<Point<spacedim>,
unsigned int>> points_and_ids(np);
5662 for (
unsigned int i = 0; i < np; ++i)
5663 points_and_ids[i] = std::make_pair(points[i], i);
5664 const auto p_tree =
pack_rtree(points_and_ids);
5667 std::vector<bool> found_points(points.size(),
false);
5670 const auto already_found = [&found_points](
const auto &id) {
5672 return found_points[
id.second];
5678 const auto store_cell_point_and_id =
5682 const unsigned int &id) {
5683 const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
5684 if (it != cells_out.rend())
5686 const auto cell_id =
5687 (cells_out.size() - 1 - (it - cells_out.rbegin()));
5688 qpoints_out[cell_id].emplace_back(ref_point);
5689 maps_out[cell_id].emplace_back(
id);
5693 cells_out.emplace_back(cell);
5694 qpoints_out.emplace_back(std::vector<
Point<dim>>({ref_point}));
5695 maps_out.emplace_back(std::vector<unsigned int>({
id}));
5700 const auto check_all_points_within_box = [&](
const auto &leaf) {
5701 const double relative_tolerance = 1
e-12;
5704 const auto &cell_hint = leaf.second;
5706 for (
const auto &point_and_id :
5707 p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
5708 bgi::intersects(box)))
5710 const auto id = point_and_id.second;
5711 const auto cell_and_ref =
5715 const auto &cell = cell_and_ref.first;
5716 const auto &ref_point = cell_and_ref.second;
5719 store_cell_point_and_id(cell, ref_point,
id);
5721 missing_points_out.emplace_back(
id);
5724 found_points[id] =
true;
5730 check_all_points_within_box(
5731 std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
5734 for (
unsigned int i = 0; i < np; ++i)
5735 if (found_points[i] ==
false)
5738 const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
5740 if (leaf != b_tree.qend())
5741 check_all_points_within_box(*leaf);
5749 for (
unsigned int i = 0; i < np; ++i)
5750 if (found_points[i] ==
false)
5751 missing_points_out.emplace_back(i);
5758 unsigned int c = cells_out.size();
5759 unsigned int qps = 0;
5765 for (
unsigned int n = 0; n < c; ++n)
5768 qps += qpoints_out[n].size();
5771 Assert(qps + missing_points_out.size() == np,
5775 return std::make_tuple(std::move(cells_out),
5776 std::move(qpoints_out),
5777 std::move(maps_out),
5778 std::move(missing_points_out));
5783 template <
int dim,
int spacedim>
5786 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5787 std::vector<std::vector<Point<dim>>>,
5788 std::vector<std::vector<unsigned int>>,
5789 std::vector<std::vector<Point<spacedim>>>,
5790 std::vector<std::vector<unsigned int>>>
5798 const double tolerance,
5799 const std::vector<bool> & marked_vertices,
5800 const bool enforce_unique_mapping)
5810 enforce_unique_mapping)
5815 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5816 std::vector<std::vector<Point<dim>>>,
5817 std::vector<std::vector<unsigned int>>,
5818 std::vector<std::vector<Point<spacedim>>>,
5819 std::vector<std::vector<unsigned int>>>
5822 std::pair<int, int> dummy{-1, -1};
5824 for (
unsigned int i = 0; i < all.size(); ++i)
5826 if (dummy != std::get<0>(all[i]))
5828 std::get<0>(result).push_back(
5831 std::get<0>(all[i]).
first,
5832 std::get<0>(all[i]).
second});
5834 const unsigned int new_size = std::get<0>(result).size();
5836 std::get<1>(result).resize(new_size);
5837 std::get<2>(result).resize(new_size);
5838 std::get<3>(result).resize(new_size);
5839 std::get<4>(result).resize(new_size);
5841 dummy = std::get<0>(all[i]);
5844 std::get<1>(result).back().push_back(
5845 std::get<3>(all[i]));
5846 std::get<2>(result).back().push_back(std::get<2>(all[i]));
5847 std::get<3>(result).back().push_back(std::get<4>(all[i]));
5848 std::get<4>(result).back().push_back(std::get<1>(all[i]));
5864 template <
int spacedim,
typename T>
5865 std::tuple<std::vector<unsigned int>,
5866 std::vector<unsigned int>,
5867 std::vector<unsigned int>>
5869 const MPI_Comm
comm,
5871 const std::vector<T> & entities,
5872 const double tolerance)
5874 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes_temp;
5875 auto *global_bboxes_to_be_used = &global_bboxes;
5877 if (global_bboxes.size() == 1)
5879 global_bboxes_temp =
5881 global_bboxes_to_be_used = &global_bboxes_temp;
5884 std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
5885 ranks_and_indices.reserve(entities.size());
5890 const auto is_valid = [](
const auto &bb) {
5891 for (
unsigned int i = 0; i < spacedim; ++i)
5892 if (bb.get_boundary_points().first[i] >
5893 bb.get_boundary_points().second[i])
5900 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>>
5903 for (
unsigned rank = 0; rank < global_bboxes_to_be_used->size();
5905 for (
const auto &box : (*global_bboxes_to_be_used)[rank])
5907 boxes_and_ranks.emplace_back(box, rank);
5910 const auto tree =
pack_rtree(boxes_and_ranks);
5913 for (
unsigned int i = 0; i < entities.size(); ++i)
5920 std::set<unsigned int> my_ranks;
5922 for (
const auto &box_and_rank :
5923 tree | boost::geometry::index::adaptors::queried(
5924 boost::geometry::index::intersects(bb)))
5925 my_ranks.insert(box_and_rank.second);
5927 for (
const auto rank : my_ranks)
5928 ranks_and_indices.emplace_back(rank, i);
5937 std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
5939 std::vector<unsigned int> ranks;
5940 std::vector<unsigned int> ptr;
5941 std::vector<unsigned int> indices;
5945 for (
const std::pair<unsigned int, unsigned int> &i : ranks_and_indices)
5947 if (current_rank != i.first)
5949 current_rank = i.first;
5950 ranks.push_back(current_rank);
5951 ptr.push_back(indices.size());
5954 indices.push_back(i.second);
5956 ptr.push_back(indices.size());
5958 return {std::move(ranks), std::move(ptr), std::move(indices)};
5963 template <
int dim,
int spacedim>
5965 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5971 const std::vector<bool> &marked_vertices,
5972 const double tolerance,
5973 const bool enforce_unique_mapping)
5976 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5978 locally_owned_active_cells_around_point;
5995 cell_hint = first_cell.first;
5998 const auto active_cells_around_point =
6006 if (enforce_unique_mapping)
6013 for (
const auto &cell : active_cells_around_point)
6014 lowes_rank =
std::min(lowes_rank, cell.first->subdomain_id());
6016 if (lowes_rank != my_rank)
6020 locally_owned_active_cells_around_point.reserve(
6021 active_cells_around_point.size());
6023 for (
const auto &cell : active_cells_around_point)
6024 if (cell.first->is_locally_owned())
6025 locally_owned_active_cells_around_point.push_back(cell);
6028 std::sort(locally_owned_active_cells_around_point.begin(),
6029 locally_owned_active_cells_around_point.end(),
6030 [](
const auto &a,
const auto &
b) { return a.first < b.first; });
6032 if (enforce_unique_mapping &&
6033 locally_owned_active_cells_around_point.size() > 1)
6035 return {locally_owned_active_cells_around_point.front()};
6037 return locally_owned_active_cells_around_point;
6040 template <
int dim,
int spacedim>
6046 template <
int dim,
int spacedim>
6054 Assert(recv_components.empty() ||
6055 std::get<1>(*std::max_element(recv_components.begin(),
6056 recv_components.end(),
6057 [](
const auto &a,
const auto &
b) {
6058 return std::get<1>(a) <
6060 })) < n_searched_points,
6072 std::sort(send_components.begin(),
6073 send_components.end(),
6074 [&](
const auto &a,
const auto &
b) {
6075 if (std::get<1>(a) != std::get<1>(b))
6076 return std::get<1>(a) < std::get<1>(b);
6078 if (std::get<2>(a) != std::get<2>(b))
6079 return std::get<2>(a) < std::get<2>(b);
6081 return std::get<0>(a) < std::get<0>(b);
6086 i < send_components.size();
6089 std::get<5>(send_components[i]) = i;
6091 if (dummy != std::get<1>(send_components[i]))
6093 dummy = std::get<1>(send_components[i]);
6094 send_ranks.push_back(dummy);
6095 send_ptrs.push_back(i);
6098 send_ptrs.push_back(send_components.size());
6102 std::sort(send_components.begin(),
6103 send_components.end(),
6104 [&](
const auto &a,
const auto &
b) {
6105 if (std::get<0>(a) != std::get<0>(b))
6106 return std::get<0>(a) < std::get<0>(b);
6108 if (std::get<1>(a) != std::get<1>(b))
6109 return std::get<1>(a) < std::get<1>(b);
6111 if (std::get<2>(a) != std::get<2>(b))
6112 return std::get<2>(a) < std::get<2>(b);
6114 return std::get<5>(a) < std::get<5>(b);
6118 if (recv_components.size() > 0)
6121 std::sort(recv_components.begin(),
6122 recv_components.end(),
6123 [&](
const auto &a,
const auto &
b) {
6124 if (std::get<0>(a) != std::get<0>(b))
6125 return std::get<0>(a) < std::get<0>(b);
6127 return std::get<1>(a) < std::get<1>(b);
6132 i < recv_components.size();
6135 std::get<2>(recv_components[i]) = i;
6137 if (dummy != std::get<0>(recv_components[i]))
6139 dummy = std::get<0>(recv_components[i]);
6140 recv_ranks.push_back(dummy);
6141 recv_ptrs.push_back(i);
6144 recv_ptrs.push_back(recv_components.size());
6148 std::sort(recv_components.begin(),
6149 recv_components.end(),
6150 [&](
const auto &a,
const auto &
b) {
6151 if (std::get<1>(a) != std::get<1>(b))
6152 return std::get<1>(a) < std::get<1>(b);
6154 if (std::get<0>(a) != std::get<0>(b))
6155 return std::get<0>(a) < std::get<0>(b);
6157 return std::get<2>(a) < std::get<2>(b);
6164 template <
int dim,
int spacedim>
6170 const std::vector<bool> & marked_vertices,
6171 const double tolerance,
6172 const bool perform_handshake,
6173 const bool enforce_unique_mapping)
6184 comm, global_bboxes, points, tolerance);
6186 const auto &potential_owners_ranks = std::get<0>(potential_owners);
6187 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
6188 const auto &potential_owners_indices = std::get<2>(potential_owners);
6192 const auto translate = [&](
const unsigned int other_rank) {
6193 const auto ptr = std::find(potential_owners_ranks.begin(),
6194 potential_owners_ranks.end(),
6199 const auto other_rank_index =
6200 std::distance(potential_owners_ranks.begin(), ptr);
6202 return other_rank_index;
6206 (marked_vertices.size() == 0) ||
6209 "The marked_vertices vector has to be either empty or its size has "
6210 "to equal the number of vertices of the triangulation."));
6212 using RequestType = std::vector<std::pair<unsigned int, Point<spacedim>>>;
6213 using AnswerType = std::vector<unsigned int>;
6219 const bool has_relevant_vertices =
6220 (marked_vertices.size() == 0) ||
6221 (std::find(marked_vertices.begin(), marked_vertices.end(),
true) !=
6222 marked_vertices.end());
6224 const auto create_request = [&](
const unsigned int other_rank) {
6225 const auto other_rank_index = translate(other_rank);
6227 RequestType request;
6228 request.reserve(potential_owners_ptrs[other_rank_index + 1] -
6229 potential_owners_ptrs[other_rank_index]);
6231 for (
unsigned int i = potential_owners_ptrs[other_rank_index];
6232 i < potential_owners_ptrs[other_rank_index + 1];
6234 request.emplace_back(potential_owners_indices[i],
6235 points[potential_owners_indices[i]]);
6240 const auto answer_request =
6241 [&](
const unsigned int &other_rank,
6242 const RequestType & request) -> AnswerType {
6243 AnswerType answer(request.size(), 0);
6245 if (has_relevant_vertices)
6249 for (
unsigned int i = 0; i < request.size(); ++i)
6251 const auto &index_and_point = request[i];
6253 const auto cells_and_reference_positions =
6256 index_and_point.second,
6260 enforce_unique_mapping);
6265 for (
const auto &cell_and_reference_position :
6266 cells_and_reference_positions)
6268 const auto cell = cell_and_reference_position.first;
6269 auto reference_position =
6270 cell_and_reference_position.second;
6274 if (cell->reference_cell().is_hyper_cube())
6275 reference_position =
6277 reference_position);
6279 send_components.emplace_back(
6280 std::pair<int, int>(cell->level(), cell->index()),
6282 index_and_point.first,
6284 index_and_point.second,
6288 answer[i] = cells_and_reference_positions.size();
6292 if (perform_handshake)
6298 const auto process_answer = [&](
const unsigned int other_rank,
6299 const AnswerType & answer) {
6300 if (perform_handshake)
6302 const auto other_rank_index = translate(other_rank);
6304 for (
unsigned int i = 0; i < answer.size(); ++i)
6305 for (
unsigned int j = 0; j < answer[i]; ++j)
6306 recv_components.emplace_back(
6308 potential_owners_indices
6309 [i + potential_owners_ptrs[other_rank_index]],
6314 Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
6315 potential_owners_ranks,
6328 template <
int structdim,
int spacedim>
6330 DistributedComputePointLocationsInternal<dim, spacedim>
6333 const unsigned int n_quadrature_points,
6336 const bool consistent_numbering_of_sender_and_receiver)
const
6338 using CellIterator =