18 #include <deal.II/base/mpi.templates.h>
63 #include <boost/random/mersenne_twister.hpp>
64 #include <boost/random/uniform_real_distribution.hpp>
75 #include <unordered_map>
82 template <
int dim,
int spacedim>
92 #if defined(DEAL_II_WITH_P4EST) && defined(DEBUG)
107 std::vector<bool> boundary_vertices(
vertices.size(),
false);
113 for (; cell != endc; ++cell)
114 for (
const unsigned int face : cell->face_indices())
115 if (cell->face(face)->at_boundary())
116 for (
unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
117 boundary_vertices[cell->face(face)->vertex_index(i)] =
true;
121 double max_distance_sqr = 0;
122 std::vector<bool>::const_iterator pi = boundary_vertices.begin();
123 const unsigned int N = boundary_vertices.size();
124 for (
unsigned int i = 0; i <
N; ++i, ++pi)
126 std::vector<bool>::const_iterator pj = pi + 1;
127 for (
unsigned int j = i + 1; j <
N; ++j, ++pj)
128 if ((*pi ==
true) && (*pj ==
true) &&
133 return std::sqrt(max_distance_sqr);
138 template <
int dim,
int spacedim>
144 unsigned int mapping_degree = 1;
146 mapping_degree = p->get_degree();
147 else if (
const auto *p =
149 mapping_degree = p->get_degree();
152 const QGauss<dim> quadrature_formula(mapping_degree + 1);
153 const unsigned int n_q_points = quadrature_formula.
size();
168 double local_volume = 0;
171 for (; cell != endc; ++cell)
172 if (cell->is_locally_owned())
175 for (
unsigned int q = 0; q < n_q_points; ++q)
176 local_volume += fe_values.
JxW(q);
179 double global_volume = 0;
181 #ifdef DEAL_II_WITH_MPI
189 global_volume = local_volume;
191 return global_volume;
213 struct TransformR2UAffine
228 [1] = {{-1.000000}, {1.000000}};
232 {1.000000, 0.000000};
243 [2] = {{-0.500000, -0.500000},
244 {0.500000, -0.500000},
245 {-0.500000, 0.500000},
246 {0.500000, 0.500000}};
256 {0.750000, 0.250000, 0.250000, -0.250000};
262 {-0.250000, -0.250000, -0.250000},
263 {0.250000, -0.250000, -0.250000},
264 {-0.250000, 0.250000, -0.250000},
265 {0.250000, 0.250000, -0.250000},
266 {-0.250000, -0.250000, 0.250000},
267 {0.250000, -0.250000, 0.250000},
268 {-0.250000, 0.250000, 0.250000},
269 {0.250000, 0.250000, 0.250000}
288 template <
int dim,
int spacedim>
297 for (
unsigned int d = 0;
d < spacedim; ++
d)
298 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
299 for (
unsigned int e = 0;
e < dim; ++
e)
304 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
307 return std::make_pair(
A,
b);
324 for (
const auto &cell :
triangulation.active_cell_iterators())
326 if (cell->is_locally_owned())
328 double aspect_ratio_cell = 0.0;
333 for (
unsigned int q = 0; q < quadrature.
size(); ++q)
344 aspect_ratio_cell = std::numeric_limits<double>::infinity();
349 for (
unsigned int i = 0; i < dim; ++i)
350 for (
unsigned int j = 0; j < dim; ++j)
351 J(i, j) = jacobian[i][j];
357 double const ar = max_sv / min_sv;
364 aspect_ratio_cell =
std::max(aspect_ratio_cell, ar);
369 aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
373 return aspect_ratio_vector;
394 template <
int dim,
int spacedim>
400 const auto predicate = [](
const iterator &) {
return true; };
403 tria, std::function<
bool(
const iterator &)>(predicate));
429 template <
int structdim>
453 "Two CellData objects with equal vertices must "
454 "have the same material/boundary ids and manifold "
479 template <
class FaceIteratorType>
483 CellData<dim - 1> face_cell_data(face->n_vertices());
484 for (
unsigned int vertex_n = 0; vertex_n < face->n_vertices();
486 face_cell_data.vertices[vertex_n] = face->vertex_index(vertex_n);
487 face_cell_data.boundary_id = face->boundary_id();
488 face_cell_data.manifold_id = face->manifold_id();
490 face_data.insert(std::move(face_cell_data));
518 template <
class FaceIteratorType>
533 template <
int dim,
int spacedim>
535 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>,
SubCellData>
539 ExcMessage(
"The input triangulation must be non-empty."));
541 std::vector<Point<spacedim>>
vertices;
542 std::vector<CellData<dim>> cells;
544 unsigned int max_level_0_vertex_n = 0;
546 for (
const unsigned int cell_vertex_n : cell->vertex_indices())
547 max_level_0_vertex_n =
548 std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
549 vertices.resize(max_level_0_vertex_n + 1);
559 for (
const unsigned int cell_vertex_n : cell->vertex_indices())
563 vertices[cell->vertex_index(cell_vertex_n)] =
564 cell->vertex(cell_vertex_n);
566 cell->vertex_index(cell_vertex_n);
570 cells.push_back(cell_data);
575 for (
const unsigned int face_n : cell->face_indices())
578 const auto face = cell->face(face_n);
587 for (
unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
589 const auto line = cell->line(line_n);
595 for (
unsigned int vertex_n : line->vertex_indices())
597 line->vertex_index(vertex_n);
600 line_data.insert(std::move(line_cell_data));
609 std::vector<bool> used_vertices(
vertices.size());
611 for (
const auto v : cell_data.vertices)
612 used_vertices[v] =
true;
613 Assert(std::find(used_vertices.begin(), used_vertices.end(),
false) ==
615 ExcMessage(
"The level zero vertices should form a contiguous "
623 for (
const CellData<1> &face_line_data : line_data)
626 return std::tuple<std::vector<Point<spacedim>>,
627 std::vector<CellData<dim>>,
630 std::move(subcell_data));
635 template <
int dim,
int spacedim>
644 "Invalid SubCellData supplied according to ::check_consistency(). "
645 "This is caused by data containing objects for the wrong dimension."));
648 std::vector<bool> vertex_used(
vertices.size(),
false);
649 for (
unsigned int c = 0; c < cells.size(); ++c)
650 for (
unsigned int v = 0; v < cells[c].vertices.size(); ++v)
653 ExcMessage(
"Invalid vertex index encountered! cells[" +
657 " is invalid, because only " +
659 " vertices were supplied."));
660 vertex_used[cells[c].vertices[v]] =
true;
667 std::vector<unsigned int> new_vertex_numbers(
vertices.size(),
669 unsigned int next_free_number = 0;
670 for (
unsigned int i = 0; i <
vertices.size(); ++i)
671 if (vertex_used[i] ==
true)
673 new_vertex_numbers[i] = next_free_number;
678 for (
unsigned int c = 0; c < cells.size(); ++c)
680 v = new_vertex_numbers[v];
685 for (
unsigned int v = 0;
690 new_vertex_numbers.size(),
692 "Invalid vertex index in subcelldata.boundary_lines. "
693 "subcelldata.boundary_lines[" +
698 " is invalid, because only " +
700 " vertices were supplied."));
707 for (
unsigned int v = 0;
712 new_vertex_numbers.size(),
714 "Invalid vertex index in subcelldata.boundary_quads. "
715 "subcelldata.boundary_quads[" +
720 " is invalid, because only " +
722 " vertices were supplied."));
730 std::vector<Point<spacedim>> tmp;
731 tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(),
true));
732 for (
unsigned int v = 0; v <
vertices.size(); ++v)
733 if (vertex_used[v] ==
true)
740 template <
int dim,
int spacedim>
745 std::vector<unsigned int> & considered_vertices,
752 std::vector<unsigned int> new_vertex_numbers(
vertices.size());
753 std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
756 if (considered_vertices.size() == 0)
757 considered_vertices = new_vertex_numbers;
772 unsigned int longest_coordinate_direction = 0;
773 double longest_coordinate_length = bbox.
side_length(0);
774 for (
unsigned int d = 1;
d < spacedim; ++
d)
777 if (longest_coordinate_length < coordinate_length)
779 longest_coordinate_length = coordinate_length;
780 longest_coordinate_direction =
d;
786 std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
787 sorted_vertices.reserve(
vertices.size());
788 for (
const unsigned int vertex_n : considered_vertices)
791 sorted_vertices.emplace_back(vertex_n,
vertices[vertex_n]);
793 std::sort(sorted_vertices.begin(),
794 sorted_vertices.end(),
797 return a.second[longest_coordinate_direction] <
798 b.second[longest_coordinate_direction];
803 for (
unsigned int d = 0;
d < spacedim; ++
d)
804 if (std::abs(a[
d] -
b[
d]) > tol)
811 auto range_start = sorted_vertices.begin();
812 while (range_start != sorted_vertices.end())
814 auto range_end = range_start + 1;
815 while (range_end != sorted_vertices.end() &&
816 std::abs(range_end->second[longest_coordinate_direction] -
817 range_start->second[longest_coordinate_direction]) <
823 std::sort(range_start,
827 return a.first <
b.first;
836 for (
auto reference = range_start; reference != range_end; ++reference)
839 for (
auto it = reference + 1; it != range_end; ++it)
841 if (within_tolerance(reference->second, it->second))
843 new_vertex_numbers[it->first] = reference->first;
849 range_start = range_end;
856 for (
auto &cell : cells)
857 for (
auto &vertex_index : cell.vertices)
858 vertex_index = new_vertex_numbers[vertex_index];
860 for (
auto &vertex_index : quad.vertices)
861 vertex_index = new_vertex_numbers[vertex_index];
863 for (
auto &vertex_index : line.vertices)
864 vertex_index = new_vertex_numbers[vertex_index];
884 for (
unsigned int i = 0; i <
vertices.size(); ++i)
885 map_point_to_local_vertex_index[
vertices[i]] = i;
888 if (map_point_to_local_vertex_index.size() ==
vertices.size())
892 vertices.resize(map_point_to_local_vertex_index.size());
895 for (
const auto &p : map_point_to_local_vertex_index)
902 template <
int dim,
int spacedim>
913 if (dim == 2 && spacedim == 3)
916 std::size_t n_negative_cells = 0;
917 std::size_t cell_no = 0;
918 for (
auto &cell : cells)
934 std::swap(cell.vertices[1], cell.vertices[2]);
939 std::swap(cell.vertices[0], cell.vertices[2]);
940 std::swap(cell.vertices[1], cell.vertices[3]);
941 std::swap(cell.vertices[4], cell.vertices[6]);
942 std::swap(cell.vertices[5], cell.vertices[7]);
950 std::swap(cell.vertices[cell.vertices.size() - 2],
951 cell.vertices[cell.vertices.size() - 1]);
956 std::swap(cell.vertices[0], cell.vertices[3]);
957 std::swap(cell.vertices[1], cell.vertices[4]);
958 std::swap(cell.vertices[2], cell.vertices[5]);
964 std::swap(cell.vertices[2], cell.vertices[3]);
978 return n_negative_cells;
982 template <
int dim,
int spacedim>
988 const std::size_t n_negative_cells =
997 AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
1000 "This function assumes that either all cells have positive "
1001 "volume, or that all cells have been specified in an "
1002 "inverted vertex order so that their volume is negative. "
1003 "(In the latter case, this class automatically inverts "
1004 "every cell.) However, the mesh you have specified "
1005 "appears to have both cells with positive and cells with "
1006 "negative volume. You need to check your mesh which "
1007 "cells these are and how they got there.\n"
1008 "As a hint, of the total ") +
1011 " appear to have a negative volume."));
1029 CheapEdge(
const unsigned int v0,
const unsigned int v1)
1041 return ((
v0 <
e.v0) || ((
v0 ==
e.v0) && (
v1 <
e.v1)));
1064 std::set<CheapEdge> edges;
1066 for (
typename std::vector<
CellData<dim>>::const_iterator c =
1076 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1078 const CheapEdge reverse_edge(
1081 if (edges.find(reverse_edge) != edges.end())
1086 const CheapEdge correct_edge(
1089 edges.insert(correct_edge);
1105 struct ParallelEdges
1119 static const unsigned int
1184 class AdjacentCells;
1192 class AdjacentCells<2>
1199 using const_iterator =
const AdjacentCell *;
1210 push_back(
const AdjacentCell &adjacent_cell)
1274 class AdjacentCells<3> :
public std::vector<AdjacentCell>
1296 Edge(
const CellData<dim> &cell,
const unsigned int edge_number)
1347 enum OrientationStatus
1377 Cell(
const CellData<dim> &c,
const std::vector<Edge<dim>> &edge_list)
1384 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1386 const Edge<dim>
e(c,
l);
1422 class EdgeDeltaSet<2>
1428 using const_iterator =
const unsigned int *;
1454 insert(
const unsigned int edge_index)
1515 class EdgeDeltaSet<3> :
public std::set<unsigned int>
1525 std::vector<Edge<dim>>
1532 std::vector<Edge<dim>> edge_list;
1534 for (
unsigned int i = 0; i < cells.size(); ++i)
1535 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1536 edge_list.emplace_back(cells[i],
l);
1539 std::sort(edge_list.begin(), edge_list.end());
1540 edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
1553 std::vector<Cell<dim>>
1554 build_cells_and_connect_edges(
const std::vector<
CellData<dim>> &cells,
1555 std::vector<Edge<dim>> & edges)
1557 std::vector<Cell<dim>> cell_list;
1558 cell_list.reserve(cells.size());
1559 for (
unsigned int i = 0; i < cells.size(); ++i)
1563 cell_list.emplace_back(cells[i], edges);
1567 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1568 edges[cell_list.back().edge_indices[
l]].adjacent_cells.push_back(
1569 AdjacentCell(i,
l));
1584 get_next_unoriented_cell(
const std::vector<Cell<dim>> &cells,
1585 const std::vector<Edge<dim>> &edges,
1586 const unsigned int current_cell)
1588 for (
unsigned int c = current_cell; c < cells.size(); ++c)
1589 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1591 Edge<dim>::not_oriented)
1606 orient_one_set_of_parallel_edges(
const std::vector<Cell<dim>> &cells,
1607 std::vector<Edge<dim>> & edges,
1608 const unsigned int cell,
1609 const unsigned int local_edge)
1636 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1637 (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
1653 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1654 (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
1663 EdgeDeltaSet<dim> Delta_k;
1664 EdgeDeltaSet<dim> Delta_k_minus_1;
1665 Delta_k_minus_1.insert(cells[cell].
edge_indices[local_edge]);
1667 while (Delta_k_minus_1.begin() !=
1668 Delta_k_minus_1.end())
1672 for (
typename EdgeDeltaSet<dim>::const_iterator delta =
1673 Delta_k_minus_1.begin();
1674 delta != Delta_k_minus_1.end();
1678 Edge<dim>::not_oriented,
1682 for (
typename AdjacentCells<dim>::const_iterator adjacent_cell =
1684 adjacent_cell != edges[*delta].adjacent_cells.end();
1687 const unsigned int K = adjacent_cell->cell_index;
1688 const unsigned int delta_is_edge_in_K =
1689 adjacent_cell->edge_within_cell;
1693 const unsigned int first_edge_vertex =
1694 (edges[*delta].orientation_status == Edge<dim>::forward ?
1695 edges[*delta].vertex_indices[0] :
1696 edges[*delta].vertex_indices[1]);
1697 const unsigned int first_edge_vertex_in_K =
1700 delta_is_edge_in_K, 0)];
1702 first_edge_vertex == first_edge_vertex_in_K ||
1703 first_edge_vertex ==
1705 dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
1710 for (
unsigned int o_e = 0;
1717 const unsigned int opposite_edge =
1718 cells[
K].edge_indices[ParallelEdges<
1720 const unsigned int first_opposite_edge_vertex =
1721 cells[
K].vertex_indices
1725 (first_edge_vertex == first_edge_vertex_in_K ? 0 :
1732 const typename Edge<dim>::OrientationStatus
1733 opposite_edge_orientation =
1734 (edges[opposite_edge].vertex_indices[0] ==
1735 first_opposite_edge_vertex ?
1736 Edge<dim>::forward :
1737 Edge<dim>::backward);
1742 Edge<dim>::not_oriented)
1746 edges[opposite_edge].orientation_status =
1747 opposite_edge_orientation;
1748 Delta_k.insert(opposite_edge);
1764 opposite_edge_orientation,
1770 opposite_edge_orientation)
1783 Delta_k_minus_1 = Delta_k;
1797 rotate_cell(
const std::vector<Cell<dim>> &cell_list,
1798 const std::vector<Edge<dim>> &edge_list,
1805 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++
e)
1812 starting_vertex_of_edge[
e] =
1816 starting_vertex_of_edge[
e] =
1832 if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
1833 (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
1834 origin_vertex_of_cell = starting_vertex_of_edge[0];
1835 else if ((starting_vertex_of_edge[1] ==
1836 starting_vertex_of_edge[2]) ||
1837 (starting_vertex_of_edge[1] ==
1838 starting_vertex_of_edge[3]))
1839 origin_vertex_of_cell = starting_vertex_of_edge[1];
1851 for (origin_vertex_of_cell = 0;
1852 origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
1853 ++origin_vertex_of_cell)
1854 if (std::count(starting_vertex_of_edge,
1855 starting_vertex_of_edge +
1860 Assert(origin_vertex_of_cell <
1885 const unsigned int tmp = raw_cells[
cell_index].vertices[0];
1908 static const unsigned int cube_permutations[8][8] = {
1909 {0, 1, 2, 3, 4, 5, 6, 7},
1910 {1, 5, 3, 7, 0, 4, 2, 6},
1911 {2, 6, 0, 4, 3, 7, 1, 5},
1912 {3, 2, 1, 0, 7, 6, 5, 4},
1913 {4, 0, 6, 2, 5, 1, 7, 3},
1914 {5, 4, 7, 6, 1, 0, 3, 2},
1915 {6, 7, 4, 5, 2, 3, 0, 1},
1916 {7, 3, 5, 1, 6, 2, 4, 0}};
1921 temp_vertex_indices[v] =
1923 .vertices[cube_permutations[origin_vertex_of_cell][v]];
1925 raw_cells[
cell_index].vertices[v] = temp_vertex_indices[v];
1949 std::vector<Edge<dim>> edge_list = build_edges(cells);
1950 std::vector<Cell<dim>> cell_list =
1951 build_cells_and_connect_edges(cells, edge_list);
1955 unsigned int next_cell_with_unoriented_edge = 0;
1956 while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
1957 cell_list, edge_list, next_cell_with_unoriented_edge)) !=
1971 for (
unsigned int l = 0;
l < dim; ++
l)
1972 if (edge_list[cell_list[next_cell_with_unoriented_edge]
1975 orient_one_set_of_parallel_edges(
1978 next_cell_with_unoriented_edge,
1983 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1984 Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
1993 for (
unsigned int c = 0; c < cells.size(); ++c)
1994 rotate_cell(cell_list, edge_list, c, cells);
2009 Assert(cells.size() != 0,
2011 "List of elements to orient must have at least one cell"));
2019 if (!is_consistent(cells))
2038 template <
int spacedim>
2096 template <
int spacedim>
2115 template <
int dim,
int spacedim>
2125 template <
int dim,
int spacedim>
2134 "GridTools::rotate() is only available for spacedim = 2."));
2169 const unsigned int axis,
2181 template <
int dim,
int spacedim>
2203 const unsigned int n_dofs = S.
n();
2215 const auto constrained_rhs =
2217 solver.
solve(SF, u, constrained_rhs, prec);
2230 const bool solve_for_absolute_positions)
2256 const auto reference_cell = ReferenceCells::get_hypercube<dim>();
2265 std::array<AffineConstraints<double>, dim> constraints;
2266 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
2274 for (
const unsigned int vertex_no : cell->vertex_indices())
2276 const unsigned int vertex_index = cell->vertex_index(vertex_no);
2277 const Point<dim> & vertex_point = cell->vertex(vertex_no);
2279 const typename std::map<unsigned int, Point<dim>>::const_iterator
2280 map_iter = new_points.find(vertex_index);
2282 if (map_iter != map_end)
2283 for (
unsigned int i = 0; i < dim; ++i)
2285 constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
2286 constraints[i].set_inhomogeneity(
2287 cell->vertex_dof_index(vertex_no, 0),
2288 (solve_for_absolute_positions ?
2289 map_iter->second(i) :
2290 map_iter->second(i) - vertex_point[i]));
2295 for (
unsigned int i = 0; i < dim; ++i)
2296 constraints[i].close();
2300 for (
unsigned int i = 0; i < dim; ++i)
2305 for (
unsigned int i = 0; i < dim; ++i)
2312 std::vector<bool> vertex_touched(
triangulation.n_vertices(),
false);
2314 for (
const unsigned int vertex_no : cell->vertex_indices())
2315 if (vertex_touched[cell->vertex_index(vertex_no)] ==
false)
2320 cell->vertex_dof_index(vertex_no, 0);
2321 for (
unsigned int i = 0; i < dim; ++i)
2322 if (solve_for_absolute_positions)
2323 v(i) = us[i](dof_index);
2325 v(i) += us[i](dof_index);
2327 vertex_touched[cell->vertex_index(vertex_no)] =
true;
2331 template <
int dim,
int spacedim>
2332 std::map<unsigned int, Point<spacedim>>
2335 std::map<unsigned int, Point<spacedim>> vertex_map;
2339 for (; cell != endc; ++cell)
2341 for (
unsigned int i : cell->face_indices())
2345 if (face->at_boundary())
2347 for (
unsigned j = 0; j < face->
n_vertices(); ++j)
2350 const unsigned int vertex_index = face->vertex_index(j);
2351 vertex_map[vertex_index] = vertex;
2363 template <
int dim,
int spacedim>
2367 const bool keep_boundary,
2368 const unsigned int seed)
2385 double almost_infinite_length = 0;
2390 almost_infinite_length += cell->diameter();
2392 std::vector<double> minimal_length(
triangulation.n_vertices(),
2393 almost_infinite_length);
2396 std::vector<bool> at_boundary(keep_boundary ?
triangulation.n_vertices() :
2401 const bool is_parallel_shared =
2404 for (
const auto &cell :
triangulation.active_cell_iterators())
2405 if (is_parallel_shared || cell->is_locally_owned())
2409 for (
unsigned int i = 0; i < cell->n_lines(); ++i)
2412 line = cell->line(i);
2414 if (keep_boundary && line->at_boundary())
2416 at_boundary[line->vertex_index(0)] =
true;
2417 at_boundary[line->vertex_index(1)] =
true;
2420 minimal_length[line->vertex_index(0)] =
2422 minimal_length[line->vertex_index(0)]);
2423 minimal_length[line->vertex_index(1)] =
2425 minimal_length[line->vertex_index(1)]);
2431 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
2432 if (cell->at_boundary(vertex) ==
true)
2433 at_boundary[cell->vertex_index(vertex)] =
true;
2435 minimal_length[cell->vertex_index(0)] =
2437 minimal_length[cell->vertex_index(0)]);
2438 minimal_length[cell->vertex_index(1)] =
2440 minimal_length[cell->vertex_index(1)]);
2445 boost::random::mt19937 rng(seed);
2446 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
2450 if (
auto distributed_triangulation =
2454 const std::vector<bool> locally_owned_vertices =
2456 std::vector<bool> vertex_moved(
triangulation.n_vertices(),
false);
2459 for (
const auto &cell :
triangulation.active_cell_iterators())
2460 if (cell->is_locally_owned())
2462 for (
const unsigned int vertex_no : cell->vertex_indices())
2464 const unsigned global_vertex_no =
2465 cell->vertex_index(vertex_no);
2470 if ((keep_boundary && at_boundary[global_vertex_no]) ||
2471 vertex_moved[global_vertex_no] ||
2472 !locally_owned_vertices[global_vertex_no])
2477 for (
unsigned int d = 0;
d < spacedim; ++
d)
2478 shift_vector(
d) = uniform_distribution(rng);
2480 shift_vector *= factor * minimal_length[global_vertex_no] /
2481 std::sqrt(shift_vector.
square());
2484 cell->vertex(vertex_no) += shift_vector;
2485 vertex_moved[global_vertex_no] =
true;
2489 distributed_triangulation->communicate_locally_moved_vertices(
2490 locally_owned_vertices);
2500 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
2501 const std::vector<Point<spacedim>> &old_vertex_locations =
2504 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2508 if (keep_boundary && at_boundary[vertex])
2509 new_vertex_locations[vertex] = old_vertex_locations[vertex];
2514 for (
unsigned int d = 0;
d < spacedim; ++
d)
2515 shift_vector(
d) = uniform_distribution(rng);
2517 shift_vector *= factor * minimal_length[vertex] /
2518 std::sqrt(shift_vector.
square());
2521 new_vertex_locations[vertex] =
2522 old_vertex_locations[vertex] + shift_vector;
2527 for (
const auto &cell :
triangulation.active_cell_iterators())
2528 for (
const unsigned int vertex_no : cell->vertex_indices())
2529 cell->vertex(vertex_no) =
2530 new_vertex_locations[cell->vertex_index(vertex_no)];
2546 for (; cell != endc; ++cell)
2547 if (!cell->is_artificial())
2548 for (
const unsigned int face : cell->face_indices())
2549 if (cell->face(face)->has_children() &&
2550 !cell->face(face)->at_boundary())
2554 cell->face(face)->child(0)->vertex(1) =
2555 (cell->face(face)->vertex(0) +
2556 cell->face(face)->vertex(1)) /
2560 cell->face(face)->child(0)->vertex(1) =
2561 .5 * (cell->face(face)->vertex(0) +
2562 cell->face(face)->vertex(1));
2563 cell->face(face)->child(0)->vertex(2) =
2564 .5 * (cell->face(face)->vertex(0) +
2565 cell->face(face)->vertex(2));
2566 cell->face(face)->child(1)->vertex(3) =
2567 .5 * (cell->face(face)->vertex(1) +
2568 cell->face(face)->vertex(3));
2569 cell->face(face)->child(2)->vertex(3) =
2570 .5 * (cell->face(face)->vertex(2) +
2571 cell->face(face)->vertex(3));
2574 cell->face(face)->child(0)->vertex(3) =
2575 .25 * (cell->face(face)->vertex(0) +
2576 cell->face(face)->vertex(1) +
2577 cell->face(face)->vertex(2) +
2578 cell->face(face)->vertex(3));
2586 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2590 const std::vector<bool> & marked_vertices)
2599 marked_vertices.size() == 0,
2601 marked_vertices.size()));
2608 marked_vertices.size() == 0 ||
2609 std::equal(marked_vertices.begin(),
2610 marked_vertices.end(),
2612 [](
bool p,
bool q) { return !p || q; }),
2614 "marked_vertices should be a subset of used vertices in the triangulation "
2615 "but marked_vertices contains one or more vertices that are not used vertices!"));
2619 const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
2625 std::vector<bool>::const_iterator
first =
2626 std::find(vertices_to_use.begin(), vertices_to_use.end(),
true);
2631 unsigned int best_vertex = std::distance(vertices_to_use.begin(),
first);
2632 double best_dist = (p -
vertices[best_vertex]).norm_square();
2636 for (
unsigned int j = best_vertex + 1; j <
vertices.size(); ++j)
2637 if (vertices_to_use[j])
2639 const double dist = (p -
vertices[j]).norm_square();
2640 if (dist < best_dist)
2652 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2655 const MeshType<dim, spacedim> &mesh,
2657 const std::vector<bool> & marked_vertices)
2670 marked_vertices.size() == 0,
2672 marked_vertices.size()));
2680 marked_vertices.size() == 0 ||
2681 std::equal(marked_vertices.begin(),
2682 marked_vertices.end(),
2684 [](
bool p,
bool q) { return !p || q; }),
2686 "marked_vertices should be a subset of used vertices in the triangulation "
2687 "but marked_vertices contains one or more vertices that are not used vertices!"));
2690 if (marked_vertices.size() != 0)
2693 if (marked_vertices[it->first] ==
false)
2708 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2710 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
2713 typename ::internal::
2714 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2717 const unsigned int vertex)
2723 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
2730 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2766 for (
const auto &cell : mesh.active_cell_iterators())
2768 for (
const unsigned int v : cell->vertex_indices())
2769 if (cell->vertex_index(v) == vertex)
2781 for (
const auto face :
2782 cell->reference_cell().faces_for_given_vertex(v))
2783 if (!cell->at_boundary(face) &&
2784 cell->neighbor(face)->is_active())
2806 for (
unsigned int e = 0;
e < cell->n_lines(); ++
e)
2807 if (cell->line(
e)->has_children())
2811 if (cell->line(
e)->child(0)->vertex_index(1) == vertex)
2835 typename ::internal::
2836 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
2842 template <
int dim,
int spacedim>
2843 std::vector<std::vector<Tensor<1, spacedim>>>
2851 const unsigned int n_vertices = vertex_to_cells.size();
2856 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
2858 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2861 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
2862 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
2864 typename std::set<typename Triangulation<dim, spacedim>::
2865 active_cell_iterator>::iterator it =
2866 vertex_to_cells[vertex].begin();
2867 for (
unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
2869 vertex_to_cell_centers[vertex][cell] =
2870 (*it)->center() -
vertices[vertex];
2871 vertex_to_cell_centers[vertex][cell] /=
2872 vertex_to_cell_centers[vertex][cell].
norm();
2875 return vertex_to_cell_centers;
2881 template <
int spacedim>
2884 const unsigned int a,
2885 const unsigned int b,
2889 const double scalar_product_a = center_directions[a] * point_direction;
2890 const double scalar_product_b = center_directions[
b] * point_direction;
2895 return (scalar_product_a > scalar_product_b);
2899 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2901 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
2903 std::pair<typename ::internal::
2904 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
2909 const MeshType<dim, spacedim> &mesh,
2912 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
2915 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
2916 const std::vector<bool> & marked_vertices,
2918 const double tolerance,
2922 *relevant_cell_bounding_boxes_rtree)
2924 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2927 cell_and_position.first = mesh.end();
2931 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2933 cell_and_position_approx;
2935 if (relevant_cell_bounding_boxes_rtree !=
nullptr &&
2936 !relevant_cell_bounding_boxes_rtree->empty())
2942 for (
unsigned int d = 0;
d < spacedim; ++
d)
2944 p1[
d] = p1[
d] - tolerance;
2945 p2[
d] = p2[
d] + tolerance;
2950 if (relevant_cell_bounding_boxes_rtree->qbegin(
2951 boost::geometry::index::intersects(bb)) ==
2952 relevant_cell_bounding_boxes_rtree->qend())
2953 return cell_and_position;
2956 bool found_cell =
false;
2957 bool approx_cell =
false;
2959 unsigned int closest_vertex_index = 0;
2964 if (marked_vertices.size() && !used_vertices_rtree.empty())
2967 std::find(marked_vertices.begin(), marked_vertices.end(),
true);
2968 Assert(itr != marked_vertices.end(),
2970 closest_vertex_index = std::distance(marked_vertices.begin(), itr);
2974 auto current_cell = cell_hint;
2977 const auto cell_marked = [&mesh, &marked_vertices](
const auto &cell) {
2978 if (marked_vertices.size() == 0)
2981 if (cell != mesh.active_cell_iterators().end())
2982 for (
unsigned int i = 0; i < cell->n_vertices(); ++i)
2983 if (marked_vertices[cell->vertex_index(i)])
2990 const auto any_cell_marked = [&cell_marked](
const auto &cells) {
2991 return std::any_of(cells.begin(),
2993 [&cell_marked](
const auto &cell) {
2994 return cell_marked(cell);
2998 while (found_cell ==
false)
3003 cell_marked(cell_hint))
3005 const auto cell_vertices = mapping.
get_vertices(current_cell);
3006 const unsigned int closest_vertex =
3007 find_closest_vertex_of_cell<dim, spacedim>(current_cell,
3010 vertex_to_point = p - cell_vertices[closest_vertex];
3011 closest_vertex_index = current_cell->vertex_index(closest_vertex);
3020 #if defined(DEAL_II_WITH_BOOST_BUNDLED) || \
3021 !(defined(__clang_major__) && __clang_major__ >= 16) || \
3022 BOOST_VERSION >= 108100
3023 if (!used_vertices_rtree.empty())
3026 using ValueType = std::pair<Point<spacedim>,
unsigned int>;
3027 std::function<
bool(
const ValueType &)> marked;
3028 if (marked_vertices.size() == mesh.n_vertices())
3029 marked = [&marked_vertices](
const ValueType &value) ->
bool {
3030 return marked_vertices[value.second];
3033 marked = [](
const ValueType &) ->
bool {
return true; };
3035 std::vector<std::pair<Point<spacedim>,
unsigned int>> res;
3036 used_vertices_rtree.query(
3037 boost::geometry::index::nearest(p, 1) &&
3038 boost::geometry::index::satisfies(marked),
3039 std::back_inserter(res));
3044 ::
ExcMessage(
"There can not be multiple results"));
3047 if (any_cell_marked(vertex_to_cells[res[0].
second]))
3048 closest_vertex_index = res[0].second;
3054 mapping, mesh, p, marked_vertices);
3056 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
3062 Assert(any_cell_marked(vertex_to_cells[closest_vertex_index]),
3067 const double vertex_point_norm = vertex_to_point.
norm();
3068 if (vertex_point_norm > 0)
3069 vertex_to_point /= vertex_point_norm;
3071 const unsigned int n_neighbor_cells =
3072 vertex_to_cells[closest_vertex_index].size();
3075 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
3077 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
3078 neighbor_permutation[i] = i;
3080 auto comp = [&](
const unsigned int a,
const unsigned int b) ->
bool {
3081 return internal::compare_point_association<spacedim>(
3085 vertex_to_cell_centers[closest_vertex_index]);
3088 std::sort(neighbor_permutation.begin(),
3089 neighbor_permutation.end(),
3094 double best_distance = tolerance;
3098 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
3102 auto cell = vertex_to_cells[closest_vertex_index].begin();
3105 if (!(*cell)->is_artificial())
3112 cell_and_position.first = *cell;
3113 cell_and_position.second = p_unit;
3115 approx_cell =
false;
3124 if (dist < best_distance)
3126 best_distance = dist;
3127 cell_and_position_approx.first = *cell;
3128 cell_and_position_approx.second = p_unit;
3137 if (found_cell ==
true)
3138 return cell_and_position;
3139 else if (approx_cell ==
true)
3140 return cell_and_position_approx;
3156 mapping, mesh, p, marked_vertices, tolerance);
3158 current_cell =
typename MeshType<dim, spacedim>::active_cell_iterator();
3160 return cell_and_position;
3165 template <
int dim,
int spacedim>
3174 unsigned int closest_vertex = 0;
3176 for (
unsigned int v = 1; v < cell->
n_vertices(); ++v)
3179 if (vertex_distance < minimum_distance)
3182 minimum_distance = vertex_distance;
3185 return closest_vertex;
3192 namespace BoundingBoxPredicate
3194 template <
class MeshType>
3195 std::tuple<BoundingBox<MeshType::space_dimension>,
bool>
3197 const typename MeshType::cell_iterator &parent_cell,
3198 const std::function<
3199 bool(
const typename MeshType::active_cell_iterator &)> &predicate)
3201 bool has_predicate =
3203 std::vector<typename MeshType::active_cell_iterator> active_cells;
3204 if (parent_cell->is_active())
3205 active_cells = {parent_cell};
3209 active_cells = get_active_child_cells<MeshType>(parent_cell);
3211 const unsigned int spacedim = MeshType::space_dimension;
3215 while (i < active_cells.size() && !predicate(active_cells[i]))
3219 if (active_cells.size() == 0 || i == active_cells.size())
3222 return std::make_tuple(bbox, has_predicate);
3229 for (; i < active_cells.size(); ++i)
3230 if (predicate(active_cells[i]))
3232 for (
unsigned int d = 0;
d < spacedim; ++
d)
3234 minp[
d] =
std::min(minp[
d], active_cells[i]->vertex(v)[
d]);
3235 maxp[
d] =
std::max(maxp[
d], active_cells[i]->vertex(v)[
d]);
3238 has_predicate =
true;
3240 return std::make_tuple(bbox, has_predicate);
3247 template <
class MeshType>
3248 std::vector<BoundingBox<MeshType::space_dimension>>
3250 const MeshType &mesh,
3251 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
3253 const unsigned int refinement_level,
3254 const bool allow_merge,
3255 const unsigned int max_boxes)
3262 refinement_level <= mesh.n_levels(),
3264 "Error: refinement level is higher then total levels in the triangulation!"));
3266 const unsigned int spacedim = MeshType::space_dimension;
3267 std::vector<BoundingBox<spacedim>> bounding_boxes;
3271 for (
unsigned int i = 0; i < refinement_level; ++i)
3272 for (
const typename MeshType::cell_iterator &cell :
3273 mesh.active_cell_iterators_on_level(i))
3275 bool has_predicate =
false;
3277 std::tie(bbox, has_predicate) =
3279 MeshType>(cell, predicate);
3281 bounding_boxes.push_back(bbox);
3285 for (
const typename MeshType::cell_iterator &cell :
3286 mesh.cell_iterators_on_level(refinement_level))
3288 bool has_predicate =
false;
3290 std::tie(bbox, has_predicate) =
3292 MeshType>(cell, predicate);
3294 bounding_boxes.push_back(bbox);
3299 return bounding_boxes;
3305 std::vector<unsigned int> merged_boxes_idx;
3306 bool found_neighbors =
true;
3311 while (found_neighbors)
3313 found_neighbors =
false;
3314 for (
unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
3316 if (std::find(merged_boxes_idx.begin(),
3317 merged_boxes_idx.end(),
3318 i) == merged_boxes_idx.end())
3319 for (
unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
3320 if (std::find(merged_boxes_idx.begin(),
3321 merged_boxes_idx.end(),
3322 j) == merged_boxes_idx.end() &&
3323 bounding_boxes[i].get_neighbor_type(
3324 bounding_boxes[j]) ==
3327 bounding_boxes[i].merge_with(bounding_boxes[j]);
3328 merged_boxes_idx.push_back(j);
3329 found_neighbors =
true;
3335 std::vector<BoundingBox<spacedim>> merged_b_boxes;
3336 for (
unsigned int i = 0; i < bounding_boxes.size(); ++i)
3337 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
3338 merged_boxes_idx.end())
3339 merged_b_boxes.push_back(bounding_boxes[i]);
3344 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
3346 std::vector<double> volumes;
3347 for (
unsigned int i = 0; i < merged_b_boxes.size(); ++i)
3348 volumes.push_back(merged_b_boxes[i].volume());
3350 while (merged_b_boxes.size() > max_boxes)
3352 unsigned int min_idx =
3353 std::min_element(volumes.begin(), volumes.end()) -
3355 volumes.erase(volumes.begin() + min_idx);
3357 bool not_removed =
true;
3358 for (
unsigned int i = 0;
3359 i < merged_b_boxes.size() && not_removed;
3364 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
3365 merged_b_boxes[min_idx]) ==
3367 merged_b_boxes[i].get_neighbor_type(
3368 merged_b_boxes[min_idx]) ==
3371 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
3372 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
3373 not_removed =
false;
3376 ExcMessage(
"Error: couldn't merge bounding boxes!"));
3379 Assert(merged_b_boxes.size() <= max_boxes,
3381 "Error: couldn't reach target number of bounding boxes!"));
3382 return merged_b_boxes;
3388 template <
int spacedim>
3390 std::tuple<std::vector<std::vector<unsigned int>>,
3391 std::map<unsigned int, unsigned int>,
3392 std::map<unsigned int, std::vector<unsigned int>>>
3400 unsigned int n_procs = global_bboxes.size();
3401 std::vector<std::vector<unsigned int>> point_owners(n_procs);
3402 std::map<unsigned int, unsigned int> map_owners_found;
3403 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3405 unsigned int n_points = points.size();
3406 for (
unsigned int pt = 0; pt < n_points; ++pt)
3409 std::vector<unsigned int> owners_found;
3411 for (
unsigned int rk = 0; rk < n_procs; ++rk)
3414 if (bbox.point_inside(points[pt]))
3416 point_owners[rk].emplace_back(pt);
3417 owners_found.emplace_back(rk);
3421 Assert(owners_found.size() > 0,
3422 ExcMessage(
"No owners found for the point " +
3424 if (owners_found.size() == 1)
3425 map_owners_found[pt] = owners_found[0];
3428 map_owners_guessed[pt] = owners_found;
3431 return std::make_tuple(std::move(point_owners),
3432 std::move(map_owners_found),
3433 std::move(map_owners_guessed));
3436 template <
int spacedim>
3438 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
3439 std::map<unsigned int, unsigned int>,
3440 std::map<unsigned int, std::vector<unsigned int>>>
3448 std::map<unsigned int, std::vector<unsigned int>> point_owners;
3449 std::map<unsigned int, unsigned int> map_owners_found;
3450 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3451 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> search_result;
3453 unsigned int n_points = points.size();
3454 for (
unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
3456 search_result.clear();
3459 covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
3460 std::back_inserter(search_result));
3463 std::set<unsigned int> owners_found;
3465 for (
const auto &rank_bbox : search_result)
3469 const bool pt_inserted = owners_found.insert(pt_n).second;
3471 point_owners[rank_bbox.second].emplace_back(pt_n);
3473 Assert(owners_found.size() > 0,
3474 ExcMessage(
"No owners found for the point " +
3476 if (owners_found.size() == 1)
3477 map_owners_found[pt_n] = *owners_found.begin();
3482 std::back_inserter(map_owners_guessed[pt_n]));
3485 return std::make_tuple(std::move(point_owners),
3486 std::move(map_owners_found),
3487 std::move(map_owners_guessed));
3491 template <
int dim,
int spacedim>
3493 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3497 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3502 for (; cell != endc; ++cell)
3503 for (
const unsigned int i : cell->vertex_indices())
3508 for (; cell != endc; ++cell)
3510 for (
unsigned int i : cell->face_indices())
3512 if ((cell->at_boundary(i) ==
false) &&
3513 (cell->neighbor(i)->is_active()))
3516 adjacent_cell = cell->neighbor(i);
3517 for (
unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
3526 for (
unsigned int i = 0; i < cell->
n_lines(); ++i)
3527 if (cell->line(i)->has_children())
3541 template <
int dim,
int spacedim>
3542 std::map<unsigned int, types::global_vertex_index>
3546 std::map<unsigned int, types::global_vertex_index>
3547 local_to_global_vertex_index;
3549 #ifndef DEAL_II_WITH_MPI
3556 ExcMessage(
"This function does not make any sense "
3557 "for parallel::distributed::Triangulation "
3558 "objects if you do not have MPI enabled."));
3562 using active_cell_iterator =
3564 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
3569 unsigned int max_cellid_size = 0;
3570 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
3572 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
3578 std::set<active_cell_iterator> missing_vert_cells;
3579 std::set<unsigned int> used_vertex_index;
3580 for (
const auto &cell :
triangulation.active_cell_iterators())
3582 if (cell->is_locally_owned())
3584 for (
const unsigned int i : cell->vertex_indices())
3587 for (
const auto &adjacent_cell :
3588 vertex_to_cell[cell->vertex_index(i)])
3589 lowest_subdomain_id =
std::min(lowest_subdomain_id,
3590 adjacent_cell->subdomain_id());
3593 if (lowest_subdomain_id == cell->subdomain_id())
3597 if (used_vertex_index.find(cell->vertex_index(i)) ==
3598 used_vertex_index.end())
3601 local_to_global_vertex_index[cell->vertex_index(i)] =
3606 for (
const auto &adjacent_cell :
3607 vertex_to_cell[cell->vertex_index(i)])
3608 if (adjacent_cell->subdomain_id() !=
3609 cell->subdomain_id())
3613 tmp(adjacent_cell->subdomain_id(),
3614 cell->vertex_index(i));
3615 if (vertices_added.find(tmp) ==
3616 vertices_added.end())
3618 vertices_to_send[adjacent_cell
3621 cell->vertex_index(i),
3622 cell->id().to_string());
3623 if (cell->id().to_string().size() >
3626 cell->id().to_string().size();
3627 vertices_added.insert(tmp);
3630 used_vertex_index.insert(cell->vertex_index(i));
3637 vertices_to_recv[lowest_subdomain_id].insert(
3638 cell->vertex_index(i));
3639 missing_vert_cells.insert(cell);
3646 if (cell->is_ghost())
3648 for (
unsigned int i : cell->face_indices())
3650 if (cell->at_boundary(i) ==
false)
3652 if (cell->neighbor(i)->is_active())
3655 spacedim>::active_cell_iterator
3656 adjacent_cell = cell->neighbor(i);
3657 if ((adjacent_cell->is_locally_owned()))
3660 adjacent_cell->subdomain_id();
3661 if (cell->subdomain_id() < adj_subdomain_id)
3662 for (
unsigned int j = 0;
3663 j < cell->face(i)->n_vertices();
3666 vertices_to_recv[cell->subdomain_id()].insert(
3667 cell->face(i)->vertex_index(j));
3668 missing_vert_cells.insert(cell);
3684 int ierr = MPI_Exscan(&next_index,
3692 for (
auto &global_index_it : local_to_global_vertex_index)
3693 global_index_it.second +=
shift;
3707 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
3708 vertices_to_send.size());
3709 std::vector<MPI_Request> first_requests(vertices_to_send.size());
3713 std::string>>>::iterator
3714 vert_to_send_it = vertices_to_send.begin(),
3715 vert_to_send_end = vertices_to_send.end();
3716 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3717 ++vert_to_send_it, ++i)
3719 int destination = vert_to_send_it->first;
3720 const unsigned int n_vertices = vert_to_send_it->second.size();
3721 const int buffer_size = 2 * n_vertices;
3722 vertices_send_buffers[i].resize(buffer_size);
3725 for (
unsigned int j = 0; j < n_vertices; ++j)
3727 vertices_send_buffers[i][2 * j] =
3728 std::get<0>(vert_to_send_it->second[j]);
3729 vertices_send_buffers[i][2 * j + 1] =
3730 local_to_global_vertex_index[std::get<1>(
3731 vert_to_send_it->second[j])];
3735 ierr = MPI_Isend(vertices_send_buffers[i].data(),
3741 &first_requests[i]);
3746 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
3747 vertices_to_recv.size());
3748 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
3749 vert_to_recv_it = vertices_to_recv.begin(),
3750 vert_to_recv_end = vertices_to_recv.end();
3751 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3752 ++vert_to_recv_it, ++i)
3754 int source = vert_to_recv_it->first;
3755 const unsigned int n_vertices = vert_to_recv_it->second.size();
3756 const int buffer_size = 2 * n_vertices;
3757 vertices_recv_buffers[i].resize(buffer_size);
3760 ierr = MPI_Recv(vertices_recv_buffers[i].data(),
3772 std::vector<std::vector<char>> cellids_send_buffers(
3773 vertices_to_send.size());
3774 std::vector<MPI_Request> second_requests(vertices_to_send.size());
3775 vert_to_send_it = vertices_to_send.begin();
3776 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3777 ++vert_to_send_it, ++i)
3779 int destination = vert_to_send_it->first;
3780 const unsigned int n_vertices = vert_to_send_it->second.size();
3781 const int buffer_size = max_cellid_size * n_vertices;
3782 cellids_send_buffers[i].resize(buffer_size);
3785 unsigned int pos = 0;
3786 for (
unsigned int j = 0; j < n_vertices; ++j)
3788 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
3789 for (
unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
3791 if (k < cell_id.size())
3792 cellids_send_buffers[i][pos] = cell_id[k];
3796 cellids_send_buffers[i][pos] =
'-';
3801 ierr = MPI_Isend(cellids_send_buffers[i].data(),
3807 &second_requests[i]);
3812 std::vector<std::vector<char>> cellids_recv_buffers(
3813 vertices_to_recv.size());
3814 vert_to_recv_it = vertices_to_recv.begin();
3815 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3816 ++vert_to_recv_it, ++i)
3818 int source = vert_to_recv_it->first;
3819 const unsigned int n_vertices = vert_to_recv_it->second.size();
3820 const int buffer_size = max_cellid_size * n_vertices;
3821 cellids_recv_buffers[i].resize(buffer_size);
3824 ierr = MPI_Recv(cellids_recv_buffers[i].data(),
3836 vert_to_recv_it = vertices_to_recv.begin();
3837 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3838 ++i, ++vert_to_recv_it)
3840 for (
unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
3842 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
3844 vertices_recv_buffers[i][2 * j + 1];
3845 const std::string cellid_recv(
3846 &cellids_recv_buffers[i][max_cellid_size * j],
3847 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
3849 typename std::set<active_cell_iterator>::iterator
3850 cell_set_it = missing_vert_cells.begin(),
3851 end_cell_set = missing_vert_cells.end();
3852 for (; (found ==
false) && (cell_set_it != end_cell_set);
3855 typename std::set<active_cell_iterator>::iterator
3857 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
3859 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
3860 for (; candidate_cell != end_cell; ++candidate_cell)
3862 std::string current_cellid =
3863 (*candidate_cell)->id().to_string();
3864 current_cellid.resize(max_cellid_size,
'-');
3865 if (current_cellid.compare(cellid_recv) == 0)
3867 local_to_global_vertex_index
3868 [(*candidate_cell)->vertex_index(local_pos_recv)] =
3880 return local_to_global_vertex_index;
3885 template <
int dim,
int spacedim>
3901 for (
const auto &cell :
triangulation.active_cell_iterators())
3903 const unsigned int index = cell->active_cell_index();
3904 cell_connectivity.
add(index, index);
3905 for (
auto f : cell->face_indices())
3906 if ((cell->at_boundary(f) ==
false) &&
3907 (cell->neighbor(f)->has_children() ==
false))
3909 const unsigned int other_index =
3910 cell->neighbor(f)->active_cell_index();
3911 cell_connectivity.
add(index, other_index);
3912 cell_connectivity.
add(other_index, index);
3919 template <
int dim,
int spacedim>
3925 std::vector<std::vector<unsigned int>> vertex_to_cell(
3927 for (
const auto &cell :
triangulation.active_cell_iterators())
3929 for (
const unsigned int v : cell->vertex_indices())
3930 vertex_to_cell[cell->vertex_index(v)].push_back(
3931 cell->active_cell_index());
3936 for (
const auto &cell :
triangulation.active_cell_iterators())
3938 for (
const unsigned int v : cell->vertex_indices())
3939 for (
unsigned int n = 0;
3940 n < vertex_to_cell[cell->vertex_index(v)].size();
3942 cell_connectivity.
add(cell->active_cell_index(),
3943 vertex_to_cell[cell->vertex_index(v)][n]);
3948 template <
int dim,
int spacedim>
3952 const unsigned int level,
3955 std::vector<std::vector<unsigned int>> vertex_to_cell(
3962 for (
const unsigned int v : cell->vertex_indices())
3963 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
3973 for (
const unsigned int v : cell->vertex_indices())
3974 for (
unsigned int n = 0;
3975 n < vertex_to_cell[cell->vertex_index(v)].size();
3977 cell_connectivity.
add(cell->index(),
3978 vertex_to_cell[cell->vertex_index(v)][n]);
3984 template <
int dim,
int spacedim>
3992 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3993 "are already partitioned implicitly and can not be "
3994 "partitioned again explicitly."));
3996 std::vector<unsigned int> cell_weights;
4006 for (
const auto &cell :
triangulation.active_cell_iterators())
4007 if (cell->is_locally_owned())
4008 cell_weights[cell->active_cell_index()] =
4018 if (
const auto shared_tria =
4022 shared_tria->get_communicator(),
4026 Assert(std::accumulate(cell_weights.begin(),
4028 std::uint64_t(0)) > 0,
4029 ExcMessage(
"The global sum of weights over all active cells "
4030 "is zero. Please verify how you generate weights."));
4042 template <
int dim,
int spacedim>
4045 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."));
4057 if (n_partitions == 1)
4059 for (
const auto &cell :
triangulation.active_cell_iterators())
4060 cell->set_subdomain_id(0);
4074 sp_cell_connectivity.
copy_from(cell_connectivity);
4077 sp_cell_connectivity,
4084 template <
int dim,
int spacedim>
4093 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4094 "are already partitioned implicitly and can not be "
4095 "partitioned again explicitly."));
4097 std::vector<unsigned int> cell_weights;
4107 for (
const auto &cell :
triangulation.active_cell_iterators() |
4109 cell_weights[cell->active_cell_index()] =
4119 if (
const auto shared_tria =
4123 shared_tria->get_communicator(),
4127 Assert(std::accumulate(cell_weights.begin(),
4129 std::uint64_t(0)) > 0,
4130 ExcMessage(
"The global sum of weights over all active cells "
4131 "is zero. Please verify how you generate weights."));
4137 cell_connection_graph,
4144 template <
int dim,
int spacedim>
4147 const std::vector<unsigned int> &cell_weights,
4154 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4155 "are already partitioned implicitly and can not be "
4156 "partitioned again explicitly."));
4159 ExcMessage(
"Connectivity graph has wrong size"));
4161 ExcMessage(
"Connectivity graph has wrong size"));
4167 if (n_partitions == 1)
4169 for (
const auto &cell :
triangulation.active_cell_iterators())
4170 cell->set_subdomain_id(0);
4178 std::vector<unsigned int> partition_indices(
triangulation.n_active_cells());
4186 for (
const auto &cell :
triangulation.active_cell_iterators())
4187 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
4199 unsigned int & current_proc_idx,
4200 unsigned int & current_cell_idx,
4202 const unsigned int n_partitions)
4204 if (cell->is_active())
4206 while (current_cell_idx >=
4208 (current_proc_idx + 1) / n_partitions))
4210 cell->set_subdomain_id(current_proc_idx);
4215 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4225 template <
int dim,
int spacedim>
4229 const bool group_siblings)
4233 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4234 "are already partitioned implicitly and can not be "
4235 "partitioned again explicitly."));
4243 if (n_partitions == 1)
4245 for (
const auto &cell :
triangulation.active_cell_iterators())
4246 cell->set_subdomain_id(0);
4252 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
4253 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
4259 coarse_cell_to_p4est_tree_permutation.resize(
triangulation.n_cells(0));
4261 coarse_cell_to_p4est_tree_permutation);
4263 p4est_tree_to_coarse_cell_permutation =
4266 unsigned int current_proc_idx = 0;
4267 unsigned int current_cell_idx = 0;
4272 for (
unsigned int idx = 0; idx <
triangulation.n_cells(0); ++idx)
4274 const unsigned int coarse_cell_idx =
4275 p4est_tree_to_coarse_cell_permutation[idx];
4298 for (; cell != endc; ++cell)
4300 if (cell->is_active())
4302 bool all_children_active =
true;
4303 std::map<unsigned int, unsigned int> map_cpu_n_cells;
4304 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4305 if (!cell->child(n)->is_active())
4307 all_children_active =
false;
4311 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
4313 if (!all_children_active)
4316 unsigned int new_owner = cell->child(0)->subdomain_id();
4317 for (std::map<unsigned int, unsigned int>::iterator it =
4318 map_cpu_n_cells.begin();
4319 it != map_cpu_n_cells.end();
4321 if (it->second > map_cpu_n_cells[new_owner])
4322 new_owner = it->first;
4324 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4325 cell->child(n)->set_subdomain_id(new_owner);
4331 template <
int dim,
int spacedim>
4336 for (
int lvl = n_levels - 1; lvl >= 0; --lvl)
4338 for (
const auto &cell :
triangulation.cell_iterators_on_level(lvl))
4340 if (cell->is_active())
4341 cell->set_level_subdomain_id(cell->subdomain_id());
4344 Assert(cell->child(0)->level_subdomain_id() !=
4347 cell->set_level_subdomain_id(
4348 cell->child(0)->level_subdomain_id());
4360 template <
int dim,
int spacedim>
4365 const std::vector<CellId> & cell_ids,
4366 std::vector<types::subdomain_id> &subdomain_ids)
4368 #ifndef DEAL_II_WITH_P4EST
4371 (void)subdomain_ids;
4375 "You are attempting to use a functionality that is only available "
4376 "if deal.II was configured to use p4est, but cmake did not find a "
4377 "valid p4est library."));
4382 for (
const auto &cell_id : cell_ids)
4385 typename ::internal::p4est::types<dim>::quadrant p4est_cell,
4388 ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
4389 for (
const auto &child_index : cell_id.get_child_indices())
4391 ::internal::p4est::init_quadrant_children<dim>(
4392 p4est_cell, p4est_children);
4394 p4est_children[
static_cast<unsigned int>(child_index)];
4400 const_cast<typename ::internal::p4est::types<dim>::forest
4402 cell_id.get_coarse_cell_id(),
4409 subdomain_ids.push_back(owner);
4416 template <
int spacedim>
4420 const std::vector<CellId> &,
4421 std::vector<types::subdomain_id> &)
4430 template <
int dim,
int spacedim>
4431 std::vector<types::subdomain_id>
4433 const std::vector<CellId> & cell_ids)
4435 std::vector<types::subdomain_id> subdomain_ids;
4436 subdomain_ids.reserve(cell_ids.size());
4445 *parallel_tria =
dynamic_cast<
4459 const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
4460 shared_tria->get_true_subdomain_ids_of_cells();
4462 for (
const auto &cell_id : cell_ids)
4464 const unsigned int active_cell_index =
4465 shared_tria->create_cell_iterator(cell_id)->active_cell_index();
4466 subdomain_ids.push_back(
4467 true_subdomain_ids_of_cells[active_cell_index]);
4474 for (
const auto &cell_id : cell_ids)
4476 subdomain_ids.push_back(
4477 triangulation.create_cell_iterator(cell_id)->subdomain_id());
4481 return subdomain_ids;
4486 template <
int dim,
int spacedim>
4489 std::vector<types::subdomain_id> & subdomain)
4494 for (
const auto &cell :
triangulation.active_cell_iterators())
4495 subdomain[cell->active_cell_index()] = cell->subdomain_id();
4500 template <
int dim,
int spacedim>
4506 unsigned int count = 0;
4507 for (
const auto &cell :
triangulation.active_cell_iterators())
4508 if (cell->subdomain_id() == subdomain)
4516 template <
int dim,
int spacedim>
4521 std::vector<bool> locally_owned_vertices =
4528 if (
const auto *tr =
dynamic_cast<
4531 for (
const auto &cell :
triangulation.active_cell_iterators())
4532 if (cell->is_artificial() ||
4533 (cell->is_ghost() &&
4534 (cell->subdomain_id() < tr->locally_owned_subdomain())))
4535 for (
const unsigned int v : cell->vertex_indices())
4536 locally_owned_vertices[cell->vertex_index(v)] =
false;
4538 return locally_owned_vertices;
4543 template <
int dim,
int spacedim>
4549 for (
const auto &cell :
triangulation.active_cell_iterators())
4550 if (!cell->is_artificial())
4551 min_diameter =
std::min(min_diameter, cell->diameter(mapping));
4553 double global_min_diameter = 0;
4555 #ifdef DEAL_II_WITH_MPI
4559 global_min_diameter =
4563 global_min_diameter = min_diameter;
4565 return global_min_diameter;
4570 template <
int dim,
int spacedim>
4575 double max_diameter = 0.;
4576 for (
const auto &cell :
triangulation.active_cell_iterators())
4577 if (!cell->is_artificial())
4578 max_diameter =
std::max(max_diameter, cell->diameter(mapping));
4580 double global_max_diameter = 0;
4582 #ifdef DEAL_II_WITH_MPI
4586 global_max_diameter =
4590 global_max_diameter = max_diameter;
4592 return global_max_diameter;
4599 namespace FixUpDistortedChildCells
4622 template <
typename Iterator,
int spacedim>
4627 const unsigned int structdim =
4628 Iterator::AccessorType::structure_dimension;
4629 Assert(spacedim == Iterator::AccessorType::dimension,
4635 Assert(object->refinement_case() ==
4643 Tensor<spacedim - structdim, spacedim>
4646 for (
const unsigned int i : object->vertex_indices())
4647 parent_vertices[i] =
object->vertex(i);
4650 parent_vertices, parent_alternating_forms);
4652 const Tensor<spacedim - structdim, spacedim>
4653 average_parent_alternating_form =
4654 std::accumulate(parent_alternating_forms,
4655 parent_alternating_forms +
4669 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4673 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4674 for (
const unsigned int i : object->child(c)->vertex_indices())
4675 child_vertices[c][i] =
object->child(c)->vertex(i);
4683 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4685 1] = object_mid_point;
4687 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4689 child_vertices[c], child_alternating_forms[c]);
4701 double objective = 0;
4702 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4703 for (
const unsigned int i : object->child(c)->vertex_indices())
4705 (child_alternating_forms[c][i] -
4706 average_parent_alternating_form / std::pow(2., 1. * structdim))
4718 template <
typename Iterator>
4721 const unsigned int f,
4722 std::integral_constant<int, 1>)
4724 return object->vertex(f);
4734 template <
typename Iterator>
4737 const unsigned int f,
4738 std::integral_constant<int, 2>)
4740 return object->line(f)->center();
4750 template <
typename Iterator>
4753 const unsigned int f,
4754 std::integral_constant<int, 3>)
4756 return object->face(f)->center();
4783 template <
typename Iterator>
4787 const unsigned int structdim =
4788 Iterator::AccessorType::structure_dimension;
4790 double diameter =
object->diameter();
4791 for (
const unsigned int f : object->face_indices())
4792 for (
unsigned int e = f + 1;
e <
object->n_faces(); ++
e)
4797 std::integral_constant<int, structdim>())
4799 object,
e, std::integral_constant<int, structdim>())));
4810 template <
typename Iterator>
4814 const unsigned int structdim =
4815 Iterator::AccessorType::structure_dimension;
4816 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
4823 Assert(object->refinement_case() ==
4833 unsigned int iteration = 0;
4838 double initial_delta = 0;
4845 const double step_length =
diameter / 4 / (iteration + 1);
4850 for (
unsigned int d = 0;
d < spacedim; ++
d)
4852 const double eps = step_length / 10;
4866 if (gradient.norm() == 0)
4875 std::min(2 * current_value / (gradient * gradient),
4876 step_length / gradient.norm()) *
4881 const double previous_value = current_value;
4885 initial_delta = (previous_value - current_value);
4888 if ((iteration >= 1) &&
4889 ((previous_value - current_value < 0) ||
4890 (
std::fabs(previous_value - current_value) <
4891 0.001 * initial_delta)))
4896 while (iteration < 20);
4912 double old_min_product, new_min_product;
4917 parent_vertices[i] = object->vertex(i);
4919 Tensor<spacedim - structdim, spacedim>
4922 parent_vertices, parent_alternating_forms);
4928 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4929 for (
const unsigned int i : object->child(c)->vertex_indices())
4930 child_vertices[c][i] =
object->child(c)->vertex(i);
4932 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4936 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4938 child_vertices[c], child_alternating_forms[c]);
4941 child_alternating_forms[0][0] * parent_alternating_forms[0];
4942 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4943 for (
const unsigned int i : object->child(c)->vertex_indices())
4944 for (
const unsigned int j : object->vertex_indices())
4945 old_min_product = std::min<double>(old_min_product,
4946 child_alternating_forms[c][i] *
4947 parent_alternating_forms[j]);
4955 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4957 1] = object_mid_point;
4959 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4961 child_vertices[c], child_alternating_forms[c]);
4964 child_alternating_forms[0][0] * parent_alternating_forms[0];
4965 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4966 for (
const unsigned int i : object->child(c)->vertex_indices())
4967 for (
const unsigned int j : object->vertex_indices())
4968 new_min_product = std::min<double>(new_min_product,
4969 child_alternating_forms[c][i] *
4970 parent_alternating_forms[j]);
4978 if (new_min_product >= old_min_product)
4979 object->child(0)->vertex(
4986 return (
std::max(new_min_product, old_min_product) > 0);
4992 template <
int dim,
int spacedim>
4995 const typename ::Triangulation<dim, spacedim>::cell_iterator
4997 std::integral_constant<int, dim>,
4998 std::integral_constant<int, spacedim>)
5008 for (
auto f : cell->face_indices())
5011 Assert(cell->face(f)->refinement_case() ==
5015 bool subface_is_more_refined =
false;
5016 for (
unsigned int g = 0;
5017 g < GeometryInfo<dim>::max_children_per_face;
5019 if (cell->face(f)->child(g)->has_children())
5021 subface_is_more_refined =
true;
5025 if (subface_is_more_refined ==
true)
5036 template <
int dim,
int spacedim>
5044 dim != 1 && spacedim != 1,
5045 "This function is only valid when dim != 1 or spacedim != 1.");
5049 for (
typename std::list<
5051 cell_ptr = distorted_cells.distorted_cells.
begin();
5052 cell_ptr != distorted_cells.distorted_cells.
end();
5058 Assert(!cell->is_active(),
5060 "This function is only valid for a list of cells that "
5061 "have children (i.e., no cell in the list may be active)."));
5065 std::integral_constant<int, dim>(),
5066 std::integral_constant<int, spacedim>());
5070 unfixable_subset.distorted_cells.push_back(cell);
5073 return unfixable_subset;
5078 template <
int dim,
int spacedim>
5081 const bool reset_boundary_ids)
5084 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
5085 auto m_it = dst_manifold_ids.begin();
5086 for (
const auto b : src_boundary_ids)
5091 const std::vector<types::boundary_id> reset_boundary_id =
5092 reset_boundary_ids ?
5093 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
5103 template <
int dim,
int spacedim>
5106 const std::vector<types::boundary_id> &src_boundary_ids,
5107 const std::vector<types::manifold_id> &dst_manifold_ids,
5109 const std::vector<types::boundary_id> &reset_boundary_ids_)
5112 const auto reset_boundary_ids =
5113 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
5122 for (
auto f : cell->face_indices())
5123 if (cell->face(f)->at_boundary())
5124 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
5126 const auto bid = cell->face(f)->line(
e)->boundary_id();
5127 const unsigned int ind = std::find(src_boundary_ids.begin(),
5128 src_boundary_ids.end(),
5130 src_boundary_ids.begin();
5131 if (ind < src_boundary_ids.size())
5132 cell->face(f)->line(
e)->set_manifold_id(
5133 dst_manifold_ids[ind]);
5138 for (
auto f : cell->face_indices())
5139 if (cell->face(f)->at_boundary())
5141 const auto bid = cell->face(f)->boundary_id();
5142 const unsigned int ind =
5143 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
5144 src_boundary_ids.begin();
5146 if (ind < src_boundary_ids.size())
5149 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
5151 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
5155 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
5157 const auto bid = cell->face(f)->line(
e)->boundary_id();
5158 const unsigned int ind = std::find(src_boundary_ids.begin(),
5159 src_boundary_ids.end(),
5161 src_boundary_ids.begin();
5162 if (ind < src_boundary_ids.size())
5163 cell->face(f)->line(
e)->set_boundary_id(
5164 reset_boundary_ids[ind]);
5170 template <
int dim,
int spacedim>
5173 const bool compute_face_ids)
5179 for (; cell != endc; ++cell)
5181 cell->set_manifold_id(cell->material_id());
5182 if (compute_face_ids ==
true)
5184 for (
auto f : cell->face_indices())
5186 if (cell->at_boundary(f) ==
false)
5187 cell->face(f)->set_manifold_id(
5189 cell->neighbor(f)->material_id()));
5191 cell->face(f)->set_manifold_id(cell->material_id());
5198 template <
int dim,
int spacedim>
5203 const std::set<types::manifold_id> &)> &disambiguation_function,
5204 bool overwrite_only_flat_manifold_ids)
5209 const unsigned int n_subobjects =
5213 std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
5214 std::vector<unsigned int> backup;
5218 unsigned next_index = 1;
5222 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
5224 if (cell->line(
l)->user_index() == 0)
5227 manifold_ids[next_index].insert(cell->manifold_id());
5228 cell->line(
l)->set_user_index(next_index++);
5231 manifold_ids[cell->line(
l)->user_index()].insert(
5232 cell->manifold_id());
5235 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
5237 if (cell->quad(
l)->user_index() == 0)
5240 manifold_ids[next_index].insert(cell->manifold_id());
5241 cell->quad(
l)->set_user_index(next_index++);
5244 manifold_ids[cell->quad(
l)->user_index()].insert(
5245 cell->manifold_id());
5251 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
5253 const auto id = cell->line(
l)->user_index();
5257 if (cell->line(
l)->manifold_id() ==
5259 overwrite_only_flat_manifold_ids ==
false)
5260 cell->line(
l)->set_manifold_id(
5261 disambiguation_function(manifold_ids[
id]));
5262 cell->line(
l)->set_user_index(0);
5266 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
5268 const auto id = cell->quad(
l)->user_index();
5272 if (cell->quad(
l)->manifold_id() ==
5274 overwrite_only_flat_manifold_ids ==
false)
5275 cell->quad(
l)->set_manifold_id(
5276 disambiguation_function(manifold_ids[
id]));
5277 cell->quad(
l)->set_user_index(0);
5286 template <
int dim,
int spacedim>
5287 std::pair<unsigned int, double>
5291 double max_ratio = 1;
5292 unsigned int index = 0;
5294 for (
unsigned int i = 0; i < dim; ++i)
5295 for (
unsigned int j = i + 1; j < dim; ++j)
5297 unsigned int ax = i % dim;
5298 unsigned int next_ax = j % dim;
5301 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
5303 if (ratio > max_ratio)
5308 else if (1.0 / ratio > max_ratio)
5310 max_ratio = 1.0 / ratio;
5314 return std::make_pair(index, max_ratio);
5318 template <
int dim,
int spacedim>
5321 const bool isotropic,
5322 const unsigned int max_iterations)
5324 unsigned int iter = 0;
5325 bool continue_refinement =
true;
5327 while (continue_refinement && (iter < max_iterations))
5331 continue_refinement =
false;
5334 for (
const unsigned int j : cell->face_indices())
5335 if (cell->at_boundary(j) ==
false &&
5336 cell->neighbor(j)->has_children())
5340 cell->set_refine_flag();
5341 continue_refinement =
true;
5344 continue_refinement |= cell->flag_for_face_refinement(j);
5351 template <
int dim,
int spacedim>
5354 const double max_ratio,
5355 const unsigned int max_iterations)
5357 unsigned int iter = 0;
5358 bool continue_refinement =
true;
5360 while (continue_refinement && (iter < max_iterations))
5363 continue_refinement =
false;
5366 std::pair<unsigned int, double> info =
5367 GridTools::get_longest_direction<dim, spacedim>(cell);
5368 if (info.second > max_ratio)
5370 cell->set_refine_flag(
5372 continue_refinement =
true;
5380 template <
int dim,
int spacedim>
5383 const double limit_angle_fraction)
5391 "have hanging nodes."));
5395 bool has_cells_with_more_than_dim_faces_on_boundary =
true;
5396 bool has_cells_with_dim_faces_on_boundary =
false;
5398 unsigned int refinement_cycles = 0;
5400 while (has_cells_with_more_than_dim_faces_on_boundary)
5402 has_cells_with_more_than_dim_faces_on_boundary =
false;
5406 unsigned int boundary_face_counter = 0;
5407 for (
auto f : cell->face_indices())
5408 if (cell->face(f)->at_boundary())
5409 boundary_face_counter++;
5410 if (boundary_face_counter > dim)
5412 has_cells_with_more_than_dim_faces_on_boundary =
true;
5415 else if (boundary_face_counter == dim)
5416 has_cells_with_dim_faces_on_boundary =
true;
5418 if (has_cells_with_more_than_dim_faces_on_boundary)
5421 refinement_cycles++;
5425 if (has_cells_with_dim_faces_on_boundary)
5428 refinement_cycles++;
5432 while (refinement_cycles > 0)
5435 cell->set_coarsen_flag();
5437 refinement_cycles--;
5447 std::vector<CellData<dim>> cells_to_add;
5451 const unsigned int v0 = 0,
v1 = 1, v2 = (dim > 1 ? 2 : 0),
5452 v3 = (dim > 1 ? 3 : 0);
5456 double angle_fraction = 0;
5462 p0[spacedim > 1 ? 1 : 0] = 1;
5466 if (cell->face(
v0)->at_boundary() && cell->face(v3)->at_boundary())
5468 p0 = cell->vertex(
v0) - cell->vertex(v2);
5469 p1 = cell->vertex(v3) - cell->vertex(v2);
5470 vertex_at_corner = v2;
5472 else if (cell->face(v3)->at_boundary() &&
5473 cell->face(
v1)->at_boundary())
5475 p0 = cell->vertex(v2) - cell->vertex(v3);
5476 p1 = cell->vertex(
v1) - cell->vertex(v3);
5477 vertex_at_corner = v3;
5479 else if (cell->face(1)->at_boundary() &&
5480 cell->face(2)->at_boundary())
5482 p0 = cell->vertex(
v0) - cell->vertex(
v1);
5483 p1 = cell->vertex(v3) - cell->vertex(
v1);
5484 vertex_at_corner =
v1;
5486 else if (cell->face(2)->at_boundary() &&
5487 cell->face(0)->at_boundary())
5489 p0 = cell->vertex(v2) - cell->vertex(
v0);
5490 p1 = cell->vertex(
v1) - cell->vertex(
v0);
5491 vertex_at_corner =
v0;
5502 if (angle_fraction > limit_angle_fraction)
5504 auto flags_removal = [&](
unsigned int f1,
5507 unsigned int n2) ->
void {
5508 cells_to_remove[cell->active_cell_index()] =
true;
5509 cells_to_remove[cell->neighbor(n1)->active_cell_index()] =
true;
5510 cells_to_remove[cell->neighbor(n2)->active_cell_index()] =
true;
5512 faces_to_remove[cell->face(f1)->index()] =
true;
5513 faces_to_remove[cell->face(f2)->index()] =
true;
5515 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] =
true;
5516 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] =
true;
5519 auto cell_creation = [&](
const unsigned int vv0,
5520 const unsigned int vv1,
5521 const unsigned int f0,
5522 const unsigned int f1,
5524 const unsigned int n0,
5525 const unsigned int v0n0,
5526 const unsigned int v1n0,
5528 const unsigned int n1,
5529 const unsigned int v0n1,
5530 const unsigned int v1n1) {
5536 c1.
vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
5537 c1.
vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
5543 c2.
vertices[
v1] = cell->neighbor(n1)->vertex_index(v0n1);
5544 c2.
vertices[v2] = cell->vertex_index(vv1);
5545 c2.
vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
5550 l1.
vertices[0] = cell->vertex_index(vv0);
5551 l1.
vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
5557 l2.
vertices[0] = cell->vertex_index(vv0);
5558 l2.
vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
5564 cells_to_add.push_back(c1);
5565 cells_to_add.push_back(c2);
5570 switch (vertex_at_corner)
5573 flags_removal(0, 2, 3, 1);
5574 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
5577 flags_removal(1, 2, 3, 0);
5578 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
5581 flags_removal(3, 0, 1, 2);
5582 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
5585 flags_removal(3, 1, 0, 2);
5586 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
5599 if (cells_to_add.size() == 0)
5601 while (refinement_cycles > 0)
5604 cell->set_coarsen_flag();
5606 refinement_cycles--;
5614 if (cells_to_remove[cell->active_cell_index()] ==
false)
5617 for (
const unsigned int v : cell->vertex_indices())
5618 c.
vertices[v] = cell->vertex_index(v);
5621 cells_to_add.push_back(c);
5629 for (; face != endf; ++face)
5630 if ((face->at_boundary() ||
5632 faces_to_remove[face->index()] ==
false)
5634 for (
unsigned int l = 0;
l < face->
n_lines(); ++
l)
5639 for (
const unsigned int v : face->vertex_indices())
5640 line.
vertices[v] = face->vertex_index(v);
5646 for (
const unsigned int v : face->line(
l)->vertex_indices())
5647 line.
vertices[v] = face->line(
l)->vertex_index(v);
5656 for (
const unsigned int v : face->vertex_indices())
5657 quad.
vertices[v] = face->vertex_index(v);
5665 subcelldata_to_add);
5670 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
5689 template <
int dim,
int spacedim>
5692 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5693 std::vector<std::vector<Point<dim>>>,
5694 std::vector<std::vector<unsigned int>>>
5706 auto &cells = std::get<0>(cqmp);
5707 auto &qpoints = std::get<1>(cqmp);
5708 auto &maps = std::get<2>(cqmp);
5710 return std::make_tuple(std::move(cells),
5717 template <
int dim,
int spacedim>
5720 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5721 std::vector<std::vector<Point<dim>>>,
5722 std::vector<std::vector<unsigned int>>,
5723 std::vector<unsigned int>>
5733 Assert((dim == spacedim),
5734 ExcMessage(
"Only implemented for dim==spacedim."));
5743 const unsigned int np = points.size();
5745 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5747 std::vector<std::vector<Point<dim>>> qpoints_out;
5748 std::vector<std::vector<unsigned int>> maps_out;
5749 std::vector<unsigned int> missing_points_out;
5753 return std::make_tuple(std::move(cells_out),
5754 std::move(qpoints_out),
5755 std::move(maps_out),
5756 std::move(missing_points_out));
5764 std::vector<std::pair<Point<spacedim>,
unsigned int>> points_and_ids(np);
5765 for (
unsigned int i = 0; i < np; ++i)
5766 points_and_ids[i] = std::make_pair(points[i], i);
5767 const auto p_tree =
pack_rtree(points_and_ids);
5770 std::vector<bool> found_points(points.size(),
false);
5773 const auto already_found = [&found_points](
const auto &id) {
5775 return found_points[
id.second];
5781 const auto store_cell_point_and_id =
5785 const unsigned int &id) {
5786 const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
5787 if (it != cells_out.rend())
5789 const auto cell_id =
5790 (cells_out.size() - 1 - (it - cells_out.rbegin()));
5791 qpoints_out[cell_id].emplace_back(ref_point);
5792 maps_out[cell_id].emplace_back(
id);
5796 cells_out.emplace_back(cell);
5797 qpoints_out.emplace_back(std::vector<
Point<dim>>({ref_point}));
5798 maps_out.emplace_back(std::vector<unsigned int>({
id}));
5803 const auto check_all_points_within_box = [&](
const auto &leaf) {
5804 const auto &box = leaf.first;
5805 const auto &cell_hint = leaf.second;
5807 for (
const auto &point_and_id :
5808 p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
5809 bgi::intersects(box)))
5811 const auto id = point_and_id.second;
5812 const auto cell_and_ref =
5816 const auto &cell = cell_and_ref.first;
5817 const auto &ref_point = cell_and_ref.second;
5820 store_cell_point_and_id(cell, ref_point,
id);
5822 missing_points_out.emplace_back(
id);
5825 found_points[id] =
true;
5831 check_all_points_within_box(
5832 std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
5835 for (
unsigned int i = 0; i < np; ++i)
5836 if (found_points[i] ==
false)
5839 const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
5841 if (leaf != b_tree.qend())
5842 check_all_points_within_box(*leaf);
5850 for (
unsigned int i = 0; i < np; ++i)
5851 if (found_points[i] ==
false)
5852 missing_points_out.emplace_back(i);
5859 unsigned int c = cells_out.size();
5860 unsigned int qps = 0;
5866 for (
unsigned int n = 0; n < c; ++n)
5869 qps += qpoints_out[n].size();
5872 Assert(qps + missing_points_out.size() == np,
5876 return std::make_tuple(std::move(cells_out),
5877 std::move(qpoints_out),
5878 std::move(maps_out),
5879 std::move(missing_points_out));
5884 template <
int dim,
int spacedim>
5887 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5888 std::vector<std::vector<Point<dim>>>,
5889 std::vector<std::vector<unsigned int>>,
5890 std::vector<std::vector<Point<spacedim>>>,
5891 std::vector<std::vector<unsigned int>>>
5899 const double tolerance)
5903 cache, points, global_bboxes, {}, tolerance,
false,
true)
5908 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5909 std::vector<std::vector<Point<dim>>>,
5910 std::vector<std::vector<unsigned int>>,
5911 std::vector<std::vector<Point<spacedim>>>,
5912 std::vector<std::vector<unsigned int>>>
5915 std::pair<int, int> dummy{-1, -1};
5917 for (
unsigned int i = 0; i < all.size(); ++i)
5919 if (dummy != std::get<0>(all[i]))
5921 std::get<0>(result).push_back(
5924 std::get<0>(all[i]).
first,
5925 std::get<0>(all[i]).
second});
5927 const unsigned int new_size = std::get<0>(result).size();
5929 std::get<1>(result).resize(new_size);
5930 std::get<2>(result).resize(new_size);
5931 std::get<3>(result).resize(new_size);
5932 std::get<4>(result).resize(new_size);
5934 dummy = std::get<0>(all[i]);
5937 std::get<1>(result).back().push_back(
5938 std::get<3>(all[i]));
5939 std::get<2>(result).back().push_back(std::get<2>(all[i]));
5940 std::get<3>(result).back().push_back(std::get<4>(all[i]));
5941 std::get<4>(result).back().push_back(std::get<1>(all[i]));
5951 template <
int spacedim>
5952 std::tuple<std::vector<unsigned int>,
5953 std::vector<unsigned int>,
5954 std::vector<unsigned int>>
5958 const double tolerance)
5960 std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
5961 ranks_and_indices.reserve(points.size());
5963 for (
unsigned int i = 0; i < points.size(); ++i)
5965 const auto &
point = points[i];
5966 for (
unsigned rank = 0; rank < global_bboxes.size(); ++rank)
5967 for (
const auto &box : global_bboxes[rank])
5968 if (box.point_inside(
point, tolerance))
5970 ranks_and_indices.emplace_back(rank, i);
5976 std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
5978 std::vector<unsigned int> ranks;
5979 std::vector<unsigned int> ptr;
5980 std::vector<unsigned int> indices;
5984 for (
const auto &i : ranks_and_indices)
5986 if (dummy_rank != i.first)
5988 dummy_rank = i.first;
5989 ranks.push_back(dummy_rank);
5990 ptr.push_back(indices.size());
5993 indices.push_back(i.second);
5995 ptr.push_back(indices.size());
5997 return std::make_tuple(std::move(ranks),
5999 std::move(indices));
6004 template <
int dim,
int spacedim>
6006 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
6012 const std::vector<bool> &marked_vertices,
6013 const double tolerance,
6014 const bool enforce_unique_mapping)
6017 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
6019 locally_owned_active_cells_around_point;
6036 cell_hint = first_cell.first;
6039 const auto active_cells_around_point =
6047 if (enforce_unique_mapping)
6054 for (
const auto &cell : active_cells_around_point)
6055 lowes_rank =
std::min(lowes_rank, cell.first->subdomain_id());
6057 if (lowes_rank != my_rank)
6061 locally_owned_active_cells_around_point.reserve(
6062 active_cells_around_point.size());
6064 for (
const auto &cell : active_cells_around_point)
6065 if (cell.first->is_locally_owned())
6066 locally_owned_active_cells_around_point.push_back(cell);
6069 std::sort(locally_owned_active_cells_around_point.begin(),
6070 locally_owned_active_cells_around_point.end(),
6071 [](
const auto &a,
const auto &
b) { return a.first < b.first; });
6073 if (enforce_unique_mapping &&
6074 locally_owned_active_cells_around_point.size() > 1)
6076 return {locally_owned_active_cells_around_point.front()};
6078 return locally_owned_active_cells_around_point;
6083 template <
int dim,
int spacedim>
6084 DistributedComputePointLocationsInternal<dim, spacedim>
6089 const std::vector<bool> & marked_vertices,
6090 const double tolerance,
6091 const bool perform_handshake,
6092 const bool enforce_unique_mapping)
6103 const auto potential_owners =
6106 const auto &potential_owners_ranks = std::get<0>(potential_owners);
6107 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
6108 const auto &potential_owners_indices = std::get<2>(potential_owners);
6112 const auto translate = [&](
const unsigned int other_rank) {
6113 const auto ptr = std::find(potential_owners_ranks.begin(),
6114 potential_owners_ranks.end(),
6119 const auto other_rank_index =
6120 std::distance(potential_owners_ranks.begin(), ptr);
6122 return other_rank_index;
6126 (marked_vertices.size() == 0) ||
6129 "The marked_vertices vector has to be either empty or its size has "
6130 "to equal the number of vertices of the triangulation."));
6132 using RequestType = std::vector<std::pair<unsigned int, Point<spacedim>>>;
6133 using AnswerType = std::vector<unsigned int>;
6139 const bool has_relevant_vertices =
6140 (marked_vertices.size() == 0) ||
6141 (std::find(marked_vertices.begin(), marked_vertices.end(),
true) !=
6142 marked_vertices.end());
6144 const auto create_request = [&](
const unsigned int other_rank) {
6145 const auto other_rank_index = translate(other_rank);
6147 RequestType request;
6148 request.reserve(potential_owners_ptrs[other_rank_index + 1] -
6149 potential_owners_ptrs[other_rank_index]);
6151 for (
unsigned int i = potential_owners_ptrs[other_rank_index];
6152 i < potential_owners_ptrs[other_rank_index + 1];
6154 request.emplace_back(potential_owners_indices[i],
6155 points[potential_owners_indices[i]]);
6160 const auto answer_request =
6161 [&](
const unsigned int &other_rank,
6162 const RequestType & request) -> AnswerType {
6163 AnswerType answer(request.size(), 0);
6165 if (has_relevant_vertices)
6169 for (
unsigned int i = 0; i < request.size(); ++i)
6171 const auto &index_and_point = request[i];
6173 const auto cells_and_reference_positions =
6176 index_and_point.second,
6180 enforce_unique_mapping);
6182 for (
const auto &cell_and_reference_position :
6183 cells_and_reference_positions)
6185 const auto cell = cell_and_reference_position.first;
6186 auto reference_position =
6187 cell_and_reference_position.second;
6191 if (cell->reference_cell().is_hyper_cube())
6192 reference_position =
6194 reference_position);
6196 send_components.emplace_back(
6197 std::pair<int, int>(cell->level(), cell->index()),
6199 index_and_point.first,
6201 index_and_point.second,
6205 answer[i] = cells_and_reference_positions.size();
6209 if (perform_handshake)
6215 const auto process_answer = [&](
const unsigned int other_rank,
6216 const AnswerType & answer) {
6217 if (perform_handshake)
6219 const auto other_rank_index = translate(other_rank);
6221 for (
unsigned int i = 0; i < answer.size(); ++i)
6222 for (
unsigned int j = 0; j < answer[i]; ++j)
6223 recv_components.emplace_back(
6225 potential_owners_indices
6226 [i + potential_owners_ptrs[other_rank_index]],
6231 Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
6232 potential_owners_ranks,
6242 std::sort(send_components.begin(),
6243 send_components.end(),
6244 [&](
const auto &a,
const auto &
b) {
6245 if (std::get<1>(a) != std::get<1>(b))
6246 return std::get<1>(a) < std::get<1>(b);
6248 if (std::get<2>(a) != std::get<2>(b))
6249 return std::get<2>(a) < std::get<2>(b);
6251 return std::get<0>(a) < std::get<0>(b);
6256 i < send_components.size();
6259 std::get<5>(send_components[i]) = i;
6261 if (dummy != std::get<1>(send_components[i]))
6263 dummy = std::get<1>(send_components[i]);
6264 send_ranks.push_back(dummy);
6265 send_ptrs.push_back(i);
6268 send_ptrs.push_back(send_components.size());
6272 std::sort(send_components.begin(),
6273 send_components.end(),
6274 [&](
const auto &a,
const auto &
b) {
6275 if (std::get<0>(a) != std::get<0>(b))
6276 return std::get<0>(a) < std::get<0>(b);
6278 if (std::get<1>(a) != std::get<1>(b))
6279 return std::get<1>(a) < std::get<1>(b);
6281 if (std::get<2>(a) != std::get<2>(b))
6282 return std::get<2>(a) < std::get<2>(b);
6284 return std::get<5>(a) < std::get<5>(b);
6288 if (perform_handshake)
6291 std::sort(recv_components.begin(),
6292 recv_components.end(),
6293 [&](
const auto &a,
const auto &
b) {
6294 if (std::get<0>(a) != std::get<0>(b))
6295 return std::get<0>(a) < std::get<0>(b);
6297 return std::get<1>(a) < std::get<1>(b);
6302 i < recv_components.size();
6305 std::get<2>(recv_components[i]) = i;
6307 if (dummy != std::get<0>(recv_components[i]))
6309 dummy = std::get<0>(recv_components[i]);
6310 recv_ranks.push_back(dummy);
6311 recv_ptrs.push_back(i);
6314 recv_ptrs.push_back(recv_components.size());
6318 std::sort(recv_components.begin(),
6319 recv_components.end(),
6320 [&](
const auto &a,
const auto &
b) {
6321 if (std::get<1>(a) != std::get<1>(b))
6322 return std::get<1>(a) < std::get<1>(b);
6324 if (std::get<0>(a) != std::get<0>(b))
6325 return std::get<0>(a) < std::get<0>(b);
6327 return std::get<2>(a) < std::get<2>(b);
6337 template <
int dim,
int spacedim>
6338 std::map<unsigned int, Point<spacedim>>
6342 std::map<unsigned int, Point<spacedim>> result;