18 #include <deal.II/base/mpi.templates.h>
63 #include <boost/random/mersenne_twister.hpp>
64 #include <boost/random/uniform_real_distribution.hpp>
74 #include <unordered_map>
81 template <
int dim,
int spacedim>
91 #if defined(DEAL_II_WITH_P4EST) && defined(DEBUG)
106 std::vector<bool> boundary_vertices(
vertices.size(),
false);
112 for (; cell != endc; ++cell)
113 for (
const unsigned int face : cell->face_indices())
114 if (cell->face(face)->at_boundary())
115 for (
unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
116 boundary_vertices[cell->face(face)->vertex_index(i)] =
true;
120 double max_distance_sqr = 0;
121 std::vector<bool>::const_iterator pi = boundary_vertices.begin();
122 const unsigned int N = boundary_vertices.size();
123 for (
unsigned int i = 0; i <
N; ++i, ++pi)
125 std::vector<bool>::const_iterator pj = pi + 1;
126 for (
unsigned int j = i + 1; j <
N; ++j, ++pj)
127 if ((*pi ==
true) && (*pj ==
true) &&
132 return std::sqrt(max_distance_sqr);
137 template <
int dim,
int spacedim>
146 reference_cell.template get_default_linear_mapping<dim, spacedim>());
151 template <
int dim,
int spacedim>
157 unsigned int mapping_degree = 1;
159 mapping_degree = p->get_degree();
160 else if (
const auto *p =
162 mapping_degree = p->get_degree();
169 reference_cell.template get_gauss_type_quadrature<dim>(mapping_degree +
171 const unsigned int n_q_points = quadrature_formula.
size();
182 double local_volume = 0;
185 for (
const auto &cell :
triangulation.active_cell_iterators())
186 if (cell->is_locally_owned())
189 for (
unsigned int q = 0; q < n_q_points; ++q)
190 local_volume += fe_values.
JxW(q);
193 double global_volume = 0;
195 #ifdef DEAL_II_WITH_MPI
203 global_volume = local_volume;
205 return global_volume;
227 struct TransformR2UAffine
242 [1] = {{-1.000000}, {1.000000}};
246 {1.000000, 0.000000};
257 [2] = {{-0.500000, -0.500000},
258 {0.500000, -0.500000},
259 {-0.500000, 0.500000},
260 {0.500000, 0.500000}};
270 {0.750000, 0.250000, 0.250000, -0.250000};
276 {-0.250000, -0.250000, -0.250000},
277 {0.250000, -0.250000, -0.250000},
278 {-0.250000, 0.250000, -0.250000},
279 {0.250000, 0.250000, -0.250000},
280 {-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}
302 template <
int dim,
int spacedim>
311 for (
unsigned int d = 0;
d < spacedim; ++
d)
312 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
313 for (
unsigned int e = 0;
e < dim; ++
e)
318 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
321 return std::make_pair(
A,
b);
338 for (
const auto &cell :
triangulation.active_cell_iterators())
340 if (cell->is_locally_owned())
342 double aspect_ratio_cell = 0.0;
347 for (
unsigned int q = 0; q < quadrature.
size(); ++q)
358 aspect_ratio_cell = std::numeric_limits<double>::infinity();
363 for (
unsigned int i = 0; i < dim; ++i)
364 for (
unsigned int j = 0; j < dim; ++j)
365 J(i, j) = jacobian[i][j];
371 double const ar = max_sv / min_sv;
378 aspect_ratio_cell =
std::max(aspect_ratio_cell, ar);
383 aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
387 return aspect_ratio_vector;
408 template <
int dim,
int spacedim>
414 const auto predicate = [](
const iterator &) {
return true; };
417 tria, std::function<
bool(
const iterator &)>(predicate));
443 template <
int structdim>
467 "Two CellData objects with equal vertices must "
468 "have the same material/boundary ids and manifold "
493 template <
class FaceIteratorType>
497 CellData<dim - 1> face_cell_data(face->n_vertices());
498 for (
unsigned int vertex_n = 0; vertex_n < face->n_vertices();
500 face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
501 face_cell_data.boundary_id = face->boundary_id();
502 face_cell_data.manifold_id = face->manifold_id();
504 face_data.insert(std::move(face_cell_data));
532 template <
class FaceIteratorType>
547 template <
int dim,
int spacedim>
549 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>,
SubCellData>
553 ExcMessage(
"The input triangulation must be non-empty."));
555 std::vector<Point<spacedim>>
vertices;
556 std::vector<CellData<dim>> cells;
558 unsigned int max_level_0_vertex_n = 0;
560 for (
const unsigned int cell_vertex_n : cell->vertex_indices())
561 max_level_0_vertex_n =
562 std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
563 vertices.resize(max_level_0_vertex_n + 1);
573 for (
const unsigned int cell_vertex_n : cell->vertex_indices())
577 vertices[cell->vertex_index(cell_vertex_n)] =
578 cell->vertex(cell_vertex_n);
580 cell->vertex_index(cell_vertex_n);
584 cells.push_back(cell_data);
589 for (
const unsigned int face_n : cell->face_indices())
592 const auto face = cell->face(face_n);
601 for (
unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
603 const auto line = cell->line(line_n);
609 for (
unsigned int vertex_n : line->vertex_indices())
611 line->vertex_index(vertex_n);
614 line_data.insert(std::move(line_cell_data));
623 std::vector<bool> used_vertices(
vertices.size());
625 for (
const auto v : cell_data.vertices)
626 used_vertices[v] =
true;
627 Assert(std::find(used_vertices.begin(), used_vertices.end(),
false) ==
629 ExcMessage(
"The level zero vertices should form a contiguous "
637 for (
const CellData<1> &face_line_data : line_data)
640 return std::tuple<std::vector<Point<spacedim>>,
641 std::vector<CellData<dim>>,
644 std::move(subcell_data));
649 template <
int dim,
int spacedim>
658 "Invalid SubCellData supplied according to ::check_consistency(). "
659 "This is caused by data containing objects for the wrong dimension."));
662 std::vector<bool> vertex_used(
vertices.size(),
false);
663 for (
unsigned int c = 0; c < cells.size(); ++c)
664 for (
unsigned int v = 0; v < cells[c].vertices.size(); ++v)
667 ExcMessage(
"Invalid vertex index encountered! cells[" +
671 " is invalid, because only " +
673 " vertices were supplied."));
674 vertex_used[cells[c].vertices[v]] =
true;
681 std::vector<unsigned int> new_vertex_numbers(
vertices.size(),
683 unsigned int next_free_number = 0;
684 for (
unsigned int i = 0; i <
vertices.size(); ++i)
685 if (vertex_used[i] ==
true)
687 new_vertex_numbers[i] = next_free_number;
692 for (
unsigned int c = 0; c < cells.size(); ++c)
694 v = new_vertex_numbers[v];
699 for (
unsigned int v = 0;
704 new_vertex_numbers.size(),
706 "Invalid vertex index in subcelldata.boundary_lines. "
707 "subcelldata.boundary_lines[" +
712 " is invalid, because only " +
714 " vertices were supplied."));
721 for (
unsigned int v = 0;
726 new_vertex_numbers.size(),
728 "Invalid vertex index in subcelldata.boundary_quads. "
729 "subcelldata.boundary_quads[" +
734 " is invalid, because only " +
736 " vertices were supplied."));
744 std::vector<Point<spacedim>> tmp;
745 tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(),
true));
746 for (
unsigned int v = 0; v <
vertices.size(); ++v)
747 if (vertex_used[v] ==
true)
754 template <
int dim,
int spacedim>
759 std::vector<unsigned int> & considered_vertices,
766 std::vector<unsigned int> new_vertex_numbers(
vertices.size());
767 std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
770 if (considered_vertices.size() == 0)
771 considered_vertices = new_vertex_numbers;
786 unsigned int longest_coordinate_direction = 0;
787 double longest_coordinate_length = bbox.
side_length(0);
788 for (
unsigned int d = 1;
d < spacedim; ++
d)
791 if (longest_coordinate_length < coordinate_length)
793 longest_coordinate_length = coordinate_length;
794 longest_coordinate_direction =
d;
800 std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
801 sorted_vertices.reserve(
vertices.size());
802 for (
const unsigned int vertex_n : considered_vertices)
805 sorted_vertices.emplace_back(vertex_n,
vertices[vertex_n]);
807 std::sort(sorted_vertices.begin(),
808 sorted_vertices.end(),
811 return a.second[longest_coordinate_direction] <
812 b.second[longest_coordinate_direction];
817 for (
unsigned int d = 0;
d < spacedim; ++
d)
818 if (std::abs(a[
d] -
b[
d]) > tol)
825 auto range_start = sorted_vertices.begin();
826 while (range_start != sorted_vertices.end())
828 auto range_end = range_start + 1;
829 while (range_end != sorted_vertices.end() &&
830 std::abs(range_end->second[longest_coordinate_direction] -
831 range_start->second[longest_coordinate_direction]) <
837 std::sort(range_start,
841 return a.first <
b.first;
850 for (
auto reference = range_start; reference != range_end; ++reference)
853 for (
auto it = reference + 1; it != range_end; ++it)
855 if (within_tolerance(reference->second, it->second))
857 new_vertex_numbers[it->first] = reference->first;
863 range_start = range_end;
870 for (
auto &cell : cells)
871 for (
auto &vertex_index : cell.vertices)
872 vertex_index = new_vertex_numbers[vertex_index];
874 for (
auto &vertex_index : quad.vertices)
875 vertex_index = new_vertex_numbers[vertex_index];
877 for (
auto &vertex_index : line.vertices)
878 vertex_index = new_vertex_numbers[vertex_index];
898 for (
unsigned int i = 0; i <
vertices.size(); ++i)
899 map_point_to_local_vertex_index[
vertices[i]] = i;
902 if (map_point_to_local_vertex_index.size() ==
vertices.size())
906 vertices.resize(map_point_to_local_vertex_index.size());
909 for (
const auto &p : map_point_to_local_vertex_index)
916 template <
int dim,
int spacedim>
927 if (dim == 2 && spacedim == 3)
930 std::size_t n_negative_cells = 0;
931 std::size_t cell_no = 0;
932 for (
auto &cell : cells)
948 std::swap(cell.vertices[1], cell.vertices[2]);
953 std::swap(cell.vertices[0], cell.vertices[2]);
954 std::swap(cell.vertices[1], cell.vertices[3]);
955 std::swap(cell.vertices[4], cell.vertices[6]);
956 std::swap(cell.vertices[5], cell.vertices[7]);
964 std::swap(cell.vertices[cell.vertices.size() - 2],
965 cell.vertices[cell.vertices.size() - 1]);
970 std::swap(cell.vertices[0], cell.vertices[3]);
971 std::swap(cell.vertices[1], cell.vertices[4]);
972 std::swap(cell.vertices[2], cell.vertices[5]);
978 std::swap(cell.vertices[2], cell.vertices[3]);
992 return n_negative_cells;
996 template <
int dim,
int spacedim>
1002 const std::size_t n_negative_cells =
1011 AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
1014 "This function assumes that either all cells have positive "
1015 "volume, or that all cells have been specified in an "
1016 "inverted vertex order so that their volume is negative. "
1017 "(In the latter case, this class automatically inverts "
1018 "every cell.) However, the mesh you have specified "
1019 "appears to have both cells with positive and cells with "
1020 "negative volume. You need to check your mesh which "
1021 "cells these are and how they got there.\n"
1022 "As a hint, of the total ") +
1025 " appear to have a negative volume."));
1043 CheapEdge(
const unsigned int v0,
const unsigned int v1)
1055 return ((
v0 <
e.v0) || ((
v0 ==
e.v0) && (
v1 <
e.v1)));
1078 std::set<CheapEdge> edges;
1080 for (
typename std::vector<
CellData<dim>>::const_iterator c =
1090 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1092 const CheapEdge reverse_edge(
1095 if (edges.find(reverse_edge) != edges.end())
1100 const CheapEdge correct_edge(
1103 edges.insert(correct_edge);
1119 struct ParallelEdges
1133 static const unsigned int
1198 class AdjacentCells;
1206 class AdjacentCells<2>
1213 using const_iterator =
const AdjacentCell *;
1224 push_back(
const AdjacentCell &adjacent_cell)
1288 class AdjacentCells<3> :
public std::vector<AdjacentCell>
1310 Edge(
const CellData<dim> &cell,
const unsigned int edge_number)
1361 enum OrientationStatus
1391 Cell(
const CellData<dim> &c,
const std::vector<Edge<dim>> &edge_list)
1398 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1400 const Edge<dim>
e(c,
l);
1436 class EdgeDeltaSet<2>
1442 using const_iterator =
const unsigned int *;
1468 insert(
const unsigned int edge_index)
1529 class EdgeDeltaSet<3> :
public std::set<unsigned int>
1539 std::vector<Edge<dim>>
1546 std::vector<Edge<dim>> edge_list;
1548 for (
unsigned int i = 0; i < cells.size(); ++i)
1549 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1550 edge_list.emplace_back(cells[i],
l);
1553 std::sort(edge_list.begin(), edge_list.end());
1554 edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
1567 std::vector<Cell<dim>>
1568 build_cells_and_connect_edges(
const std::vector<
CellData<dim>> &cells,
1569 std::vector<Edge<dim>> & edges)
1571 std::vector<Cell<dim>> cell_list;
1572 cell_list.reserve(cells.size());
1573 for (
unsigned int i = 0; i < cells.size(); ++i)
1577 cell_list.emplace_back(cells[i], edges);
1581 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1582 edges[cell_list.back().edge_indices[
l]].adjacent_cells.push_back(
1583 AdjacentCell(i,
l));
1598 get_next_unoriented_cell(
const std::vector<Cell<dim>> &cells,
1599 const std::vector<Edge<dim>> &edges,
1600 const unsigned int current_cell)
1602 for (
unsigned int c = current_cell; c < cells.size(); ++c)
1603 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1605 Edge<dim>::not_oriented)
1620 orient_one_set_of_parallel_edges(
const std::vector<Cell<dim>> &cells,
1621 std::vector<Edge<dim>> & edges,
1622 const unsigned int cell,
1623 const unsigned int local_edge)
1650 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1651 (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
1667 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1668 (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
1677 EdgeDeltaSet<dim> Delta_k;
1678 EdgeDeltaSet<dim> Delta_k_minus_1;
1679 Delta_k_minus_1.insert(cells[cell].
edge_indices[local_edge]);
1681 while (Delta_k_minus_1.begin() !=
1682 Delta_k_minus_1.end())
1686 for (
typename EdgeDeltaSet<dim>::const_iterator delta =
1687 Delta_k_minus_1.begin();
1688 delta != Delta_k_minus_1.end();
1692 Edge<dim>::not_oriented,
1696 for (
typename AdjacentCells<dim>::const_iterator adjacent_cell =
1698 adjacent_cell != edges[*delta].adjacent_cells.end();
1701 const unsigned int K = adjacent_cell->cell_index;
1702 const unsigned int delta_is_edge_in_K =
1703 adjacent_cell->edge_within_cell;
1707 const unsigned int first_edge_vertex =
1708 (edges[*delta].orientation_status == Edge<dim>::forward ?
1709 edges[*delta].vertex_indices[0] :
1710 edges[*delta].vertex_indices[1]);
1711 const unsigned int first_edge_vertex_in_K =
1714 delta_is_edge_in_K, 0)];
1716 first_edge_vertex == first_edge_vertex_in_K ||
1717 first_edge_vertex ==
1719 dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
1724 for (
unsigned int o_e = 0;
1731 const unsigned int opposite_edge =
1732 cells[
K].edge_indices[ParallelEdges<
1734 const unsigned int first_opposite_edge_vertex =
1735 cells[
K].vertex_indices
1739 (first_edge_vertex == first_edge_vertex_in_K ? 0 :
1746 const typename Edge<dim>::OrientationStatus
1747 opposite_edge_orientation =
1748 (edges[opposite_edge].vertex_indices[0] ==
1749 first_opposite_edge_vertex ?
1750 Edge<dim>::forward :
1751 Edge<dim>::backward);
1756 Edge<dim>::not_oriented)
1760 edges[opposite_edge].orientation_status =
1761 opposite_edge_orientation;
1762 Delta_k.insert(opposite_edge);
1778 opposite_edge_orientation,
1784 opposite_edge_orientation)
1797 Delta_k_minus_1 = Delta_k;
1811 rotate_cell(
const std::vector<Cell<dim>> &cell_list,
1812 const std::vector<Edge<dim>> &edge_list,
1819 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++
e)
1826 starting_vertex_of_edge[
e] =
1830 starting_vertex_of_edge[
e] =
1846 if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
1847 (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
1848 origin_vertex_of_cell = starting_vertex_of_edge[0];
1849 else if ((starting_vertex_of_edge[1] ==
1850 starting_vertex_of_edge[2]) ||
1851 (starting_vertex_of_edge[1] ==
1852 starting_vertex_of_edge[3]))
1853 origin_vertex_of_cell = starting_vertex_of_edge[1];
1865 for (origin_vertex_of_cell = 0;
1866 origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
1867 ++origin_vertex_of_cell)
1868 if (std::count(starting_vertex_of_edge,
1869 starting_vertex_of_edge +
1874 Assert(origin_vertex_of_cell <
1899 const unsigned int tmp = raw_cells[
cell_index].vertices[0];
1922 static const unsigned int cube_permutations[8][8] = {
1923 {0, 1, 2, 3, 4, 5, 6, 7},
1924 {1, 5, 3, 7, 0, 4, 2, 6},
1925 {2, 6, 0, 4, 3, 7, 1, 5},
1926 {3, 2, 1, 0, 7, 6, 5, 4},
1927 {4, 0, 6, 2, 5, 1, 7, 3},
1928 {5, 4, 7, 6, 1, 0, 3, 2},
1929 {6, 7, 4, 5, 2, 3, 0, 1},
1930 {7, 3, 5, 1, 6, 2, 4, 0}};
1935 temp_vertex_indices[v] =
1937 .vertices[cube_permutations[origin_vertex_of_cell][v]];
1939 raw_cells[
cell_index].vertices[v] = temp_vertex_indices[v];
1963 std::vector<Edge<dim>> edge_list = build_edges(cells);
1964 std::vector<Cell<dim>> cell_list =
1965 build_cells_and_connect_edges(cells, edge_list);
1969 unsigned int next_cell_with_unoriented_edge = 0;
1970 while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
1971 cell_list, edge_list, next_cell_with_unoriented_edge)) !=
1985 for (
unsigned int l = 0;
l < dim; ++
l)
1986 if (edge_list[cell_list[next_cell_with_unoriented_edge]
1989 orient_one_set_of_parallel_edges(
1992 next_cell_with_unoriented_edge,
1997 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1998 Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
2007 for (
unsigned int c = 0; c < cells.size(); ++c)
2008 rotate_cell(cell_list, edge_list, c, cells);
2023 Assert(cells.size() != 0,
2025 "List of elements to orient must have at least one cell"));
2033 if (!is_consistent(cells))
2052 template <
int spacedim>
2110 template <
int spacedim>
2129 template <
int dim,
int spacedim>
2139 template <
int dim,
int spacedim>
2148 "GridTools::rotate() is only available for spacedim = 2."));
2183 const unsigned int axis,
2195 template <
int dim,
int spacedim>
2217 const unsigned int n_dofs = S.
n();
2229 const auto constrained_rhs =
2231 solver.
solve(SF, u, constrained_rhs, prec);
2244 const bool solve_for_absolute_positions)
2270 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
2279 std::array<AffineConstraints<double>, dim> constraints;
2280 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
2288 for (
const unsigned int vertex_no : cell->vertex_indices())
2290 const unsigned int vertex_index = cell->vertex_index(vertex_no);
2291 const Point<dim> & vertex_point = cell->vertex(vertex_no);
2293 const typename std::map<unsigned int, Point<dim>>::const_iterator
2294 map_iter = new_points.find(vertex_index);
2296 if (map_iter != map_end)
2297 for (
unsigned int i = 0; i < dim; ++i)
2299 constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
2300 constraints[i].set_inhomogeneity(
2301 cell->vertex_dof_index(vertex_no, 0),
2302 (solve_for_absolute_positions ?
2303 map_iter->second(i) :
2304 map_iter->second(i) - vertex_point[i]));
2309 for (
unsigned int i = 0; i < dim; ++i)
2310 constraints[i].close();
2314 for (
unsigned int i = 0; i < dim; ++i)
2319 for (
unsigned int i = 0; i < dim; ++i)
2326 std::vector<bool> vertex_touched(
triangulation.n_vertices(),
false);
2328 for (
const unsigned int vertex_no : cell->vertex_indices())
2329 if (vertex_touched[cell->vertex_index(vertex_no)] ==
false)
2334 cell->vertex_dof_index(vertex_no, 0);
2335 for (
unsigned int i = 0; i < dim; ++i)
2336 if (solve_for_absolute_positions)
2337 v(i) = us[i](dof_index);
2339 v(i) += us[i](dof_index);
2341 vertex_touched[cell->vertex_index(vertex_no)] =
true;
2345 template <
int dim,
int spacedim>
2346 std::map<unsigned int, Point<spacedim>>
2349 std::map<unsigned int, Point<spacedim>> vertex_map;
2353 for (; cell != endc; ++cell)
2355 for (
unsigned int i : cell->face_indices())
2359 if (face->at_boundary())
2361 for (
unsigned j = 0; j < face->
n_vertices(); ++j)
2364 const unsigned int vertex_index = face->vertex_index(j);
2365 vertex_map[vertex_index] = vertex;
2377 template <
int dim,
int spacedim>
2381 const bool keep_boundary,
2382 const unsigned int seed)
2399 double almost_infinite_length = 0;
2404 almost_infinite_length += cell->diameter();
2406 std::vector<double> minimal_length(
triangulation.n_vertices(),
2407 almost_infinite_length);
2410 std::vector<bool> at_boundary(keep_boundary ?
triangulation.n_vertices() :
2415 const bool is_parallel_shared =
2418 for (
const auto &cell :
triangulation.active_cell_iterators())
2419 if (is_parallel_shared || cell->is_locally_owned())
2423 for (
unsigned int i = 0; i < cell->n_lines(); ++i)
2426 line = cell->line(i);
2428 if (keep_boundary && line->at_boundary())
2430 at_boundary[line->vertex_index(0)] =
true;
2431 at_boundary[line->vertex_index(1)] =
true;
2434 minimal_length[line->vertex_index(0)] =
2436 minimal_length[line->vertex_index(0)]);
2437 minimal_length[line->vertex_index(1)] =
2439 minimal_length[line->vertex_index(1)]);
2445 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
2446 if (cell->at_boundary(vertex) ==
true)
2447 at_boundary[cell->vertex_index(vertex)] =
true;
2449 minimal_length[cell->vertex_index(0)] =
2451 minimal_length[cell->vertex_index(0)]);
2452 minimal_length[cell->vertex_index(1)] =
2454 minimal_length[cell->vertex_index(1)]);
2459 boost::random::mt19937 rng(seed);
2460 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
2464 if (
auto distributed_triangulation =
2468 const std::vector<bool> locally_owned_vertices =
2470 std::vector<bool> vertex_moved(
triangulation.n_vertices(),
false);
2473 for (
const auto &cell :
triangulation.active_cell_iterators())
2474 if (cell->is_locally_owned())
2476 for (
const unsigned int vertex_no : cell->vertex_indices())
2478 const unsigned global_vertex_no =
2479 cell->vertex_index(vertex_no);
2484 if ((keep_boundary && at_boundary[global_vertex_no]) ||
2485 vertex_moved[global_vertex_no] ||
2486 !locally_owned_vertices[global_vertex_no])
2491 for (
unsigned int d = 0;
d < spacedim; ++
d)
2492 shift_vector(
d) = uniform_distribution(rng);
2494 shift_vector *= factor * minimal_length[global_vertex_no] /
2495 std::sqrt(shift_vector.
square());
2498 cell->vertex(vertex_no) += shift_vector;
2499 vertex_moved[global_vertex_no] =
true;
2503 distributed_triangulation->communicate_locally_moved_vertices(
2504 locally_owned_vertices);
2514 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
2515 const std::vector<Point<spacedim>> &old_vertex_locations =
2518 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2522 if (keep_boundary && at_boundary[vertex])
2523 new_vertex_locations[vertex] = old_vertex_locations[vertex];
2528 for (
unsigned int d = 0;
d < spacedim; ++
d)
2529 shift_vector(
d) = uniform_distribution(rng);
2531 shift_vector *= factor * minimal_length[vertex] /
2532 std::sqrt(shift_vector.
square());
2535 new_vertex_locations[vertex] =
2536 old_vertex_locations[vertex] + shift_vector;
2541 for (
const auto &cell :
triangulation.active_cell_iterators())
2542 for (
const unsigned int vertex_no : cell->vertex_indices())
2543 cell->vertex(vertex_no) =
2544 new_vertex_locations[cell->vertex_index(vertex_no)];
2560 for (; cell != endc; ++cell)
2561 if (!cell->is_artificial())
2562 for (
const unsigned int face : cell->face_indices())
2563 if (cell->face(face)->has_children() &&
2564 !cell->face(face)->at_boundary())
2568 cell->face(face)->child(0)->vertex(1) =
2569 (cell->face(face)->vertex(0) +
2570 cell->face(face)->vertex(1)) /
2574 cell->face(face)->child(0)->vertex(1) =
2575 .5 * (cell->face(face)->vertex(0) +
2576 cell->face(face)->vertex(1));
2577 cell->face(face)->child(0)->vertex(2) =
2578 .5 * (cell->face(face)->vertex(0) +
2579 cell->face(face)->vertex(2));
2580 cell->face(face)->child(1)->vertex(3) =
2581 .5 * (cell->face(face)->vertex(1) +
2582 cell->face(face)->vertex(3));
2583 cell->face(face)->child(2)->vertex(3) =
2584 .5 * (cell->face(face)->vertex(2) +
2585 cell->face(face)->vertex(3));
2588 cell->face(face)->child(0)->vertex(3) =
2589 .25 * (cell->face(face)->vertex(0) +
2590 cell->face(face)->vertex(1) +
2591 cell->face(face)->vertex(2) +
2592 cell->face(face)->vertex(3));
2600 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2604 const
Point<spacedim> & p,
2605 const std::vector<
bool> &marked_vertices)
2614 marked_vertices.size() == 0,
2616 marked_vertices.size()));
2623 marked_vertices.size() == 0 ||
2624 std::equal(marked_vertices.begin(),
2625 marked_vertices.end(),
2627 [](
bool p,
bool q) { return !p || q; }),
2629 "marked_vertices should be a subset of used vertices in the triangulation "
2630 "but marked_vertices contains one or more vertices that are not used vertices!"));
2634 const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
2640 std::vector<bool>::const_iterator
first =
2641 std::find(vertices_to_use.begin(), vertices_to_use.end(),
true);
2646 unsigned int best_vertex = std::distance(vertices_to_use.begin(),
first);
2647 double best_dist = (p -
vertices[best_vertex]).norm_square();
2651 for (
unsigned int j = best_vertex + 1; j <
vertices.size(); ++j)
2652 if (vertices_to_use[j])
2654 const double dist = (p -
vertices[j]).norm_square();
2655 if (dist < best_dist)
2667 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2671 const MeshType<dim, spacedim> &mesh,
2672 const
Point<spacedim> & p,
2673 const std::vector<
bool> &marked_vertices)
2676 if (mapping.preserves_vertex_locations() ==
true)
2686 marked_vertices.size() == 0,
2688 marked_vertices.size()));
2696 marked_vertices.size() == 0 ||
2697 std::equal(marked_vertices.begin(),
2698 marked_vertices.end(),
2700 [](
bool p,
bool q) { return !p || q; }),
2702 "marked_vertices should be a subset of used vertices in the triangulation "
2703 "but marked_vertices contains one or more vertices that are not used vertices!"));
2706 if (marked_vertices.size() != 0)
2709 if (marked_vertices[it->first] ==
false)
2724 template <
int dim,
int spacedim>
2725 std::vector<std::vector<Tensor<1, spacedim>>>
2733 const unsigned int n_vertices = vertex_to_cells.size();
2738 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
2740 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2743 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
2744 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
2746 typename std::set<typename Triangulation<dim, spacedim>::
2747 active_cell_iterator>::iterator it =
2748 vertex_to_cells[vertex].begin();
2749 for (
unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
2751 vertex_to_cell_centers[vertex][cell] =
2752 (*it)->center() -
vertices[vertex];
2753 vertex_to_cell_centers[vertex][cell] /=
2754 vertex_to_cell_centers[vertex][cell].
norm();
2757 return vertex_to_cell_centers;
2763 template <
int spacedim>
2766 const unsigned int a,
2767 const unsigned int b,
2771 const double scalar_product_a = center_directions[a] * point_direction;
2772 const double scalar_product_b = center_directions[
b] * point_direction;
2777 return (scalar_product_a > scalar_product_b);
2781 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2785 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
2787 std::pair<typename ::internal::
2788 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
2793 const MeshType<dim, spacedim> &mesh,
2796 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
2799 &vertex_to_cell_centers,
2800 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
2801 const std::vector<bool> &marked_vertices,
2803 & used_vertices_rtree,
2804 const double tolerance,
2808 *relevant_cell_bounding_boxes_rtree)
2810 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2813 cell_and_position.first = mesh.end();
2817 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2819 cell_and_position_approx;
2821 if (relevant_cell_bounding_boxes_rtree !=
nullptr &&
2822 !relevant_cell_bounding_boxes_rtree->empty())
2828 for (
unsigned int d = 0;
d < spacedim; ++
d)
2830 p1[
d] = p1[
d] - tolerance;
2831 p2[
d] = p2[
d] + tolerance;
2836 if (relevant_cell_bounding_boxes_rtree->qbegin(
2837 boost::geometry::index::intersects(bb)) ==
2838 relevant_cell_bounding_boxes_rtree->qend())
2839 return cell_and_position;
2842 bool found_cell =
false;
2843 bool approx_cell =
false;
2845 unsigned int closest_vertex_index = 0;
2850 if (marked_vertices.size() && !used_vertices_rtree.empty())
2853 std::find(marked_vertices.begin(), marked_vertices.end(),
true);
2854 Assert(itr != marked_vertices.end(),
2856 closest_vertex_index = std::distance(marked_vertices.begin(), itr);
2860 auto current_cell = cell_hint;
2863 const auto cell_marked = [&mesh, &marked_vertices](
const auto &cell) {
2864 if (marked_vertices.size() == 0)
2867 if (cell != mesh.active_cell_iterators().end())
2868 for (
unsigned int i = 0; i < cell->n_vertices(); ++i)
2869 if (marked_vertices[cell->vertex_index(i)])
2876 const auto any_cell_marked = [&cell_marked](
const auto &cells) {
2877 return std::any_of(cells.begin(),
2879 [&cell_marked](
const auto &cell) {
2880 return cell_marked(cell);
2884 while (found_cell ==
false)
2889 cell_marked(cell_hint))
2891 const auto cell_vertices = mapping.
get_vertices(current_cell);
2892 const unsigned int closest_vertex =
2893 find_closest_vertex_of_cell<dim, spacedim>(current_cell,
2896 vertex_to_point = p - cell_vertices[closest_vertex];
2897 closest_vertex_index = current_cell->vertex_index(closest_vertex);
2906 #if defined(DEAL_II_WITH_BOOST_BUNDLED) || \
2907 !(defined(__clang_major__) && __clang_major__ >= 16) || \
2908 BOOST_VERSION >= 108100
2909 if (!used_vertices_rtree.empty())
2912 using ValueType = std::pair<Point<spacedim>,
unsigned int>;
2913 std::function<
bool(
const ValueType &)> marked;
2914 if (marked_vertices.size() == mesh.n_vertices())
2915 marked = [&marked_vertices](
const ValueType &value) ->
bool {
2916 return marked_vertices[value.second];
2919 marked = [](
const ValueType &) ->
bool {
return true; };
2921 std::vector<std::pair<Point<spacedim>,
unsigned int>> res;
2922 used_vertices_rtree.query(
2923 boost::geometry::index::nearest(p, 1) &&
2924 boost::geometry::index::satisfies(marked),
2925 std::back_inserter(res));
2930 ::
ExcMessage(
"There can not be multiple results"));
2933 if (any_cell_marked(vertex_to_cells[res[0].
second]))
2934 closest_vertex_index = res[0].second;
2940 mapping, mesh, p, marked_vertices);
2942 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
2948 Assert(any_cell_marked(vertex_to_cells[closest_vertex_index]),
2953 const double vertex_point_norm = vertex_to_point.
norm();
2954 if (vertex_point_norm > 0)
2955 vertex_to_point /= vertex_point_norm;
2957 const unsigned int n_neighbor_cells =
2958 vertex_to_cells[closest_vertex_index].size();
2961 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
2963 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
2964 neighbor_permutation[i] = i;
2966 auto comp = [&](
const unsigned int a,
const unsigned int b) ->
bool {
2967 return internal::compare_point_association<spacedim>(
2971 vertex_to_cell_centers[closest_vertex_index]);
2974 std::sort(neighbor_permutation.begin(),
2975 neighbor_permutation.end(),
2980 double best_distance = tolerance;
2984 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
2988 auto cell = vertex_to_cells[closest_vertex_index].begin();
2991 if (!(*cell)->is_artificial())
2998 cell_and_position.first = *cell;
2999 cell_and_position.second = p_unit;
3001 approx_cell =
false;
3010 if (dist < best_distance)
3012 best_distance = dist;
3013 cell_and_position_approx.first = *cell;
3014 cell_and_position_approx.second = p_unit;
3023 if (found_cell ==
true)
3024 return cell_and_position;
3025 else if (approx_cell ==
true)
3026 return cell_and_position_approx;
3042 mapping, mesh, p, marked_vertices, tolerance);
3044 current_cell =
typename MeshType<dim, spacedim>::active_cell_iterator();
3046 return cell_and_position;
3051 template <
int dim,
int spacedim>
3060 unsigned int closest_vertex = 0;
3062 for (
unsigned int v = 1; v < cell->
n_vertices(); ++v)
3065 if (vertex_distance < minimum_distance)
3068 minimum_distance = vertex_distance;
3071 return closest_vertex;
3078 namespace BoundingBoxPredicate
3080 template <
class MeshType>
3082 concepts::is_triangulation_or_dof_handler<MeshType>)
3086 cell_iterator &parent_cell,
3087 const std::function<
bool(
3088 const typename MeshType::
3089 active_cell_iterator &)>
3092 bool has_predicate =
3094 std::vector<typename MeshType::active_cell_iterator> active_cells;
3095 if (parent_cell->is_active())
3096 active_cells = {parent_cell};
3100 active_cells = get_active_child_cells<MeshType>(parent_cell);
3102 const unsigned int spacedim = MeshType::space_dimension;
3106 while (i < active_cells.size() && !predicate(active_cells[i]))
3110 if (active_cells.size() == 0 || i == active_cells.size())
3113 return std::make_tuple(bbox, has_predicate);
3120 for (; i < active_cells.size(); ++i)
3121 if (predicate(active_cells[i]))
3123 for (
unsigned int d = 0;
d < spacedim; ++
d)
3125 minp[
d] =
std::min(minp[
d], active_cells[i]->vertex(v)[
d]);
3126 maxp[
d] =
std::max(maxp[
d], active_cells[i]->vertex(v)[
d]);
3129 has_predicate =
true;
3131 return std::make_tuple(bbox, has_predicate);
3138 template <
class MeshType>
3142 const MeshType &mesh,
3143 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
3145 const unsigned int refinement_level,
3146 const bool allow_merge,
3147 const unsigned int max_boxes)
3154 refinement_level <= mesh.n_levels(),
3156 "Error: refinement level is higher then total levels in the triangulation!"));
3158 const unsigned int spacedim = MeshType::space_dimension;
3159 std::vector<BoundingBox<spacedim>> bounding_boxes;
3163 for (
unsigned int i = 0; i < refinement_level; ++i)
3164 for (
const typename MeshType::cell_iterator &cell :
3165 mesh.active_cell_iterators_on_level(i))
3167 bool has_predicate =
false;
3169 std::tie(bbox, has_predicate) =
3171 MeshType>(cell, predicate);
3173 bounding_boxes.push_back(bbox);
3177 for (
const typename MeshType::cell_iterator &cell :
3178 mesh.cell_iterators_on_level(refinement_level))
3180 bool has_predicate =
false;
3182 std::tie(bbox, has_predicate) =
3184 MeshType>(cell, predicate);
3186 bounding_boxes.push_back(bbox);
3191 return bounding_boxes;
3197 std::vector<unsigned int> merged_boxes_idx;
3198 bool found_neighbors =
true;
3203 while (found_neighbors)
3205 found_neighbors =
false;
3206 for (
unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
3208 if (std::find(merged_boxes_idx.begin(),
3209 merged_boxes_idx.end(),
3210 i) == merged_boxes_idx.end())
3211 for (
unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
3212 if (std::find(merged_boxes_idx.begin(),
3213 merged_boxes_idx.end(),
3214 j) == merged_boxes_idx.end() &&
3215 bounding_boxes[i].get_neighbor_type(
3216 bounding_boxes[j]) ==
3219 bounding_boxes[i].merge_with(bounding_boxes[j]);
3220 merged_boxes_idx.push_back(j);
3221 found_neighbors =
true;
3227 std::vector<BoundingBox<spacedim>> merged_b_boxes;
3228 for (
unsigned int i = 0; i < bounding_boxes.size(); ++i)
3229 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
3230 merged_boxes_idx.end())
3231 merged_b_boxes.push_back(bounding_boxes[i]);
3236 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
3238 std::vector<double> volumes;
3239 for (
unsigned int i = 0; i < merged_b_boxes.size(); ++i)
3240 volumes.push_back(merged_b_boxes[i].volume());
3242 while (merged_b_boxes.size() > max_boxes)
3244 unsigned int min_idx =
3245 std::min_element(volumes.begin(), volumes.end()) -
3247 volumes.erase(volumes.begin() + min_idx);
3249 bool not_removed =
true;
3250 for (
unsigned int i = 0;
3251 i < merged_b_boxes.size() && not_removed;
3256 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
3257 merged_b_boxes[min_idx]) ==
3259 merged_b_boxes[i].get_neighbor_type(
3260 merged_b_boxes[min_idx]) ==
3263 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
3264 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
3265 not_removed =
false;
3268 ExcMessage(
"Error: couldn't merge bounding boxes!"));
3271 Assert(merged_b_boxes.size() <= max_boxes,
3273 "Error: couldn't reach target number of bounding boxes!"));
3274 return merged_b_boxes;
3280 template <
int spacedim>
3282 std::tuple<std::vector<std::vector<unsigned int>>,
3283 std::map<unsigned int, unsigned int>,
3284 std::map<unsigned int, std::vector<unsigned int>>>
3292 unsigned int n_procs = global_bboxes.size();
3293 std::vector<std::vector<unsigned int>> point_owners(n_procs);
3294 std::map<unsigned int, unsigned int> map_owners_found;
3295 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3297 unsigned int n_points = points.size();
3298 for (
unsigned int pt = 0; pt < n_points; ++pt)
3301 std::vector<unsigned int> owners_found;
3303 for (
unsigned int rk = 0; rk < n_procs; ++rk)
3306 if (bbox.point_inside(points[pt]))
3308 point_owners[rk].emplace_back(pt);
3309 owners_found.emplace_back(rk);
3313 Assert(owners_found.size() > 0,
3314 ExcMessage(
"No owners found for the point " +
3316 if (owners_found.size() == 1)
3317 map_owners_found[pt] = owners_found[0];
3320 map_owners_guessed[pt] = owners_found;
3323 return std::make_tuple(std::move(point_owners),
3324 std::move(map_owners_found),
3325 std::move(map_owners_guessed));
3328 template <
int spacedim>
3330 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
3331 std::map<unsigned int, unsigned int>,
3332 std::map<unsigned int, std::vector<unsigned int>>>
3340 std::map<unsigned int, std::vector<unsigned int>> point_owners;
3341 std::map<unsigned int, unsigned int> map_owners_found;
3342 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3343 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> search_result;
3345 unsigned int n_points = points.size();
3346 for (
unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
3348 search_result.clear();
3351 covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
3352 std::back_inserter(search_result));
3355 std::set<unsigned int> owners_found;
3357 for (
const auto &rank_bbox : search_result)
3361 const bool pt_inserted = owners_found.insert(pt_n).second;
3363 point_owners[rank_bbox.second].emplace_back(pt_n);
3365 Assert(owners_found.size() > 0,
3366 ExcMessage(
"No owners found for the point " +
3368 if (owners_found.size() == 1)
3369 map_owners_found[pt_n] = *owners_found.begin();
3374 std::back_inserter(map_owners_guessed[pt_n]));
3377 return std::make_tuple(std::move(point_owners),
3378 std::move(map_owners_found),
3379 std::move(map_owners_guessed));
3383 template <
int dim,
int spacedim>
3385 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3389 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3394 for (; cell != endc; ++cell)
3395 for (
const unsigned int i : cell->vertex_indices())
3400 for (; cell != endc; ++cell)
3402 for (
unsigned int i : cell->face_indices())
3404 if ((cell->at_boundary(i) ==
false) &&
3405 (cell->neighbor(i)->is_active()))
3408 adjacent_cell = cell->neighbor(i);
3409 for (
unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
3418 for (
unsigned int i = 0; i < cell->
n_lines(); ++i)
3419 if (cell->line(i)->has_children())
3433 template <
int dim,
int spacedim>
3434 std::map<unsigned int, types::global_vertex_index>
3438 std::map<unsigned int, types::global_vertex_index>
3439 local_to_global_vertex_index;
3441 #ifndef DEAL_II_WITH_MPI
3448 ExcMessage(
"This function does not make any sense "
3449 "for parallel::distributed::Triangulation "
3450 "objects if you do not have MPI enabled."));
3454 using active_cell_iterator =
3456 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
3461 unsigned int max_cellid_size = 0;
3462 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
3464 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
3470 std::set<active_cell_iterator> missing_vert_cells;
3471 std::set<unsigned int> used_vertex_index;
3472 for (
const auto &cell :
triangulation.active_cell_iterators())
3474 if (cell->is_locally_owned())
3476 for (
const unsigned int i : cell->vertex_indices())
3479 for (
const auto &adjacent_cell :
3480 vertex_to_cell[cell->vertex_index(i)])
3481 lowest_subdomain_id =
std::min(lowest_subdomain_id,
3482 adjacent_cell->subdomain_id());
3485 if (lowest_subdomain_id == cell->subdomain_id())
3489 if (used_vertex_index.find(cell->vertex_index(i)) ==
3490 used_vertex_index.end())
3493 local_to_global_vertex_index[cell->vertex_index(i)] =
3498 for (
const auto &adjacent_cell :
3499 vertex_to_cell[cell->vertex_index(i)])
3500 if (adjacent_cell->subdomain_id() !=
3501 cell->subdomain_id())
3505 tmp(adjacent_cell->subdomain_id(),
3506 cell->vertex_index(i));
3507 if (vertices_added.find(tmp) ==
3508 vertices_added.end())
3510 vertices_to_send[adjacent_cell
3513 cell->vertex_index(i),
3514 cell->id().to_string());
3515 if (cell->id().to_string().size() >
3518 cell->id().to_string().size();
3519 vertices_added.insert(tmp);
3522 used_vertex_index.insert(cell->vertex_index(i));
3529 vertices_to_recv[lowest_subdomain_id].insert(
3530 cell->vertex_index(i));
3531 missing_vert_cells.insert(cell);
3538 if (cell->is_ghost())
3540 for (
unsigned int i : cell->face_indices())
3542 if (cell->at_boundary(i) ==
false)
3544 if (cell->neighbor(i)->is_active())
3547 spacedim>::active_cell_iterator
3548 adjacent_cell = cell->neighbor(i);
3549 if ((adjacent_cell->is_locally_owned()))
3552 adjacent_cell->subdomain_id();
3553 if (cell->subdomain_id() < adj_subdomain_id)
3554 for (
unsigned int j = 0;
3555 j < cell->face(i)->n_vertices();
3558 vertices_to_recv[cell->subdomain_id()].insert(
3559 cell->face(i)->vertex_index(j));
3560 missing_vert_cells.insert(cell);
3576 int ierr = MPI_Exscan(&next_index,
3584 for (
auto &global_index_it : local_to_global_vertex_index)
3585 global_index_it.second +=
shift;
3599 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
3600 vertices_to_send.size());
3601 std::vector<MPI_Request> first_requests(vertices_to_send.size());
3605 std::string>>>::iterator
3606 vert_to_send_it = vertices_to_send.begin(),
3607 vert_to_send_end = vertices_to_send.end();
3608 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3609 ++vert_to_send_it, ++i)
3611 int destination = vert_to_send_it->first;
3612 const unsigned int n_vertices = vert_to_send_it->second.size();
3613 const int buffer_size = 2 * n_vertices;
3614 vertices_send_buffers[i].resize(buffer_size);
3617 for (
unsigned int j = 0; j < n_vertices; ++j)
3619 vertices_send_buffers[i][2 * j] =
3620 std::get<0>(vert_to_send_it->second[j]);
3621 vertices_send_buffers[i][2 * j + 1] =
3622 local_to_global_vertex_index[std::get<1>(
3623 vert_to_send_it->second[j])];
3627 ierr = MPI_Isend(vertices_send_buffers[i].data(),
3633 &first_requests[i]);
3638 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
3639 vertices_to_recv.size());
3640 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
3641 vert_to_recv_it = vertices_to_recv.begin(),
3642 vert_to_recv_end = vertices_to_recv.end();
3643 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3644 ++vert_to_recv_it, ++i)
3646 int source = vert_to_recv_it->first;
3647 const unsigned int n_vertices = vert_to_recv_it->second.size();
3648 const int buffer_size = 2 * n_vertices;
3649 vertices_recv_buffers[i].resize(buffer_size);
3652 ierr = MPI_Recv(vertices_recv_buffers[i].data(),
3664 std::vector<std::vector<char>> cellids_send_buffers(
3665 vertices_to_send.size());
3666 std::vector<MPI_Request> second_requests(vertices_to_send.size());
3667 vert_to_send_it = vertices_to_send.begin();
3668 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3669 ++vert_to_send_it, ++i)
3671 int destination = vert_to_send_it->first;
3672 const unsigned int n_vertices = vert_to_send_it->second.size();
3673 const int buffer_size = max_cellid_size * n_vertices;
3674 cellids_send_buffers[i].resize(buffer_size);
3677 unsigned int pos = 0;
3678 for (
unsigned int j = 0; j < n_vertices; ++j)
3680 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
3681 for (
unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
3683 if (k < cell_id.size())
3684 cellids_send_buffers[i][pos] = cell_id[k];
3688 cellids_send_buffers[i][pos] =
'-';
3693 ierr = MPI_Isend(cellids_send_buffers[i].data(),
3699 &second_requests[i]);
3704 std::vector<std::vector<char>> cellids_recv_buffers(
3705 vertices_to_recv.size());
3706 vert_to_recv_it = vertices_to_recv.begin();
3707 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3708 ++vert_to_recv_it, ++i)
3710 int source = vert_to_recv_it->first;
3711 const unsigned int n_vertices = vert_to_recv_it->second.size();
3712 const int buffer_size = max_cellid_size * n_vertices;
3713 cellids_recv_buffers[i].resize(buffer_size);
3716 ierr = MPI_Recv(cellids_recv_buffers[i].data(),
3728 vert_to_recv_it = vertices_to_recv.begin();
3729 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3730 ++i, ++vert_to_recv_it)
3732 for (
unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
3734 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
3736 vertices_recv_buffers[i][2 * j + 1];
3737 const std::string cellid_recv(
3738 &cellids_recv_buffers[i][max_cellid_size * j],
3739 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
3741 typename std::set<active_cell_iterator>::iterator
3742 cell_set_it = missing_vert_cells.begin(),
3743 end_cell_set = missing_vert_cells.end();
3744 for (; (found ==
false) && (cell_set_it != end_cell_set);
3747 typename std::set<active_cell_iterator>::iterator
3749 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
3751 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
3752 for (; candidate_cell != end_cell; ++candidate_cell)
3754 std::string current_cellid =
3755 (*candidate_cell)->id().to_string();
3756 current_cellid.resize(max_cellid_size,
'-');
3757 if (current_cellid.compare(cellid_recv) == 0)
3759 local_to_global_vertex_index
3760 [(*candidate_cell)->vertex_index(local_pos_recv)] =
3772 return local_to_global_vertex_index;
3777 template <
int dim,
int spacedim>
3793 for (
const auto &cell :
triangulation.active_cell_iterators())
3795 const unsigned int index = cell->active_cell_index();
3796 cell_connectivity.
add(index, index);
3797 for (
auto f : cell->face_indices())
3798 if ((cell->at_boundary(f) ==
false) &&
3799 (cell->neighbor(f)->has_children() ==
false))
3801 const unsigned int other_index =
3802 cell->neighbor(f)->active_cell_index();
3803 cell_connectivity.
add(index, other_index);
3804 cell_connectivity.
add(other_index, index);
3811 template <
int dim,
int spacedim>
3817 std::vector<std::vector<unsigned int>> vertex_to_cell(
3819 for (
const auto &cell :
triangulation.active_cell_iterators())
3821 for (
const unsigned int v : cell->vertex_indices())
3822 vertex_to_cell[cell->vertex_index(v)].push_back(
3823 cell->active_cell_index());
3828 for (
const auto &cell :
triangulation.active_cell_iterators())
3830 for (
const unsigned int v : cell->vertex_indices())
3831 for (
unsigned int n = 0;
3832 n < vertex_to_cell[cell->vertex_index(v)].size();
3834 cell_connectivity.
add(cell->active_cell_index(),
3835 vertex_to_cell[cell->vertex_index(v)][n]);
3840 template <
int dim,
int spacedim>
3844 const unsigned int level,
3847 std::vector<std::vector<unsigned int>> vertex_to_cell(
3854 for (
const unsigned int v : cell->vertex_indices())
3855 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
3865 for (
const unsigned int v : cell->vertex_indices())
3866 for (
unsigned int n = 0;
3867 n < vertex_to_cell[cell->vertex_index(v)].size();
3869 cell_connectivity.
add(cell->index(),
3870 vertex_to_cell[cell->vertex_index(v)][n]);
3876 template <
int dim,
int spacedim>
3884 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3885 "are already partitioned implicitly and can not be "
3886 "partitioned again explicitly."));
3888 std::vector<unsigned int> cell_weights;
3898 for (
const auto &cell :
triangulation.active_cell_iterators())
3899 if (cell->is_locally_owned())
3900 cell_weights[cell->active_cell_index()] =
3910 if (
const auto shared_tria =
3914 shared_tria->get_communicator(),
3918 Assert(std::accumulate(cell_weights.begin(),
3920 std::uint64_t(0)) > 0,
3921 ExcMessage(
"The global sum of weights over all active cells "
3922 "is zero. Please verify how you generate weights."));
3934 template <
int dim,
int spacedim>
3937 const std::vector<unsigned int> &cell_weights,
3943 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3944 "are already partitioned implicitly and can not be "
3945 "partitioned again explicitly."));
3949 if (n_partitions == 1)
3951 for (
const auto &cell :
triangulation.active_cell_iterators())
3952 cell->set_subdomain_id(0);
3966 sp_cell_connectivity.
copy_from(cell_connectivity);
3969 sp_cell_connectivity,
3976 template <
int dim,
int spacedim>
3985 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3986 "are already partitioned implicitly and can not be "
3987 "partitioned again explicitly."));
3989 std::vector<unsigned int> cell_weights;
3999 for (
const auto &cell :
triangulation.active_cell_iterators() |
4001 cell_weights[cell->active_cell_index()] =
4011 if (
const auto shared_tria =
4015 shared_tria->get_communicator(),
4019 Assert(std::accumulate(cell_weights.begin(),
4021 std::uint64_t(0)) > 0,
4022 ExcMessage(
"The global sum of weights over all active cells "
4023 "is zero. Please verify how you generate weights."));
4029 cell_connection_graph,
4036 template <
int dim,
int spacedim>
4039 const std::vector<unsigned int> &cell_weights,
4046 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4047 "are already partitioned implicitly and can not be "
4048 "partitioned again explicitly."));
4051 ExcMessage(
"Connectivity graph has wrong size"));
4053 ExcMessage(
"Connectivity graph has wrong size"));
4059 if (n_partitions == 1)
4061 for (
const auto &cell :
triangulation.active_cell_iterators())
4062 cell->set_subdomain_id(0);
4070 std::vector<unsigned int> partition_indices(
triangulation.n_active_cells());
4078 for (
const auto &cell :
triangulation.active_cell_iterators())
4079 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
4091 unsigned int & current_proc_idx,
4092 unsigned int & current_cell_idx,
4094 const unsigned int n_partitions)
4096 if (cell->is_active())
4098 while (current_cell_idx >=
4100 (current_proc_idx + 1) / n_partitions))
4102 cell->set_subdomain_id(current_proc_idx);
4107 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4117 template <
int dim,
int spacedim>
4121 const bool group_siblings)
4125 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4126 "are already partitioned implicitly and can not be "
4127 "partitioned again explicitly."));
4135 if (n_partitions == 1)
4137 for (
const auto &cell :
triangulation.active_cell_iterators())
4138 cell->set_subdomain_id(0);
4144 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
4145 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
4151 coarse_cell_to_p4est_tree_permutation.resize(
triangulation.n_cells(0));
4153 coarse_cell_to_p4est_tree_permutation);
4155 p4est_tree_to_coarse_cell_permutation =
4158 unsigned int current_proc_idx = 0;
4159 unsigned int current_cell_idx = 0;
4164 for (
unsigned int idx = 0; idx <
triangulation.n_cells(0); ++idx)
4166 const unsigned int coarse_cell_idx =
4167 p4est_tree_to_coarse_cell_permutation[idx];
4190 for (; cell != endc; ++cell)
4192 if (cell->is_active())
4194 bool all_children_active =
true;
4195 std::map<unsigned int, unsigned int> map_cpu_n_cells;
4196 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4197 if (!cell->child(n)->is_active())
4199 all_children_active =
false;
4203 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
4205 if (!all_children_active)
4208 unsigned int new_owner = cell->child(0)->subdomain_id();
4209 for (std::map<unsigned int, unsigned int>::iterator it =
4210 map_cpu_n_cells.begin();
4211 it != map_cpu_n_cells.end();
4213 if (it->second > map_cpu_n_cells[new_owner])
4214 new_owner = it->first;
4216 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4217 cell->child(n)->set_subdomain_id(new_owner);
4223 template <
int dim,
int spacedim>
4228 for (
int lvl = n_levels - 1; lvl >= 0; --lvl)
4230 for (
const auto &cell :
triangulation.cell_iterators_on_level(lvl))
4232 if (cell->is_active())
4233 cell->set_level_subdomain_id(cell->subdomain_id());
4236 Assert(cell->child(0)->level_subdomain_id() !=
4239 cell->set_level_subdomain_id(
4240 cell->child(0)->level_subdomain_id());
4252 template <
int dim,
int spacedim>
4257 const std::vector<CellId> & cell_ids,
4258 std::vector<types::subdomain_id> &subdomain_ids)
4260 #ifndef DEAL_II_WITH_P4EST
4263 (void)subdomain_ids;
4267 "You are attempting to use a functionality that is only available "
4268 "if deal.II was configured to use p4est, but cmake did not find a "
4269 "valid p4est library."));
4274 for (
const auto &cell_id : cell_ids)
4277 typename ::internal::p4est::types<dim>::quadrant p4est_cell,
4280 ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
4281 for (
const auto &child_index : cell_id.get_child_indices())
4283 ::internal::p4est::init_quadrant_children<dim>(
4284 p4est_cell, p4est_children);
4286 p4est_children[
static_cast<unsigned int>(child_index)];
4292 const_cast<typename ::internal::p4est::types<dim>::forest
4294 cell_id.get_coarse_cell_id(),
4301 subdomain_ids.push_back(owner);
4308 template <
int spacedim>
4312 const std::vector<CellId> &,
4313 std::vector<types::subdomain_id> &)
4322 template <
int dim,
int spacedim>
4323 std::vector<types::subdomain_id>
4325 const std::vector<CellId> & cell_ids)
4327 std::vector<types::subdomain_id> subdomain_ids;
4328 subdomain_ids.reserve(cell_ids.size());
4337 *parallel_tria =
dynamic_cast<
4351 const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
4352 shared_tria->get_true_subdomain_ids_of_cells();
4354 for (
const auto &cell_id : cell_ids)
4356 const unsigned int active_cell_index =
4357 shared_tria->create_cell_iterator(cell_id)->active_cell_index();
4358 subdomain_ids.push_back(
4359 true_subdomain_ids_of_cells[active_cell_index]);
4366 for (
const auto &cell_id : cell_ids)
4368 subdomain_ids.push_back(
4369 triangulation.create_cell_iterator(cell_id)->subdomain_id());
4373 return subdomain_ids;
4378 template <
int dim,
int spacedim>
4381 std::vector<types::subdomain_id> & subdomain)
4386 for (
const auto &cell :
triangulation.active_cell_iterators())
4387 subdomain[cell->active_cell_index()] = cell->subdomain_id();
4392 template <
int dim,
int spacedim>
4398 unsigned int count = 0;
4399 for (
const auto &cell :
triangulation.active_cell_iterators())
4400 if (cell->subdomain_id() == subdomain)
4408 template <
int dim,
int spacedim>
4413 std::vector<bool> locally_owned_vertices =
4420 if (
const auto *tr =
dynamic_cast<
4423 for (
const auto &cell :
triangulation.active_cell_iterators())
4424 if (cell->is_artificial() ||
4425 (cell->is_ghost() &&
4426 (cell->subdomain_id() < tr->locally_owned_subdomain())))
4427 for (
const unsigned int v : cell->vertex_indices())
4428 locally_owned_vertices[cell->vertex_index(v)] =
false;
4430 return locally_owned_vertices;
4435 template <
int dim,
int spacedim>
4441 for (
const auto &cell :
triangulation.active_cell_iterators())
4442 if (!cell->is_artificial())
4443 min_diameter =
std::min(min_diameter, cell->diameter(mapping));
4445 double global_min_diameter = 0;
4447 #ifdef DEAL_II_WITH_MPI
4451 global_min_diameter =
4455 global_min_diameter = min_diameter;
4457 return global_min_diameter;
4462 template <
int dim,
int spacedim>
4467 double max_diameter = 0.;
4468 for (
const auto &cell :
triangulation.active_cell_iterators())
4469 if (!cell->is_artificial())
4470 max_diameter =
std::max(max_diameter, cell->diameter(mapping));
4472 double global_max_diameter = 0;
4474 #ifdef DEAL_II_WITH_MPI
4478 global_max_diameter =
4482 global_max_diameter = max_diameter;
4484 return global_max_diameter;
4491 namespace FixUpDistortedChildCells
4514 template <
typename Iterator,
int spacedim>
4519 const unsigned int structdim =
4520 Iterator::AccessorType::structure_dimension;
4521 Assert(spacedim == Iterator::AccessorType::dimension,
4527 Assert(object->refinement_case() ==
4535 Tensor<spacedim - structdim, spacedim>
4538 for (
const unsigned int i : object->vertex_indices())
4539 parent_vertices[i] =
object->vertex(i);
4542 parent_vertices, parent_alternating_forms);
4544 const Tensor<spacedim - structdim, spacedim>
4545 average_parent_alternating_form =
4546 std::accumulate(parent_alternating_forms,
4547 parent_alternating_forms +
4561 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4565 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4566 for (
const unsigned int i : object->child(c)->vertex_indices())
4567 child_vertices[c][i] =
object->child(c)->vertex(i);
4575 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4577 1] = object_mid_point;
4579 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4581 child_vertices[c], child_alternating_forms[c]);
4593 double objective = 0;
4594 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4595 for (
const unsigned int i : object->child(c)->vertex_indices())
4597 (child_alternating_forms[c][i] -
4598 average_parent_alternating_form / std::pow(2., 1. * structdim))
4610 template <
typename Iterator>
4613 const unsigned int f,
4614 std::integral_constant<int, 1>)
4616 return object->vertex(f);
4626 template <
typename Iterator>
4629 const unsigned int f,
4630 std::integral_constant<int, 2>)
4632 return object->line(f)->center();
4642 template <
typename Iterator>
4645 const unsigned int f,
4646 std::integral_constant<int, 3>)
4648 return object->face(f)->center();
4675 template <
typename Iterator>
4679 const unsigned int structdim =
4680 Iterator::AccessorType::structure_dimension;
4682 double diameter =
object->diameter();
4683 for (
const unsigned int f : object->face_indices())
4684 for (
unsigned int e = f + 1;
e <
object->n_faces(); ++
e)
4689 std::integral_constant<int, structdim>())
4691 object,
e, std::integral_constant<int, structdim>())));
4702 template <
typename Iterator>
4706 const unsigned int structdim =
4707 Iterator::AccessorType::structure_dimension;
4708 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
4715 Assert(object->refinement_case() ==
4725 unsigned int iteration = 0;
4730 double initial_delta = 0;
4737 const double step_length =
diameter / 4 / (iteration + 1);
4742 for (
unsigned int d = 0;
d < spacedim; ++
d)
4744 const double eps = step_length / 10;
4758 if (gradient.norm() == 0)
4767 std::min(2 * current_value / (gradient * gradient),
4768 step_length / gradient.norm()) *
4773 const double previous_value = current_value;
4777 initial_delta = (previous_value - current_value);
4780 if ((iteration >= 1) &&
4781 ((previous_value - current_value < 0) ||
4782 (
std::fabs(previous_value - current_value) <
4783 0.001 * initial_delta)))
4788 while (iteration < 20);
4804 double old_min_product, new_min_product;
4809 parent_vertices[i] = object->vertex(i);
4811 Tensor<spacedim - structdim, spacedim>
4814 parent_vertices, parent_alternating_forms);
4820 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4821 for (
const unsigned int i : object->child(c)->vertex_indices())
4822 child_vertices[c][i] =
object->child(c)->vertex(i);
4824 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4828 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4830 child_vertices[c], child_alternating_forms[c]);
4833 child_alternating_forms[0][0] * parent_alternating_forms[0];
4834 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4835 for (
const unsigned int i : object->child(c)->vertex_indices())
4836 for (
const unsigned int j : object->vertex_indices())
4837 old_min_product = std::min<double>(old_min_product,
4838 child_alternating_forms[c][i] *
4839 parent_alternating_forms[j]);
4847 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4849 1] = object_mid_point;
4851 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4853 child_vertices[c], child_alternating_forms[c]);
4856 child_alternating_forms[0][0] * parent_alternating_forms[0];
4857 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4858 for (
const unsigned int i : object->child(c)->vertex_indices())
4859 for (
const unsigned int j : object->vertex_indices())
4860 new_min_product = std::min<double>(new_min_product,
4861 child_alternating_forms[c][i] *
4862 parent_alternating_forms[j]);
4870 if (new_min_product >= old_min_product)
4871 object->child(0)->vertex(
4878 return (
std::max(new_min_product, old_min_product) > 0);
4884 template <
int dim,
int spacedim>
4887 const typename ::Triangulation<dim, spacedim>::cell_iterator
4889 std::integral_constant<int, dim>,
4890 std::integral_constant<int, spacedim>)
4900 for (
auto f : cell->face_indices())
4903 Assert(cell->face(f)->refinement_case() ==
4907 bool subface_is_more_refined =
false;
4908 for (
unsigned int g = 0;
4909 g < GeometryInfo<dim>::max_children_per_face;
4911 if (cell->face(f)->child(g)->has_children())
4913 subface_is_more_refined =
true;
4917 if (subface_is_more_refined ==
true)
4928 template <
int dim,
int spacedim>
4936 dim != 1 && spacedim != 1,
4937 "This function is only valid when dim != 1 or spacedim != 1.");
4941 for (
typename std::list<
4943 cell_ptr = distorted_cells.distorted_cells.
begin();
4944 cell_ptr != distorted_cells.distorted_cells.
end();
4950 Assert(!cell->is_active(),
4952 "This function is only valid for a list of cells that "
4953 "have children (i.e., no cell in the list may be active)."));
4957 std::integral_constant<int, dim>(),
4958 std::integral_constant<int, spacedim>());
4962 unfixable_subset.distorted_cells.push_back(cell);
4965 return unfixable_subset;
4970 template <
int dim,
int spacedim>
4973 const bool reset_boundary_ids)
4976 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
4977 auto m_it = dst_manifold_ids.begin();
4978 for (
const auto b : src_boundary_ids)
4983 const std::vector<types::boundary_id> reset_boundary_id =
4984 reset_boundary_ids ?
4985 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
4995 template <
int dim,
int spacedim>
4998 const std::vector<types::boundary_id> &src_boundary_ids,
4999 const std::vector<types::manifold_id> &dst_manifold_ids,
5001 const std::vector<types::boundary_id> &reset_boundary_ids_)
5004 const auto reset_boundary_ids =
5005 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
5014 for (
auto f : cell->face_indices())
5015 if (cell->face(f)->at_boundary())
5016 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
5018 const auto bid = cell->face(f)->line(
e)->boundary_id();
5019 const unsigned int ind = std::find(src_boundary_ids.begin(),
5020 src_boundary_ids.end(),
5022 src_boundary_ids.begin();
5023 if (ind < src_boundary_ids.size())
5024 cell->face(f)->line(
e)->set_manifold_id(
5025 dst_manifold_ids[ind]);
5030 for (
auto f : cell->face_indices())
5031 if (cell->face(f)->at_boundary())
5033 const auto bid = cell->face(f)->boundary_id();
5034 const unsigned int ind =
5035 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
5036 src_boundary_ids.begin();
5038 if (ind < src_boundary_ids.size())
5041 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
5043 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
5047 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
5049 const auto bid = cell->face(f)->line(
e)->boundary_id();
5050 const unsigned int ind = std::find(src_boundary_ids.begin(),
5051 src_boundary_ids.end(),
5053 src_boundary_ids.begin();
5054 if (ind < src_boundary_ids.size())
5055 cell->face(f)->line(
e)->set_boundary_id(
5056 reset_boundary_ids[ind]);
5062 template <
int dim,
int spacedim>
5065 const bool compute_face_ids)
5071 for (; cell != endc; ++cell)
5073 cell->set_manifold_id(cell->material_id());
5074 if (compute_face_ids ==
true)
5076 for (
auto f : cell->face_indices())
5078 if (cell->at_boundary(f) ==
false)
5079 cell->face(f)->set_manifold_id(
5081 cell->neighbor(f)->material_id()));
5083 cell->face(f)->set_manifold_id(cell->material_id());
5090 template <
int dim,
int spacedim>
5095 const std::set<types::manifold_id> &)> &disambiguation_function,
5096 bool overwrite_only_flat_manifold_ids)
5101 const unsigned int n_subobjects =
5105 std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
5106 std::vector<unsigned int> backup;
5110 unsigned next_index = 1;
5114 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
5116 if (cell->line(
l)->user_index() == 0)
5119 manifold_ids[next_index].insert(cell->manifold_id());
5120 cell->line(
l)->set_user_index(next_index++);
5123 manifold_ids[cell->line(
l)->user_index()].insert(
5124 cell->manifold_id());
5127 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
5129 if (cell->quad(
l)->user_index() == 0)
5132 manifold_ids[next_index].insert(cell->manifold_id());
5133 cell->quad(
l)->set_user_index(next_index++);
5136 manifold_ids[cell->quad(
l)->user_index()].insert(
5137 cell->manifold_id());
5143 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
5145 const auto id = cell->line(
l)->user_index();
5149 if (cell->line(
l)->manifold_id() ==
5151 overwrite_only_flat_manifold_ids ==
false)
5152 cell->line(
l)->set_manifold_id(
5153 disambiguation_function(manifold_ids[
id]));
5154 cell->line(
l)->set_user_index(0);
5158 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
5160 const auto id = cell->quad(
l)->user_index();
5164 if (cell->quad(
l)->manifold_id() ==
5166 overwrite_only_flat_manifold_ids ==
false)
5167 cell->quad(
l)->set_manifold_id(
5168 disambiguation_function(manifold_ids[
id]));
5169 cell->quad(
l)->set_user_index(0);
5178 template <
int dim,
int spacedim>
5179 std::pair<unsigned int, double>
5183 double max_ratio = 1;
5184 unsigned int index = 0;
5186 for (
unsigned int i = 0; i < dim; ++i)
5187 for (
unsigned int j = i + 1; j < dim; ++j)
5189 unsigned int ax = i % dim;
5190 unsigned int next_ax = j % dim;
5193 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
5195 if (ratio > max_ratio)
5200 else if (1.0 / ratio > max_ratio)
5202 max_ratio = 1.0 / ratio;
5206 return std::make_pair(index, max_ratio);
5210 template <
int dim,
int spacedim>
5213 const bool isotropic,
5214 const unsigned int max_iterations)
5216 unsigned int iter = 0;
5217 bool continue_refinement =
true;
5219 while (continue_refinement && (iter < max_iterations))
5223 continue_refinement =
false;
5226 for (
const unsigned int j : cell->face_indices())
5227 if (cell->at_boundary(j) ==
false &&
5228 cell->neighbor(j)->has_children())
5232 cell->set_refine_flag();
5233 continue_refinement =
true;
5236 continue_refinement |= cell->flag_for_face_refinement(j);
5243 template <
int dim,
int spacedim>
5246 const double max_ratio,
5247 const unsigned int max_iterations)
5249 unsigned int iter = 0;
5250 bool continue_refinement =
true;
5252 while (continue_refinement && (iter < max_iterations))
5255 continue_refinement =
false;
5258 std::pair<unsigned int, double> info =
5259 GridTools::get_longest_direction<dim, spacedim>(cell);
5260 if (info.second > max_ratio)
5262 cell->set_refine_flag(
5264 continue_refinement =
true;
5272 template <
int dim,
int spacedim>
5275 const double limit_angle_fraction)
5283 "have hanging nodes."));
5287 bool has_cells_with_more_than_dim_faces_on_boundary =
true;
5288 bool has_cells_with_dim_faces_on_boundary =
false;
5290 unsigned int refinement_cycles = 0;
5292 while (has_cells_with_more_than_dim_faces_on_boundary)
5294 has_cells_with_more_than_dim_faces_on_boundary =
false;
5298 unsigned int boundary_face_counter = 0;
5299 for (
auto f : cell->face_indices())
5300 if (cell->face(f)->at_boundary())
5301 boundary_face_counter++;
5302 if (boundary_face_counter > dim)
5304 has_cells_with_more_than_dim_faces_on_boundary =
true;
5307 else if (boundary_face_counter == dim)
5308 has_cells_with_dim_faces_on_boundary =
true;
5310 if (has_cells_with_more_than_dim_faces_on_boundary)
5313 refinement_cycles++;
5317 if (has_cells_with_dim_faces_on_boundary)
5320 refinement_cycles++;
5324 while (refinement_cycles > 0)
5327 cell->set_coarsen_flag();
5329 refinement_cycles--;
5339 std::vector<CellData<dim>> cells_to_add;
5343 const unsigned int v0 = 0,
v1 = 1, v2 = (dim > 1 ? 2 : 0),
5344 v3 = (dim > 1 ? 3 : 0);
5348 double angle_fraction = 0;
5354 p0[spacedim > 1 ? 1 : 0] = 1;
5358 if (cell->face(
v0)->at_boundary() && cell->face(v3)->at_boundary())
5360 p0 = cell->vertex(
v0) - cell->vertex(v2);
5361 p1 = cell->vertex(v3) - cell->vertex(v2);
5362 vertex_at_corner = v2;
5364 else if (cell->face(v3)->at_boundary() &&
5365 cell->face(
v1)->at_boundary())
5367 p0 = cell->vertex(v2) - cell->vertex(v3);
5368 p1 = cell->vertex(
v1) - cell->vertex(v3);
5369 vertex_at_corner = v3;
5371 else if (cell->face(1)->at_boundary() &&
5372 cell->face(2)->at_boundary())
5374 p0 = cell->vertex(
v0) - cell->vertex(
v1);
5375 p1 = cell->vertex(v3) - cell->vertex(
v1);
5376 vertex_at_corner =
v1;
5378 else if (cell->face(2)->at_boundary() &&
5379 cell->face(0)->at_boundary())
5381 p0 = cell->vertex(v2) - cell->vertex(
v0);
5382 p1 = cell->vertex(
v1) - cell->vertex(
v0);
5383 vertex_at_corner =
v0;
5394 if (angle_fraction > limit_angle_fraction)
5396 auto flags_removal = [&](
unsigned int f1,
5399 unsigned int n2) ->
void {
5400 cells_to_remove[cell->active_cell_index()] =
true;
5401 cells_to_remove[cell->neighbor(n1)->active_cell_index()] =
true;
5402 cells_to_remove[cell->neighbor(n2)->active_cell_index()] =
true;
5404 faces_to_remove[cell->face(f1)->index()] =
true;
5405 faces_to_remove[cell->face(f2)->index()] =
true;
5407 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] =
true;
5408 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] =
true;
5411 auto cell_creation = [&](
const unsigned int vv0,
5412 const unsigned int vv1,
5413 const unsigned int f0,
5414 const unsigned int f1,
5416 const unsigned int n0,
5417 const unsigned int v0n0,
5418 const unsigned int v1n0,
5420 const unsigned int n1,
5421 const unsigned int v0n1,
5422 const unsigned int v1n1) {
5428 c1.
vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
5429 c1.
vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
5435 c2.
vertices[
v1] = cell->neighbor(n1)->vertex_index(v0n1);
5436 c2.
vertices[v2] = cell->vertex_index(vv1);
5437 c2.
vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
5442 l1.
vertices[0] = cell->vertex_index(vv0);
5443 l1.
vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
5449 l2.
vertices[0] = cell->vertex_index(vv0);
5450 l2.
vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
5456 cells_to_add.push_back(c1);
5457 cells_to_add.push_back(c2);
5462 switch (vertex_at_corner)
5465 flags_removal(0, 2, 3, 1);
5466 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
5469 flags_removal(1, 2, 3, 0);
5470 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
5473 flags_removal(3, 0, 1, 2);
5474 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
5477 flags_removal(3, 1, 0, 2);
5478 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
5491 if (cells_to_add.size() == 0)
5493 while (refinement_cycles > 0)
5496 cell->set_coarsen_flag();
5498 refinement_cycles--;
5506 if (cells_to_remove[cell->active_cell_index()] ==
false)
5509 for (
const unsigned int v : cell->vertex_indices())
5510 c.
vertices[v] = cell->vertex_index(v);
5513 cells_to_add.push_back(c);
5521 for (; face != endf; ++face)
5522 if ((face->at_boundary() ||
5524 faces_to_remove[face->index()] ==
false)
5526 for (
unsigned int l = 0;
l < face->
n_lines(); ++
l)
5531 for (
const unsigned int v : face->vertex_indices())
5532 line.
vertices[v] = face->vertex_index(v);
5538 for (
const unsigned int v : face->line(
l)->vertex_indices())
5539 line.
vertices[v] = face->line(
l)->vertex_index(v);
5548 for (
const unsigned int v : face->vertex_indices())
5549 quad.
vertices[v] = face->vertex_index(v);
5557 subcelldata_to_add);
5562 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
5581 template <
int dim,
int spacedim>
5584 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5585 std::vector<std::vector<Point<dim>>>,
5586 std::vector<std::vector<unsigned int>>>
5598 auto &cells = std::get<0>(cqmp);
5599 auto &qpoints = std::get<1>(cqmp);
5600 auto &maps = std::get<2>(cqmp);
5602 return std::make_tuple(std::move(cells),
5609 template <
int dim,
int spacedim>
5612 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5613 std::vector<std::vector<Point<dim>>>,
5614 std::vector<std::vector<unsigned int>>,
5615 std::vector<unsigned int>>
5625 Assert((dim == spacedim),
5626 ExcMessage(
"Only implemented for dim==spacedim."));
5635 const unsigned int np = points.size();
5637 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5639 std::vector<std::vector<Point<dim>>> qpoints_out;
5640 std::vector<std::vector<unsigned int>> maps_out;
5641 std::vector<unsigned int> missing_points_out;
5645 return std::make_tuple(std::move(cells_out),
5646 std::move(qpoints_out),
5647 std::move(maps_out),
5648 std::move(missing_points_out));
5656 std::vector<std::pair<Point<spacedim>,
unsigned int>> points_and_ids(np);
5657 for (
unsigned int i = 0; i < np; ++i)
5658 points_and_ids[i] = std::make_pair(points[i], i);
5659 const auto p_tree =
pack_rtree(points_and_ids);
5662 std::vector<bool> found_points(points.size(),
false);
5665 const auto already_found = [&found_points](
const auto &id) {
5667 return found_points[
id.second];
5673 const auto store_cell_point_and_id =
5677 const unsigned int &id) {
5678 const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
5679 if (it != cells_out.rend())
5681 const auto cell_id =
5682 (cells_out.size() - 1 - (it - cells_out.rbegin()));
5683 qpoints_out[cell_id].emplace_back(ref_point);
5684 maps_out[cell_id].emplace_back(
id);
5688 cells_out.emplace_back(cell);
5689 qpoints_out.emplace_back(std::vector<
Point<dim>>({ref_point}));
5690 maps_out.emplace_back(std::vector<unsigned int>({
id}));
5695 const auto check_all_points_within_box = [&](
const auto &leaf) {
5696 const double relative_tolerance = 1
e-12;
5699 const auto &cell_hint = leaf.second;
5701 for (
const auto &point_and_id :
5702 p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
5703 bgi::intersects(box)))
5705 const auto id = point_and_id.second;
5706 const auto cell_and_ref =
5710 const auto &cell = cell_and_ref.first;
5711 const auto &ref_point = cell_and_ref.second;
5714 store_cell_point_and_id(cell, ref_point,
id);
5716 missing_points_out.emplace_back(
id);
5719 found_points[id] =
true;
5725 check_all_points_within_box(
5726 std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
5729 for (
unsigned int i = 0; i < np; ++i)
5730 if (found_points[i] ==
false)
5733 const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
5735 if (leaf != b_tree.qend())
5736 check_all_points_within_box(*leaf);
5744 for (
unsigned int i = 0; i < np; ++i)
5745 if (found_points[i] ==
false)
5746 missing_points_out.emplace_back(i);
5753 unsigned int c = cells_out.size();
5754 unsigned int qps = 0;
5760 for (
unsigned int n = 0; n < c; ++n)
5763 qps += qpoints_out[n].size();
5766 Assert(qps + missing_points_out.size() == np,
5770 return std::make_tuple(std::move(cells_out),
5771 std::move(qpoints_out),
5772 std::move(maps_out),
5773 std::move(missing_points_out));
5778 template <
int dim,
int spacedim>
5781 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5782 std::vector<std::vector<Point<dim>>>,
5783 std::vector<std::vector<unsigned int>>,
5784 std::vector<std::vector<Point<spacedim>>>,
5785 std::vector<std::vector<unsigned int>>>
5793 const double tolerance,
5794 const std::vector<bool> & marked_vertices,
5795 const bool enforce_unique_mapping)
5805 enforce_unique_mapping)
5810 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5811 std::vector<std::vector<Point<dim>>>,
5812 std::vector<std::vector<unsigned int>>,
5813 std::vector<std::vector<Point<spacedim>>>,
5814 std::vector<std::vector<unsigned int>>>
5817 std::pair<int, int> dummy{-1, -1};
5819 for (
unsigned int i = 0; i < all.size(); ++i)
5821 if (dummy != std::get<0>(all[i]))
5823 std::get<0>(result).push_back(
5826 std::get<0>(all[i]).
first,
5827 std::get<0>(all[i]).
second});
5829 const unsigned int new_size = std::get<0>(result).size();
5831 std::get<1>(result).resize(new_size);
5832 std::get<2>(result).resize(new_size);
5833 std::get<3>(result).resize(new_size);
5834 std::get<4>(result).resize(new_size);
5836 dummy = std::get<0>(all[i]);
5839 std::get<1>(result).back().push_back(
5840 std::get<3>(all[i]));
5841 std::get<2>(result).back().push_back(std::get<2>(all[i]));
5842 std::get<3>(result).back().push_back(std::get<4>(all[i]));
5843 std::get<4>(result).back().push_back(std::get<1>(all[i]));
5859 template <
int spacedim,
typename T>
5860 std::tuple<std::vector<unsigned int>,
5861 std::vector<unsigned int>,
5862 std::vector<unsigned int>>
5864 const MPI_Comm
comm,
5866 const std::vector<T> & entities,
5867 const double tolerance)
5869 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes_temp;
5870 auto *global_bboxes_to_be_used = &global_bboxes;
5872 if (global_bboxes.size() == 1)
5874 global_bboxes_temp =
5876 global_bboxes_to_be_used = &global_bboxes_temp;
5879 std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
5880 ranks_and_indices.reserve(entities.size());
5885 const auto is_valid = [](
const auto &bb) {
5886 for (
unsigned int i = 0; i < spacedim; ++i)
5887 if (bb.get_boundary_points().first[i] >
5888 bb.get_boundary_points().second[i])
5895 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>>
5898 for (
unsigned rank = 0; rank < global_bboxes_to_be_used->size();
5900 for (
const auto &box : (*global_bboxes_to_be_used)[rank])
5902 boxes_and_ranks.emplace_back(box, rank);
5905 const auto tree =
pack_rtree(boxes_and_ranks);
5908 for (
unsigned int i = 0; i < entities.size(); ++i)
5915 std::set<unsigned int> my_ranks;
5917 for (
const auto &box_and_rank :
5918 tree | boost::geometry::index::adaptors::queried(
5919 boost::geometry::index::intersects(bb)))
5920 my_ranks.insert(box_and_rank.second);
5922 for (
const auto rank : my_ranks)
5923 ranks_and_indices.emplace_back(rank, i);
5932 std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
5934 std::vector<unsigned int> ranks;
5935 std::vector<unsigned int> ptr;
5936 std::vector<unsigned int> indices;
5940 for (
const std::pair<unsigned int, unsigned int> &i : ranks_and_indices)
5942 if (current_rank != i.first)
5944 current_rank = i.first;
5945 ranks.push_back(current_rank);
5946 ptr.push_back(indices.size());
5949 indices.push_back(i.second);
5951 ptr.push_back(indices.size());
5953 return {std::move(ranks), std::move(ptr), std::move(indices)};
5958 template <
int dim,
int spacedim>
5960 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5966 const std::vector<bool> &marked_vertices,
5967 const double tolerance,
5968 const bool enforce_unique_mapping)
5971 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5973 locally_owned_active_cells_around_point;
5990 cell_hint = first_cell.first;
5993 const auto active_cells_around_point =
6001 if (enforce_unique_mapping)
6008 for (
const auto &cell : active_cells_around_point)
6009 lowes_rank =
std::min(lowes_rank, cell.first->subdomain_id());
6011 if (lowes_rank != my_rank)
6015 locally_owned_active_cells_around_point.reserve(
6016 active_cells_around_point.size());
6018 for (
const auto &cell : active_cells_around_point)
6019 if (cell.first->is_locally_owned())
6020 locally_owned_active_cells_around_point.push_back(cell);
6023 std::sort(locally_owned_active_cells_around_point.begin(),
6024 locally_owned_active_cells_around_point.end(),
6025 [](
const auto &a,
const auto &
b) { return a.first < b.first; });
6027 if (enforce_unique_mapping &&
6028 locally_owned_active_cells_around_point.size() > 1)
6030 return {locally_owned_active_cells_around_point.front()};
6032 return locally_owned_active_cells_around_point;
6035 template <
int dim,
int spacedim>
6041 template <
int dim,
int spacedim>
6049 Assert(recv_components.empty() ||
6050 std::get<1>(*std::max_element(recv_components.begin(),
6051 recv_components.end(),
6052 [](
const auto &a,
const auto &
b) {
6053 return std::get<1>(a) <
6055 })) < n_searched_points,
6067 std::sort(send_components.begin(),
6068 send_components.end(),
6069 [&](
const auto &a,
const auto &
b) {
6070 if (std::get<1>(a) != std::get<1>(b))
6071 return std::get<1>(a) < std::get<1>(b);
6073 if (std::get<2>(a) != std::get<2>(b))
6074 return std::get<2>(a) < std::get<2>(b);
6076 return std::get<0>(a) < std::get<0>(b);
6081 i < send_components.size();
6084 std::get<5>(send_components[i]) = i;
6086 if (dummy != std::get<1>(send_components[i]))
6088 dummy = std::get<1>(send_components[i]);
6089 send_ranks.push_back(dummy);
6090 send_ptrs.push_back(i);
6093 send_ptrs.push_back(send_components.size());
6097 std::sort(send_components.begin(),
6098 send_components.end(),
6099 [&](
const auto &a,
const auto &
b) {
6100 if (std::get<0>(a) != std::get<0>(b))
6101 return std::get<0>(a) < std::get<0>(b);
6103 if (std::get<1>(a) != std::get<1>(b))
6104 return std::get<1>(a) < std::get<1>(b);
6106 if (std::get<2>(a) != std::get<2>(b))
6107 return std::get<2>(a) < std::get<2>(b);
6109 return std::get<5>(a) < std::get<5>(b);
6113 if (recv_components.size() > 0)
6116 std::sort(recv_components.begin(),
6117 recv_components.end(),
6118 [&](
const auto &a,
const auto &
b) {
6119 if (std::get<0>(a) != std::get<0>(b))
6120 return std::get<0>(a) < std::get<0>(b);
6122 return std::get<1>(a) < std::get<1>(b);
6127 i < recv_components.size();
6130 std::get<2>(recv_components[i]) = i;
6132 if (dummy != std::get<0>(recv_components[i]))
6134 dummy = std::get<0>(recv_components[i]);
6135 recv_ranks.push_back(dummy);
6136 recv_ptrs.push_back(i);
6139 recv_ptrs.push_back(recv_components.size());
6143 std::sort(recv_components.begin(),
6144 recv_components.end(),
6145 [&](
const auto &a,
const auto &
b) {
6146 if (std::get<1>(a) != std::get<1>(b))
6147 return std::get<1>(a) < std::get<1>(b);
6149 if (std::get<0>(a) != std::get<0>(b))
6150 return std::get<0>(a) < std::get<0>(b);
6152 return std::get<2>(a) < std::get<2>(b);
6158 template <
int dim,
int spacedim>
6164 const std::vector<bool> & marked_vertices,
6165 const double tolerance,
6166 const bool perform_handshake,
6167 const bool enforce_unique_mapping)
6178 comm, global_bboxes, points, tolerance);
6180 const auto &potential_owners_ranks = std::get<0>(potential_owners);
6181 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
6182 const auto &potential_owners_indices = std::get<2>(potential_owners);
6186 const auto translate = [&](
const unsigned int other_rank) {
6187 const auto ptr = std::find(potential_owners_ranks.begin(),
6188 potential_owners_ranks.end(),
6193 const auto other_rank_index =
6194 std::distance(potential_owners_ranks.begin(), ptr);
6196 return other_rank_index;
6200 (marked_vertices.size() == 0) ||
6203 "The marked_vertices vector has to be either empty or its size has "
6204 "to equal the number of vertices of the triangulation."));
6206 using RequestType = std::vector<std::pair<unsigned int, Point<spacedim>>>;
6207 using AnswerType = std::vector<unsigned int>;
6213 const bool has_relevant_vertices =
6214 (marked_vertices.size() == 0) ||
6215 (std::find(marked_vertices.begin(), marked_vertices.end(),
true) !=
6216 marked_vertices.end());
6218 const auto create_request = [&](
const unsigned int other_rank) {
6219 const auto other_rank_index = translate(other_rank);
6221 RequestType request;
6222 request.reserve(potential_owners_ptrs[other_rank_index + 1] -
6223 potential_owners_ptrs[other_rank_index]);
6225 for (
unsigned int i = potential_owners_ptrs[other_rank_index];
6226 i < potential_owners_ptrs[other_rank_index + 1];
6228 request.emplace_back(potential_owners_indices[i],
6229 points[potential_owners_indices[i]]);
6234 const auto answer_request =
6235 [&](
const unsigned int &other_rank,
6236 const RequestType & request) -> AnswerType {
6237 AnswerType answer(request.size(), 0);
6239 if (has_relevant_vertices)
6243 for (
unsigned int i = 0; i < request.size(); ++i)
6245 const auto &index_and_point = request[i];
6247 const auto cells_and_reference_positions =
6250 index_and_point.second,
6254 enforce_unique_mapping);
6259 for (
const auto &cell_and_reference_position :
6260 cells_and_reference_positions)
6262 const auto cell = cell_and_reference_position.first;
6263 auto reference_position =
6264 cell_and_reference_position.second;
6268 if (cell->reference_cell().is_hyper_cube())
6269 reference_position =
6271 reference_position);
6273 send_components.emplace_back(
6274 std::pair<int, int>(cell->level(), cell->index()),
6276 index_and_point.first,
6278 index_and_point.second,
6282 answer[i] = cells_and_reference_positions.size();
6286 if (perform_handshake)
6292 const auto process_answer = [&](
const unsigned int other_rank,
6293 const AnswerType & answer) {
6294 if (perform_handshake)
6296 const auto other_rank_index = translate(other_rank);
6298 for (
unsigned int i = 0; i < answer.size(); ++i)
6299 for (
unsigned int j = 0; j < answer[i]; ++j)
6300 recv_components.emplace_back(
6302 potential_owners_indices
6303 [i + potential_owners_ptrs[other_rank_index]],
6308 Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
6309 potential_owners_ranks,
6323 template <
int dim,
int spacedim>
6324 std::map<unsigned int, Point<spacedim>>
6328 std::map<unsigned int, Point<spacedim>> result;
6331 if (!cell->is_artificial())
6334 for (
unsigned int i = 0; i < vs.size(); ++i)
6335 result[cell->vertex_index(i)] = vs[i];
6342 template <
int spacedim>