17 #include <deal.II/base/mpi.templates.h>
19 #include <deal.II/base/mpi_consensus_algorithms.templates.h>
64 #include <boost/random/mersenne_twister.hpp>
65 #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,
749 std::vector<unsigned int> new_vertex_numbers(
vertices.size());
750 std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
753 if (considered_vertices.size() == 0)
754 considered_vertices = new_vertex_numbers;
769 unsigned int longest_coordinate_direction = 0;
770 double longest_coordinate_length = bbox.
side_length(0);
771 for (
unsigned int d = 1;
d < spacedim; ++
d)
774 if (longest_coordinate_length < coordinate_length)
776 longest_coordinate_length = coordinate_length;
777 longest_coordinate_direction =
d;
783 std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
784 sorted_vertices.reserve(
vertices.size());
785 for (
const unsigned int vertex_n : considered_vertices)
788 sorted_vertices.emplace_back(vertex_n,
vertices[vertex_n]);
790 std::sort(sorted_vertices.begin(),
791 sorted_vertices.end(),
794 return a.second[longest_coordinate_direction] <
795 b.second[longest_coordinate_direction];
800 for (
unsigned int d = 0;
d < spacedim; ++
d)
801 if (std::abs(a[
d] -
b[
d]) > tol)
808 auto range_start = sorted_vertices.begin();
809 while (range_start != sorted_vertices.end())
811 auto range_end = range_start + 1;
812 while (range_end != sorted_vertices.end() &&
813 std::abs(range_end->second[longest_coordinate_direction] -
814 range_start->second[longest_coordinate_direction]) <
820 std::sort(range_start,
824 return a.first <
b.first;
833 for (
auto reference = range_start; reference != range_end; ++reference)
836 for (
auto it = reference + 1; it != range_end; ++it)
838 if (within_tolerance(reference->second, it->second))
840 new_vertex_numbers[it->first] = reference->first;
846 range_start = range_end;
853 for (
auto &cell : cells)
854 for (
auto &vertex_index : cell.vertices)
855 vertex_index = new_vertex_numbers[vertex_index];
857 for (
auto &vertex_index : quad.vertices)
858 vertex_index = new_vertex_numbers[vertex_index];
860 for (
auto &vertex_index : line.vertices)
861 vertex_index = new_vertex_numbers[vertex_index];
868 template <
int dim,
int spacedim>
879 if (dim == 2 && spacedim == 3)
882 std::size_t n_negative_cells = 0;
883 for (
auto &cell : cells)
888 const unsigned int n_vertices =
vertices.size();
898 std::swap(cell.vertices[1], cell.vertices[2]);
903 std::swap(cell.vertices[0], cell.vertices[2]);
904 std::swap(cell.vertices[1], cell.vertices[3]);
905 std::swap(cell.vertices[4], cell.vertices[6]);
906 std::swap(cell.vertices[5], cell.vertices[7]);
918 cell.vertices[n_vertices - 1]);
931 return n_negative_cells;
935 template <
int dim,
int spacedim>
941 const std::size_t n_negative_cells =
950 AssertThrow(n_negative_cells == 0 || n_negative_cells == cells.size(),
953 "This function assumes that either all cells have positive "
954 "volume, or that all cells have been specified in an "
955 "inverted vertex order so that their volume is negative. "
956 "(In the latter case, this class automatically inverts "
957 "every cell.) However, the mesh you have specified "
958 "appears to have both cells with positive and cells with "
959 "negative volume. You need to check your mesh which "
960 "cells these are and how they got there.\n"
961 "As a hint, of the total ") +
964 " appear to have a negative volume."));
982 CheapEdge(
const unsigned int v0,
const unsigned int v1)
994 return ((
v0 <
e.v0) || ((
v0 ==
e.v0) && (
v1 <
e.v1)));
1017 std::set<CheapEdge> edges;
1019 for (
typename std::vector<
CellData<dim>>::const_iterator c =
1029 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1031 const CheapEdge reverse_edge(
1034 if (edges.find(reverse_edge) != edges.end())
1039 const CheapEdge correct_edge(
1042 edges.insert(correct_edge);
1058 struct ParallelEdges
1072 static const unsigned int
1137 class AdjacentCells;
1145 class AdjacentCells<2>
1152 using const_iterator =
const AdjacentCell *;
1163 push_back(
const AdjacentCell &adjacent_cell)
1227 class AdjacentCells<3> :
public std::vector<AdjacentCell>
1249 Edge(
const CellData<dim> &cell,
const unsigned int edge_number)
1300 enum OrientationStatus
1330 Cell(
const CellData<dim> &c,
const std::vector<Edge<dim>> &edge_list)
1337 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1339 const Edge<dim>
e(c,
l);
1375 class EdgeDeltaSet<2>
1381 using const_iterator =
const unsigned int *;
1407 insert(
const unsigned int edge_index)
1468 class EdgeDeltaSet<3> :
public std::set<unsigned int>
1478 std::vector<Edge<dim>>
1485 std::vector<Edge<dim>> edge_list;
1487 for (
unsigned int i = 0; i < cells.size(); ++i)
1488 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1489 edge_list.emplace_back(cells[i],
l);
1492 std::sort(edge_list.begin(), edge_list.end());
1493 edge_list.erase(std::unique(edge_list.begin(), edge_list.end()),
1506 std::vector<Cell<dim>>
1507 build_cells_and_connect_edges(
const std::vector<
CellData<dim>> &cells,
1508 std::vector<Edge<dim>> & edges)
1510 std::vector<Cell<dim>> cell_list;
1511 cell_list.reserve(cells.size());
1512 for (
unsigned int i = 0; i < cells.size(); ++i)
1516 cell_list.emplace_back(cells[i], edges);
1520 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1521 edges[cell_list.back().edge_indices[
l]].adjacent_cells.push_back(
1522 AdjacentCell(i,
l));
1537 get_next_unoriented_cell(
const std::vector<Cell<dim>> &cells,
1538 const std::vector<Edge<dim>> &edges,
1539 const unsigned int current_cell)
1541 for (
unsigned int c = current_cell; c < cells.size(); ++c)
1542 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1544 Edge<dim>::not_oriented)
1559 orient_one_set_of_parallel_edges(
const std::vector<Cell<dim>> &cells,
1560 std::vector<Edge<dim>> & edges,
1561 const unsigned int cell,
1562 const unsigned int local_edge)
1589 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1590 (dim == 2 ? Edge<dim>::backward : Edge<dim>::forward);
1606 edges[cells[cell].edge_indices[local_edge]].orientation_status =
1607 (dim == 2 ? Edge<dim>::forward : Edge<dim>::backward);
1616 EdgeDeltaSet<dim> Delta_k;
1617 EdgeDeltaSet<dim> Delta_k_minus_1;
1618 Delta_k_minus_1.insert(cells[cell].
edge_indices[local_edge]);
1620 while (Delta_k_minus_1.begin() !=
1621 Delta_k_minus_1.end())
1625 for (
typename EdgeDeltaSet<dim>::const_iterator delta =
1626 Delta_k_minus_1.begin();
1627 delta != Delta_k_minus_1.end();
1631 Edge<dim>::not_oriented,
1635 for (
typename AdjacentCells<dim>::const_iterator adjacent_cell =
1637 adjacent_cell != edges[*delta].adjacent_cells.end();
1640 const unsigned int K = adjacent_cell->cell_index;
1641 const unsigned int delta_is_edge_in_K =
1642 adjacent_cell->edge_within_cell;
1646 const unsigned int first_edge_vertex =
1647 (edges[*delta].orientation_status == Edge<dim>::forward ?
1648 edges[*delta].vertex_indices[0] :
1649 edges[*delta].vertex_indices[1]);
1650 const unsigned int first_edge_vertex_in_K =
1653 delta_is_edge_in_K, 0)];
1655 first_edge_vertex == first_edge_vertex_in_K ||
1656 first_edge_vertex ==
1658 dim>::line_to_cell_vertices(delta_is_edge_in_K, 1)],
1663 for (
unsigned int o_e = 0;
1670 const unsigned int opposite_edge =
1671 cells[K].edge_indices[ParallelEdges<
1673 const unsigned int first_opposite_edge_vertex =
1674 cells[K].vertex_indices
1678 (first_edge_vertex == first_edge_vertex_in_K ? 0 :
1685 const typename Edge<dim>::OrientationStatus
1686 opposite_edge_orientation =
1687 (edges[opposite_edge].vertex_indices[0] ==
1688 first_opposite_edge_vertex ?
1689 Edge<dim>::forward :
1690 Edge<dim>::backward);
1695 Edge<dim>::not_oriented)
1699 edges[opposite_edge].orientation_status =
1700 opposite_edge_orientation;
1701 Delta_k.insert(opposite_edge);
1717 opposite_edge_orientation,
1723 opposite_edge_orientation)
1736 Delta_k_minus_1 = Delta_k;
1750 rotate_cell(
const std::vector<Cell<dim>> &cell_list,
1751 const std::vector<Edge<dim>> &edge_list,
1758 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++
e)
1765 starting_vertex_of_edge[
e] =
1769 starting_vertex_of_edge[
e] =
1785 if ((starting_vertex_of_edge[0] == starting_vertex_of_edge[2]) ||
1786 (starting_vertex_of_edge[0] == starting_vertex_of_edge[3]))
1787 origin_vertex_of_cell = starting_vertex_of_edge[0];
1788 else if ((starting_vertex_of_edge[1] ==
1789 starting_vertex_of_edge[2]) ||
1790 (starting_vertex_of_edge[1] ==
1791 starting_vertex_of_edge[3]))
1792 origin_vertex_of_cell = starting_vertex_of_edge[1];
1804 for (origin_vertex_of_cell = 0;
1805 origin_vertex_of_cell < GeometryInfo<dim>::vertices_per_cell;
1806 ++origin_vertex_of_cell)
1807 if (std::count(starting_vertex_of_edge,
1808 starting_vertex_of_edge +
1813 Assert(origin_vertex_of_cell <
1838 const unsigned int tmp = raw_cells[
cell_index].vertices[0];
1861 static const unsigned int cube_permutations[8][8] = {
1862 {0, 1, 2, 3, 4, 5, 6, 7},
1863 {1, 5, 3, 7, 0, 4, 2, 6},
1864 {2, 6, 0, 4, 3, 7, 1, 5},
1865 {3, 2, 1, 0, 7, 6, 5, 4},
1866 {4, 0, 6, 2, 5, 1, 7, 3},
1867 {5, 4, 7, 6, 1, 0, 3, 2},
1868 {6, 7, 4, 5, 2, 3, 0, 1},
1869 {7, 3, 5, 1, 6, 2, 4, 0}};
1874 temp_vertex_indices[v] =
1876 .vertices[cube_permutations[origin_vertex_of_cell][v]];
1878 raw_cells[
cell_index].vertices[v] = temp_vertex_indices[v];
1902 std::vector<Edge<dim>> edge_list = build_edges(cells);
1903 std::vector<Cell<dim>> cell_list =
1904 build_cells_and_connect_edges(cells, edge_list);
1908 unsigned int next_cell_with_unoriented_edge = 0;
1909 while ((next_cell_with_unoriented_edge = get_next_unoriented_cell(
1910 cell_list, edge_list, next_cell_with_unoriented_edge)) !=
1924 for (
unsigned int l = 0;
l < dim; ++
l)
1925 if (edge_list[cell_list[next_cell_with_unoriented_edge]
1928 orient_one_set_of_parallel_edges(
1931 next_cell_with_unoriented_edge,
1936 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1937 Assert(edge_list[cell_list[next_cell_with_unoriented_edge]
1946 for (
unsigned int c = 0; c < cells.size(); ++c)
1947 rotate_cell(cell_list, edge_list, c, cells);
1962 Assert(cells.size() != 0,
1964 "List of elements to orient must have at least one cell"));
1972 if (!is_consistent(cells))
1991 template <
int spacedim>
2030 template <
int spacedim>
2049 template <
int dim,
int spacedim>
2071 const unsigned int axis,
2083 template <
int dim,
int spacedim>
2105 const unsigned int n_dofs = S.
n();
2117 const auto constrained_rhs =
2119 solver.
solve(SF, u, constrained_rhs, prec);
2132 const bool solve_for_absolute_positions)
2160 std::array<AffineConstraints<double>, dim> constraints;
2161 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
2169 for (
const unsigned int vertex_no : cell->vertex_indices())
2171 const unsigned int vertex_index = cell->vertex_index(vertex_no);
2172 const Point<dim> & vertex_point = cell->vertex(vertex_no);
2174 const typename std::map<unsigned int, Point<dim>>::const_iterator
2175 map_iter = new_points.find(vertex_index);
2177 if (map_iter != map_end)
2178 for (
unsigned int i = 0; i < dim; ++i)
2180 constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
2181 constraints[i].set_inhomogeneity(
2182 cell->vertex_dof_index(vertex_no, 0),
2183 (solve_for_absolute_positions ?
2184 map_iter->second(i) :
2185 map_iter->second(i) - vertex_point[i]));
2190 for (
unsigned int i = 0; i < dim; ++i)
2191 constraints[i].close();
2195 for (
unsigned int i = 0; i < dim; ++i)
2200 for (
unsigned int i = 0; i < dim; ++i)
2207 std::vector<bool> vertex_touched(
triangulation.n_vertices(),
false);
2209 for (
const unsigned int vertex_no : cell->vertex_indices())
2210 if (vertex_touched[cell->vertex_index(vertex_no)] ==
false)
2215 cell->vertex_dof_index(vertex_no, 0);
2216 for (
unsigned int i = 0; i < dim; ++i)
2217 if (solve_for_absolute_positions)
2218 v(i) = us[i](dof_index);
2220 v(i) += us[i](dof_index);
2222 vertex_touched[cell->vertex_index(vertex_no)] =
true;
2226 template <
int dim,
int spacedim>
2227 std::map<unsigned int, Point<spacedim>>
2230 std::map<unsigned int, Point<spacedim>> vertex_map;
2234 for (; cell != endc; ++cell)
2236 for (
unsigned int i : cell->face_indices())
2240 if (face->at_boundary())
2242 for (
unsigned j = 0; j < face->
n_vertices(); ++j)
2245 const unsigned int vertex_index = face->vertex_index(j);
2246 vertex_map[vertex_index] = vertex;
2258 template <
int dim,
int spacedim>
2262 const bool keep_boundary,
2263 const unsigned int seed)
2280 double almost_infinite_length = 0;
2285 almost_infinite_length += cell->diameter();
2287 std::vector<double> minimal_length(
triangulation.n_vertices(),
2288 almost_infinite_length);
2291 std::vector<bool> at_boundary(keep_boundary ?
triangulation.n_vertices() :
2296 const bool is_parallel_shared =
2299 for (
const auto &cell :
triangulation.active_cell_iterators())
2300 if (is_parallel_shared || cell->is_locally_owned())
2304 for (
unsigned int i = 0; i < cell->n_lines(); ++i)
2307 line = cell->line(i);
2309 if (keep_boundary && line->at_boundary())
2311 at_boundary[line->vertex_index(0)] =
true;
2312 at_boundary[line->vertex_index(1)] =
true;
2315 minimal_length[line->vertex_index(0)] =
2317 minimal_length[line->vertex_index(0)]);
2318 minimal_length[line->vertex_index(1)] =
2320 minimal_length[line->vertex_index(1)]);
2326 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
2327 if (cell->at_boundary(vertex) ==
true)
2328 at_boundary[cell->vertex_index(vertex)] =
true;
2330 minimal_length[cell->vertex_index(0)] =
2332 minimal_length[cell->vertex_index(0)]);
2333 minimal_length[cell->vertex_index(1)] =
2335 minimal_length[cell->vertex_index(1)]);
2340 boost::random::mt19937 rng(seed);
2341 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
2345 if (
auto distributed_triangulation =
2349 const std::vector<bool> locally_owned_vertices =
2351 std::vector<bool> vertex_moved(
triangulation.n_vertices(),
false);
2354 for (
const auto &cell :
triangulation.active_cell_iterators())
2355 if (cell->is_locally_owned())
2357 for (
const unsigned int vertex_no : cell->vertex_indices())
2359 const unsigned global_vertex_no =
2360 cell->vertex_index(vertex_no);
2365 if ((keep_boundary && at_boundary[global_vertex_no]) ||
2366 vertex_moved[global_vertex_no] ||
2367 !locally_owned_vertices[global_vertex_no])
2372 for (
unsigned int d = 0;
d < spacedim; ++
d)
2373 shift_vector(
d) = uniform_distribution(rng);
2375 shift_vector *= factor * minimal_length[global_vertex_no] /
2376 std::sqrt(shift_vector.
square());
2379 cell->vertex(vertex_no) += shift_vector;
2380 vertex_moved[global_vertex_no] =
true;
2384 distributed_triangulation->communicate_locally_moved_vertices(
2385 locally_owned_vertices);
2395 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
2396 const std::vector<Point<spacedim>> &old_vertex_locations =
2399 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2403 if (keep_boundary && at_boundary[vertex])
2404 new_vertex_locations[vertex] = old_vertex_locations[vertex];
2409 for (
unsigned int d = 0;
d < spacedim; ++
d)
2410 shift_vector(
d) = uniform_distribution(rng);
2412 shift_vector *= factor * minimal_length[vertex] /
2413 std::sqrt(shift_vector.
square());
2416 new_vertex_locations[vertex] =
2417 old_vertex_locations[vertex] + shift_vector;
2422 for (
const auto &cell :
triangulation.active_cell_iterators())
2423 for (
const unsigned int vertex_no : cell->vertex_indices())
2424 cell->vertex(vertex_no) =
2425 new_vertex_locations[cell->vertex_index(vertex_no)];
2441 for (; cell != endc; ++cell)
2442 if (!cell->is_artificial())
2443 for (
const unsigned int face : cell->face_indices())
2444 if (cell->face(face)->has_children() &&
2445 !cell->face(face)->at_boundary())
2449 cell->face(face)->child(0)->vertex(1) =
2450 (cell->face(face)->vertex(0) +
2451 cell->face(face)->vertex(1)) /
2455 cell->face(face)->child(0)->vertex(1) =
2456 .5 * (cell->face(face)->vertex(0) +
2457 cell->face(face)->vertex(1));
2458 cell->face(face)->child(0)->vertex(2) =
2459 .5 * (cell->face(face)->vertex(0) +
2460 cell->face(face)->vertex(2));
2461 cell->face(face)->child(1)->vertex(3) =
2462 .5 * (cell->face(face)->vertex(1) +
2463 cell->face(face)->vertex(3));
2464 cell->face(face)->child(2)->vertex(3) =
2465 .5 * (cell->face(face)->vertex(2) +
2466 cell->face(face)->vertex(3));
2469 cell->face(face)->child(0)->vertex(3) =
2470 .25 * (cell->face(face)->vertex(0) +
2471 cell->face(face)->vertex(1) +
2472 cell->face(face)->vertex(2) +
2473 cell->face(face)->vertex(3));
2481 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2485 const std::vector<bool> & marked_vertices)
2494 marked_vertices.size() == 0,
2496 marked_vertices.size()));
2503 marked_vertices.size() == 0 ||
2504 std::equal(marked_vertices.begin(),
2505 marked_vertices.end(),
2507 [](
bool p,
bool q) { return !p || q; }),
2509 "marked_vertices should be a subset of used vertices in the triangulation "
2510 "but marked_vertices contains one or more vertices that are not used vertices!"));
2514 const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
2520 std::vector<bool>::const_iterator
first =
2521 std::find(vertices_to_use.begin(), vertices_to_use.end(),
true);
2526 unsigned int best_vertex = std::distance(vertices_to_use.begin(),
first);
2527 double best_dist = (p -
vertices[best_vertex]).norm_square();
2531 for (
unsigned int j = best_vertex + 1; j <
vertices.size(); ++j)
2532 if (vertices_to_use[j])
2534 const double dist = (p -
vertices[j]).norm_square();
2535 if (dist < best_dist)
2547 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2550 const MeshType<dim, spacedim> &mesh,
2552 const std::vector<bool> & marked_vertices)
2565 marked_vertices.size() == 0,
2567 marked_vertices.size()));
2575 marked_vertices.size() == 0 ||
2576 std::equal(marked_vertices.begin(),
2577 marked_vertices.end(),
2579 [](
bool p,
bool q) { return !p || q; }),
2581 "marked_vertices should be a subset of used vertices in the triangulation "
2582 "but marked_vertices contains one or more vertices that are not used vertices!"));
2585 if (marked_vertices.size() != 0)
2588 if (marked_vertices[it->first] ==
false)
2603 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2605 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
2608 typename ::internal::
2609 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2612 const unsigned int vertex)
2618 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
2625 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
2661 for (
const auto &cell : mesh.active_cell_iterators())
2663 for (
const unsigned int v : cell->vertex_indices())
2664 if (cell->vertex_index(v) == vertex)
2676 for (
const auto face :
2677 cell->reference_cell().faces_for_given_vertex(v))
2678 if (!cell->at_boundary(face) &&
2679 cell->neighbor(face)->is_active())
2701 for (
unsigned int e = 0;
e < cell->n_lines(); ++
e)
2702 if (cell->line(
e)->has_children())
2706 if (cell->line(
e)->child(0)->vertex_index(1) == vertex)
2730 typename ::internal::
2731 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
2737 template <
int dim,
int spacedim>
2738 std::vector<std::vector<Tensor<1, spacedim>>>
2746 const unsigned int n_vertices = vertex_to_cells.size();
2751 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
2753 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
2756 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
2757 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
2759 typename std::set<typename Triangulation<dim, spacedim>::
2760 active_cell_iterator>::iterator it =
2761 vertex_to_cells[vertex].begin();
2762 for (
unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
2764 vertex_to_cell_centers[vertex][cell] =
2765 (*it)->center() -
vertices[vertex];
2766 vertex_to_cell_centers[vertex][cell] /=
2767 vertex_to_cell_centers[vertex][cell].
norm();
2770 return vertex_to_cell_centers;
2776 template <
int spacedim>
2779 const unsigned int a,
2780 const unsigned int b,
2784 const double scalar_product_a = center_directions[a] * point_direction;
2785 const double scalar_product_b = center_directions[
b] * point_direction;
2790 return (scalar_product_a > scalar_product_b);
2794 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
2796 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
2798 std::pair<typename ::internal::
2799 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
2804 const MeshType<dim, spacedim> &mesh,
2807 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
2810 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
2811 const std::vector<bool> & marked_vertices,
2813 const double tolerance,
2817 *relevant_cell_bounding_boxes_rtree)
2819 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2822 cell_and_position.first = mesh.end();
2826 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
2828 cell_and_position_approx;
2830 if (relevant_cell_bounding_boxes_rtree !=
nullptr &&
2831 !relevant_cell_bounding_boxes_rtree->empty())
2837 for (
unsigned int d = 0;
d < spacedim; ++
d)
2839 p1[
d] = p1[
d] - tolerance;
2840 p2[
d] = p2[
d] + tolerance;
2845 if (relevant_cell_bounding_boxes_rtree->qbegin(
2846 boost::geometry::index::intersects(bb)) ==
2847 relevant_cell_bounding_boxes_rtree->qend())
2848 return cell_and_position;
2851 bool found_cell =
false;
2852 bool approx_cell =
false;
2854 unsigned int closest_vertex_index = 0;
2856 auto current_cell = cell_hint;
2858 while (found_cell ==
false)
2864 const auto cell_vertices = mapping.
get_vertices(current_cell);
2865 const unsigned int closest_vertex =
2866 find_closest_vertex_of_cell<dim, spacedim>(current_cell,
2869 vertex_to_point = p - cell_vertices[closest_vertex];
2870 closest_vertex_index = current_cell->vertex_index(closest_vertex);
2874 if (!used_vertices_rtree.empty())
2877 using ValueType = std::pair<Point<spacedim>,
unsigned int>;
2878 std::function<
bool(
const ValueType &)> marked;
2879 if (marked_vertices.size() == mesh.n_vertices())
2880 marked = [&marked_vertices](
const ValueType &value) ->
bool {
2881 return marked_vertices[value.second];
2884 marked = [](
const ValueType &) ->
bool {
return true; };
2886 std::vector<std::pair<Point<spacedim>,
unsigned int>> res;
2887 used_vertices_rtree.query(
2888 boost::geometry::index::nearest(p, 1) &&
2889 boost::geometry::index::satisfies(marked),
2890 std::back_inserter(res));
2894 closest_vertex_index = res[0].second;
2899 mapping, mesh, p, marked_vertices);
2901 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
2904 const double vertex_point_norm = vertex_to_point.
norm();
2905 if (vertex_point_norm > 0)
2906 vertex_to_point /= vertex_point_norm;
2908 const unsigned int n_neighbor_cells =
2909 vertex_to_cells[closest_vertex_index].size();
2912 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
2914 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
2915 neighbor_permutation[i] = i;
2917 auto comp = [&](
const unsigned int a,
const unsigned int b) ->
bool {
2918 return internal::compare_point_association<spacedim>(
2922 vertex_to_cell_centers[closest_vertex_index]);
2925 std::sort(neighbor_permutation.begin(),
2926 neighbor_permutation.end(),
2931 double best_distance = tolerance;
2935 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
2939 auto cell = vertex_to_cells[closest_vertex_index].begin();
2942 if (!(*cell)->is_artificial())
2949 cell_and_position.first = *cell;
2950 cell_and_position.second = p_unit;
2952 approx_cell =
false;
2961 if (dist < best_distance)
2963 best_distance = dist;
2964 cell_and_position_approx.first = *cell;
2965 cell_and_position_approx.second = p_unit;
2974 if (found_cell ==
true)
2975 return cell_and_position;
2976 else if (approx_cell ==
true)
2977 return cell_and_position_approx;
2993 mapping, mesh, p, marked_vertices, tolerance);
2995 current_cell =
typename MeshType<dim, spacedim>::active_cell_iterator();
2997 return cell_and_position;
3002 template <
int dim,
int spacedim>
3011 unsigned int closest_vertex = 0;
3013 for (
unsigned int v = 1; v < cell->
n_vertices(); ++v)
3016 if (vertex_distance < minimum_distance)
3019 minimum_distance = vertex_distance;
3022 return closest_vertex;
3029 namespace BoundingBoxPredicate
3031 template <
class MeshType>
3032 std::tuple<BoundingBox<MeshType::space_dimension>,
bool>
3034 const typename MeshType::cell_iterator &parent_cell,
3035 const std::function<
3036 bool(
const typename MeshType::active_cell_iterator &)> &predicate)
3038 bool has_predicate =
3040 std::vector<typename MeshType::active_cell_iterator> active_cells;
3041 if (parent_cell->is_active())
3042 active_cells = {parent_cell};
3046 active_cells = get_active_child_cells<MeshType>(parent_cell);
3048 const unsigned int spacedim = MeshType::space_dimension;
3052 while (i < active_cells.size() && !predicate(active_cells[i]))
3056 if (active_cells.size() == 0 || i == active_cells.size())
3059 return std::make_tuple(bbox, has_predicate);
3066 for (; i < active_cells.size(); ++i)
3067 if (predicate(active_cells[i]))
3069 for (
unsigned int d = 0;
d < spacedim; ++
d)
3071 minp[
d] =
std::min(minp[
d], active_cells[i]->vertex(v)[
d]);
3072 maxp[
d] =
std::max(maxp[
d], active_cells[i]->vertex(v)[
d]);
3075 has_predicate =
true;
3077 return std::make_tuple(bbox, has_predicate);
3084 template <
class MeshType>
3085 std::vector<BoundingBox<MeshType::space_dimension>>
3087 const MeshType &mesh,
3088 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
3090 const unsigned int refinement_level,
3091 const bool allow_merge,
3092 const unsigned int max_boxes)
3099 refinement_level <= mesh.n_levels(),
3101 "Error: refinement level is higher then total levels in the triangulation!"));
3103 const unsigned int spacedim = MeshType::space_dimension;
3104 std::vector<BoundingBox<spacedim>> bounding_boxes;
3108 for (
unsigned int i = 0; i < refinement_level; ++i)
3109 for (
const typename MeshType::cell_iterator &cell :
3110 mesh.active_cell_iterators_on_level(i))
3112 bool has_predicate =
false;
3114 std::tie(bbox, has_predicate) =
3116 MeshType>(cell, predicate);
3118 bounding_boxes.push_back(bbox);
3122 for (
const typename MeshType::cell_iterator &cell :
3123 mesh.cell_iterators_on_level(refinement_level))
3125 bool has_predicate =
false;
3127 std::tie(bbox, has_predicate) =
3129 MeshType>(cell, predicate);
3131 bounding_boxes.push_back(bbox);
3136 return bounding_boxes;
3142 std::vector<unsigned int> merged_boxes_idx;
3143 bool found_neighbors =
true;
3148 while (found_neighbors)
3150 found_neighbors =
false;
3151 for (
unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
3153 if (std::find(merged_boxes_idx.begin(),
3154 merged_boxes_idx.end(),
3155 i) == merged_boxes_idx.end())
3156 for (
unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
3157 if (std::find(merged_boxes_idx.begin(),
3158 merged_boxes_idx.end(),
3159 j) == merged_boxes_idx.end() &&
3160 bounding_boxes[i].get_neighbor_type(
3161 bounding_boxes[j]) ==
3164 bounding_boxes[i].merge_with(bounding_boxes[j]);
3165 merged_boxes_idx.push_back(j);
3166 found_neighbors =
true;
3172 std::vector<BoundingBox<spacedim>> merged_b_boxes;
3173 for (
unsigned int i = 0; i < bounding_boxes.size(); ++i)
3174 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
3175 merged_boxes_idx.end())
3176 merged_b_boxes.push_back(bounding_boxes[i]);
3181 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
3183 std::vector<double> volumes;
3184 for (
unsigned int i = 0; i < merged_b_boxes.size(); ++i)
3185 volumes.push_back(merged_b_boxes[i].volume());
3187 while (merged_b_boxes.size() > max_boxes)
3189 unsigned int min_idx =
3190 std::min_element(volumes.begin(), volumes.end()) -
3192 volumes.erase(volumes.begin() + min_idx);
3194 bool not_removed =
true;
3195 for (
unsigned int i = 0;
3196 i < merged_b_boxes.size() && not_removed;
3201 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
3202 merged_b_boxes[min_idx]) ==
3204 merged_b_boxes[i].get_neighbor_type(
3205 merged_b_boxes[min_idx]) ==
3208 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
3209 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
3210 not_removed =
false;
3213 ExcMessage(
"Error: couldn't merge bounding boxes!"));
3216 Assert(merged_b_boxes.size() <= max_boxes,
3218 "Error: couldn't reach target number of bounding boxes!"));
3219 return merged_b_boxes;
3225 template <
int spacedim>
3227 std::tuple<std::vector<std::vector<unsigned int>>,
3228 std::map<unsigned int, unsigned int>,
3229 std::map<unsigned int, std::vector<unsigned int>>>
3237 unsigned int n_procs = global_bboxes.size();
3238 std::vector<std::vector<unsigned int>> point_owners(n_procs);
3239 std::map<unsigned int, unsigned int> map_owners_found;
3240 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3242 unsigned int n_points = points.size();
3243 for (
unsigned int pt = 0; pt < n_points; ++pt)
3246 std::vector<unsigned int> owners_found;
3248 for (
unsigned int rk = 0; rk < n_procs; ++rk)
3251 if (bbox.point_inside(points[pt]))
3253 point_owners[rk].emplace_back(pt);
3254 owners_found.emplace_back(rk);
3258 Assert(owners_found.size() > 0,
3259 ExcMessage(
"No owners found for the point " +
3261 if (owners_found.size() == 1)
3262 map_owners_found[pt] = owners_found[0];
3265 map_owners_guessed[pt] = owners_found;
3268 return std::make_tuple(std::move(point_owners),
3269 std::move(map_owners_found),
3270 std::move(map_owners_guessed));
3273 template <
int spacedim>
3275 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
3276 std::map<unsigned int, unsigned int>,
3277 std::map<unsigned int, std::vector<unsigned int>>>
3285 std::map<unsigned int, std::vector<unsigned int>> point_owners;
3286 std::map<unsigned int, unsigned int> map_owners_found;
3287 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
3288 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> search_result;
3290 unsigned int n_points = points.size();
3291 for (
unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
3293 search_result.clear();
3296 covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
3297 std::back_inserter(search_result));
3300 std::set<unsigned int> owners_found;
3302 for (
const auto &rank_bbox : search_result)
3306 const bool pt_inserted = owners_found.insert(pt_n).second;
3308 point_owners[rank_bbox.second].emplace_back(pt_n);
3310 Assert(owners_found.size() > 0,
3311 ExcMessage(
"No owners found for the point " +
3313 if (owners_found.size() == 1)
3314 map_owners_found[pt_n] = *owners_found.begin();
3319 std::back_inserter(map_owners_guessed[pt_n]));
3322 return std::make_tuple(std::move(point_owners),
3323 std::move(map_owners_found),
3324 std::move(map_owners_guessed));
3328 template <
int dim,
int spacedim>
3330 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3334 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
3339 for (; cell != endc; ++cell)
3340 for (
const unsigned int i : cell->vertex_indices())
3345 for (; cell != endc; ++cell)
3347 for (
unsigned int i : cell->face_indices())
3349 if ((cell->at_boundary(i) ==
false) &&
3350 (cell->neighbor(i)->is_active()))
3353 adjacent_cell = cell->neighbor(i);
3354 for (
unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
3363 for (
unsigned int i = 0; i < cell->
n_lines(); ++i)
3364 if (cell->line(i)->has_children())
3378 template <
int dim,
int spacedim>
3379 std::map<unsigned int, types::global_vertex_index>
3383 std::map<unsigned int, types::global_vertex_index>
3384 local_to_global_vertex_index;
3386 #ifndef DEAL_II_WITH_MPI
3393 ExcMessage(
"This function does not make any sense "
3394 "for parallel::distributed::Triangulation "
3395 "objects if you do not have MPI enabled."));
3399 using active_cell_iterator =
3401 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
3406 unsigned int max_cellid_size = 0;
3407 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
3409 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
3417 std::set<active_cell_iterator> missing_vert_cells;
3418 std::set<unsigned int> used_vertex_index;
3419 for (; cell != endc; ++cell)
3421 if (cell->is_locally_owned())
3423 for (
const unsigned int i : cell->vertex_indices())
3426 typename std::set<active_cell_iterator>::iterator
3427 adjacent_cell = vertex_to_cell[cell->vertex_index(i)].begin(),
3428 end_adj_cell = vertex_to_cell[cell->vertex_index(i)].end();
3429 for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
3430 lowest_subdomain_id =
3432 (*adjacent_cell)->subdomain_id());
3435 if (lowest_subdomain_id == cell->subdomain_id())
3439 if (used_vertex_index.find(cell->vertex_index(i)) ==
3440 used_vertex_index.end())
3443 local_to_global_vertex_index[cell->vertex_index(i)] =
3449 vertex_to_cell[cell->vertex_index(i)].begin();
3450 for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
3451 if ((*adjacent_cell)->subdomain_id() !=
3452 cell->subdomain_id())
3456 tmp((*adjacent_cell)->subdomain_id(),
3457 cell->vertex_index(i));
3458 if (vertices_added.find(tmp) ==
3459 vertices_added.end())
3461 vertices_to_send[(*adjacent_cell)
3464 cell->vertex_index(i),
3465 cell->id().to_string());
3466 if (cell->id().to_string().size() >
3469 cell->id().to_string().size();
3470 vertices_added.insert(tmp);
3473 used_vertex_index.insert(cell->vertex_index(i));
3480 vertices_to_recv[lowest_subdomain_id].insert(
3481 cell->vertex_index(i));
3482 missing_vert_cells.insert(cell);
3489 if (cell->is_ghost())
3491 for (
unsigned int i : cell->face_indices())
3493 if (cell->at_boundary(i) ==
false)
3495 if (cell->neighbor(i)->is_active())
3498 spacedim>::active_cell_iterator
3499 adjacent_cell = cell->neighbor(i);
3500 if ((adjacent_cell->is_locally_owned()))
3503 adjacent_cell->subdomain_id();
3504 if (cell->subdomain_id() < adj_subdomain_id)
3505 for (
unsigned int j = 0;
3506 j < cell->face(i)->n_vertices();
3509 vertices_to_recv[cell->subdomain_id()].insert(
3510 cell->face(i)->vertex_index(j));
3511 missing_vert_cells.insert(cell);
3527 int ierr = MPI_Exscan(&next_index,
3535 std::map<unsigned int, types::global_vertex_index>::iterator
3536 global_index_it = local_to_global_vertex_index.begin(),
3537 global_index_end = local_to_global_vertex_index.end();
3538 for (; global_index_it != global_index_end; ++global_index_it)
3539 global_index_it->second +=
shift;
3553 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
3554 vertices_to_send.size());
3555 std::vector<MPI_Request> first_requests(vertices_to_send.size());
3559 std::string>>>::iterator
3560 vert_to_send_it = vertices_to_send.begin(),
3561 vert_to_send_end = vertices_to_send.end();
3562 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3563 ++vert_to_send_it, ++i)
3565 int destination = vert_to_send_it->first;
3566 const unsigned int n_vertices = vert_to_send_it->second.size();
3567 const int buffer_size = 2 * n_vertices;
3568 vertices_send_buffers[i].resize(buffer_size);
3571 for (
unsigned int j = 0; j < n_vertices; ++j)
3573 vertices_send_buffers[i][2 * j] =
3574 std::get<0>(vert_to_send_it->second[j]);
3575 vertices_send_buffers[i][2 * j + 1] =
3576 local_to_global_vertex_index[std::get<1>(
3577 vert_to_send_it->second[j])];
3581 ierr = MPI_Isend(vertices_send_buffers[i].data(),
3587 &first_requests[i]);
3592 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
3593 vertices_to_recv.size());
3594 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
3595 vert_to_recv_it = vertices_to_recv.begin(),
3596 vert_to_recv_end = vertices_to_recv.end();
3597 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3598 ++vert_to_recv_it, ++i)
3600 int source = vert_to_recv_it->first;
3601 const unsigned int n_vertices = vert_to_recv_it->second.size();
3602 const int buffer_size = 2 * n_vertices;
3603 vertices_recv_buffers[i].resize(buffer_size);
3606 ierr = MPI_Recv(vertices_recv_buffers[i].data(),
3618 std::vector<std::vector<char>> cellids_send_buffers(
3619 vertices_to_send.size());
3620 std::vector<MPI_Request> second_requests(vertices_to_send.size());
3621 vert_to_send_it = vertices_to_send.begin();
3622 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
3623 ++vert_to_send_it, ++i)
3625 int destination = vert_to_send_it->first;
3626 const unsigned int n_vertices = vert_to_send_it->second.size();
3627 const int buffer_size = max_cellid_size * n_vertices;
3628 cellids_send_buffers[i].resize(buffer_size);
3631 unsigned int pos = 0;
3632 for (
unsigned int j = 0; j < n_vertices; ++j)
3634 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
3635 for (
unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
3637 if (k < cell_id.size())
3638 cellids_send_buffers[i][pos] = cell_id[k];
3642 cellids_send_buffers[i][pos] =
'-';
3647 ierr = MPI_Isend(cellids_send_buffers[i].data(),
3653 &second_requests[i]);
3658 std::vector<std::vector<char>> cellids_recv_buffers(
3659 vertices_to_recv.size());
3660 vert_to_recv_it = vertices_to_recv.begin();
3661 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3662 ++vert_to_recv_it, ++i)
3664 int source = vert_to_recv_it->first;
3665 const unsigned int n_vertices = vert_to_recv_it->second.size();
3666 const int buffer_size = max_cellid_size * n_vertices;
3667 cellids_recv_buffers[i].resize(buffer_size);
3670 ierr = MPI_Recv(cellids_recv_buffers[i].data(),
3682 vert_to_recv_it = vertices_to_recv.begin();
3683 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
3684 ++i, ++vert_to_recv_it)
3686 for (
unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
3688 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
3690 vertices_recv_buffers[i][2 * j + 1];
3691 const std::string cellid_recv(
3692 &cellids_recv_buffers[i][max_cellid_size * j],
3693 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
3695 typename std::set<active_cell_iterator>::iterator
3696 cell_set_it = missing_vert_cells.begin(),
3697 end_cell_set = missing_vert_cells.end();
3698 for (; (found ==
false) && (cell_set_it != end_cell_set);
3701 typename std::set<active_cell_iterator>::iterator
3703 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
3705 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
3706 for (; candidate_cell != end_cell; ++candidate_cell)
3708 std::string current_cellid =
3709 (*candidate_cell)->id().to_string();
3710 current_cellid.resize(max_cellid_size,
'-');
3711 if (current_cellid.compare(cellid_recv) == 0)
3713 local_to_global_vertex_index
3714 [(*candidate_cell)->vertex_index(local_pos_recv)] =
3726 return local_to_global_vertex_index;
3731 template <
int dim,
int spacedim>
3747 for (
const auto &cell :
triangulation.active_cell_iterators())
3749 const unsigned int index = cell->active_cell_index();
3750 cell_connectivity.
add(index, index);
3751 for (
auto f : cell->face_indices())
3752 if ((cell->at_boundary(f) ==
false) &&
3753 (cell->neighbor(f)->has_children() ==
false))
3755 const unsigned int other_index =
3756 cell->neighbor(f)->active_cell_index();
3757 cell_connectivity.
add(index, other_index);
3758 cell_connectivity.
add(other_index, index);
3765 template <
int dim,
int spacedim>
3771 std::vector<std::vector<unsigned int>> vertex_to_cell(
3773 for (
const auto &cell :
triangulation.active_cell_iterators())
3775 for (
const unsigned int v : cell->vertex_indices())
3776 vertex_to_cell[cell->vertex_index(v)].push_back(
3777 cell->active_cell_index());
3782 for (
const auto &cell :
triangulation.active_cell_iterators())
3784 for (
const unsigned int v : cell->vertex_indices())
3785 for (
unsigned int n = 0;
3786 n < vertex_to_cell[cell->vertex_index(v)].size();
3788 cell_connectivity.
add(cell->active_cell_index(),
3789 vertex_to_cell[cell->vertex_index(v)][n]);
3794 template <
int dim,
int spacedim>
3798 const unsigned int level,
3801 std::vector<std::vector<unsigned int>> vertex_to_cell(
3808 for (
const unsigned int v : cell->vertex_indices())
3809 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
3819 for (
const unsigned int v : cell->vertex_indices())
3820 for (
unsigned int n = 0;
3821 n < vertex_to_cell[cell->vertex_index(v)].size();
3823 cell_connectivity.
add(cell->index(),
3824 vertex_to_cell[cell->vertex_index(v)][n]);
3830 template <
int dim,
int spacedim>
3838 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3839 "are already partitioned implicitly and can not be "
3840 "partitioned again explicitly."));
3842 std::vector<unsigned int> cell_weights;
3852 for (
const auto &cell :
triangulation.active_cell_iterators())
3853 if (cell->is_locally_owned())
3854 cell_weights[cell->active_cell_index()] =
3864 if (
const auto shared_tria =
3868 shared_tria->get_communicator(),
3872 Assert(std::accumulate(cell_weights.begin(),
3874 std::uint64_t(0)) > 0,
3875 ExcMessage(
"The global sum of weights over all active cells "
3876 "is zero. Please verify how you generate weights."));
3888 template <
int dim,
int spacedim>
3891 const std::vector<unsigned int> &cell_weights,
3897 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3898 "are already partitioned implicitly and can not be "
3899 "partitioned again explicitly."));
3903 if (n_partitions == 1)
3905 for (
const auto &cell :
triangulation.active_cell_iterators())
3906 cell->set_subdomain_id(0);
3920 sp_cell_connectivity.
copy_from(cell_connectivity);
3923 sp_cell_connectivity,
3930 template <
int dim,
int spacedim>
3939 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
3940 "are already partitioned implicitly and can not be "
3941 "partitioned again explicitly."));
3943 std::vector<unsigned int> cell_weights;
3953 for (
const auto &cell :
triangulation.active_cell_iterators() |
3955 cell_weights[cell->active_cell_index()] =
3965 if (
const auto shared_tria =
3969 shared_tria->get_communicator(),
3973 Assert(std::accumulate(cell_weights.begin(),
3975 std::uint64_t(0)) > 0,
3976 ExcMessage(
"The global sum of weights over all active cells "
3977 "is zero. Please verify how you generate weights."));
3983 cell_connection_graph,
3990 template <
int dim,
int spacedim>
3993 const std::vector<unsigned int> &cell_weights,
4000 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4001 "are already partitioned implicitly and can not be "
4002 "partitioned again explicitly."));
4005 ExcMessage(
"Connectivity graph has wrong size"));
4007 ExcMessage(
"Connectivity graph has wrong size"));
4013 if (n_partitions == 1)
4015 for (
const auto &cell :
triangulation.active_cell_iterators())
4016 cell->set_subdomain_id(0);
4024 std::vector<unsigned int> partition_indices(
triangulation.n_active_cells());
4032 for (
const auto &cell :
triangulation.active_cell_iterators())
4033 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
4045 unsigned int & current_proc_idx,
4046 unsigned int & current_cell_idx,
4048 const unsigned int n_partitions)
4050 if (cell->is_active())
4052 while (current_cell_idx >=
4054 (current_proc_idx + 1) / n_partitions))
4056 cell->set_subdomain_id(current_proc_idx);
4061 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4071 template <
int dim,
int spacedim>
4075 const bool group_siblings)
4079 ExcMessage(
"Objects of type parallel::distributed::Triangulation "
4080 "are already partitioned implicitly and can not be "
4081 "partitioned again explicitly."));
4089 if (n_partitions == 1)
4091 for (
const auto &cell :
triangulation.active_cell_iterators())
4092 cell->set_subdomain_id(0);
4098 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
4099 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
4105 coarse_cell_to_p4est_tree_permutation.resize(
triangulation.n_cells(0));
4107 coarse_cell_to_p4est_tree_permutation);
4109 p4est_tree_to_coarse_cell_permutation =
4112 unsigned int current_proc_idx = 0;
4113 unsigned int current_cell_idx = 0;
4118 for (
unsigned int idx = 0; idx <
triangulation.n_cells(0); ++idx)
4120 const unsigned int coarse_cell_idx =
4121 p4est_tree_to_coarse_cell_permutation[idx];
4144 for (; cell != endc; ++cell)
4146 if (cell->is_active())
4148 bool all_children_active =
true;
4149 std::map<unsigned int, unsigned int> map_cpu_n_cells;
4150 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4151 if (!cell->child(n)->is_active())
4153 all_children_active =
false;
4157 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
4159 if (!all_children_active)
4162 unsigned int new_owner = cell->child(0)->subdomain_id();
4163 for (std::map<unsigned int, unsigned int>::iterator it =
4164 map_cpu_n_cells.begin();
4165 it != map_cpu_n_cells.end();
4167 if (it->second > map_cpu_n_cells[new_owner])
4168 new_owner = it->first;
4170 for (
unsigned int n = 0; n < cell->n_children(); ++n)
4171 cell->child(n)->set_subdomain_id(new_owner);
4177 template <
int dim,
int spacedim>
4182 for (
int lvl = n_levels - 1; lvl >= 0; --lvl)
4184 for (
const auto &cell :
triangulation.cell_iterators_on_level(lvl))
4186 if (cell->is_active())
4187 cell->set_level_subdomain_id(cell->subdomain_id());
4190 Assert(cell->child(0)->level_subdomain_id() !=
4193 cell->set_level_subdomain_id(
4194 cell->child(0)->level_subdomain_id());
4206 template <
int dim,
int spacedim>
4211 const std::vector<CellId> & cell_ids,
4212 std::vector<types::subdomain_id> &subdomain_ids)
4214 #ifndef DEAL_II_WITH_P4EST
4217 (void)subdomain_ids;
4221 "You are attempting to use a functionality that is only available "
4222 "if deal.II was configured to use p4est, but cmake did not find a "
4223 "valid p4est library."));
4228 for (
const auto &cell_id : cell_ids)
4231 typename ::internal::p4est::types<dim>::quadrant p4est_cell,
4234 ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
4235 for (
const auto &child_index : cell_id.get_child_indices())
4237 ::internal::p4est::init_quadrant_children<dim>(
4238 p4est_cell, p4est_children);
4240 p4est_children[
static_cast<unsigned int>(child_index)];
4246 const_cast<typename ::internal::p4est::types<dim>::forest
4248 cell_id.get_coarse_cell_id(),
4255 subdomain_ids.push_back(owner);
4262 template <
int spacedim>
4266 const std::vector<CellId> &,
4267 std::vector<types::subdomain_id> &)
4276 template <
int dim,
int spacedim>
4277 std::vector<types::subdomain_id>
4279 const std::vector<CellId> & cell_ids)
4281 std::vector<types::subdomain_id> subdomain_ids;
4282 subdomain_ids.reserve(cell_ids.size());
4291 *parallel_tria =
dynamic_cast<
4305 const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
4306 shared_tria->get_true_subdomain_ids_of_cells();
4308 for (
const auto &cell_id : cell_ids)
4310 const unsigned int active_cell_index =
4311 shared_tria->create_cell_iterator(cell_id)->active_cell_index();
4312 subdomain_ids.push_back(
4313 true_subdomain_ids_of_cells[active_cell_index]);
4320 for (
const auto &cell_id : cell_ids)
4322 subdomain_ids.push_back(
4323 triangulation.create_cell_iterator(cell_id)->subdomain_id());
4327 return subdomain_ids;
4332 template <
int dim,
int spacedim>
4335 std::vector<types::subdomain_id> & subdomain)
4340 for (
const auto &cell :
triangulation.active_cell_iterators())
4341 subdomain[cell->active_cell_index()] = cell->subdomain_id();
4346 template <
int dim,
int spacedim>
4352 unsigned int count = 0;
4353 for (
const auto &cell :
triangulation.active_cell_iterators())
4354 if (cell->subdomain_id() == subdomain)
4362 template <
int dim,
int spacedim>
4367 std::vector<bool> locally_owned_vertices =
4374 if (
const auto *tr =
dynamic_cast<
4377 for (
const auto &cell :
triangulation.active_cell_iterators())
4378 if (cell->is_artificial() ||
4379 (cell->is_ghost() &&
4380 (cell->subdomain_id() < tr->locally_owned_subdomain())))
4381 for (
const unsigned int v : cell->vertex_indices())
4382 locally_owned_vertices[cell->vertex_index(v)] =
false;
4384 return locally_owned_vertices;
4389 template <
int dim,
int spacedim>
4395 for (
const auto &cell :
triangulation.active_cell_iterators())
4396 if (!cell->is_artificial())
4397 min_diameter =
std::min(min_diameter, cell->diameter(mapping));
4399 double global_min_diameter = 0;
4401 #ifdef DEAL_II_WITH_MPI
4405 global_min_diameter =
4409 global_min_diameter = min_diameter;
4411 return global_min_diameter;
4416 template <
int dim,
int spacedim>
4421 double max_diameter = 0.;
4422 for (
const auto &cell :
triangulation.active_cell_iterators())
4423 if (!cell->is_artificial())
4424 max_diameter =
std::max(max_diameter, cell->diameter(mapping));
4426 double global_max_diameter = 0;
4428 #ifdef DEAL_II_WITH_MPI
4432 global_max_diameter =
4436 global_max_diameter = max_diameter;
4438 return global_max_diameter;
4445 namespace FixUpDistortedChildCells
4468 template <
typename Iterator,
int spacedim>
4473 const unsigned int structdim =
4474 Iterator::AccessorType::structure_dimension;
4475 Assert(spacedim == Iterator::AccessorType::dimension,
4481 Assert(object->refinement_case() ==
4489 Tensor<spacedim - structdim, spacedim>
4492 for (
const unsigned int i : object->vertex_indices())
4493 parent_vertices[i] =
object->vertex(i);
4496 parent_vertices, parent_alternating_forms);
4498 const Tensor<spacedim - structdim, spacedim>
4499 average_parent_alternating_form =
4500 std::accumulate(parent_alternating_forms,
4501 parent_alternating_forms +
4515 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4519 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4520 for (
const unsigned int i : object->child(c)->vertex_indices())
4521 child_vertices[c][i] =
object->child(c)->vertex(i);
4529 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4531 1] = object_mid_point;
4533 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4535 child_vertices[c], child_alternating_forms[c]);
4547 double objective = 0;
4548 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4549 for (
const unsigned int i : object->child(c)->vertex_indices())
4551 (child_alternating_forms[c][i] -
4552 average_parent_alternating_form / std::pow(2., 1. * structdim))
4564 template <
typename Iterator>
4567 const unsigned int f,
4568 std::integral_constant<int, 1>)
4570 return object->vertex(f);
4580 template <
typename Iterator>
4583 const unsigned int f,
4584 std::integral_constant<int, 2>)
4586 return object->line(f)->center();
4596 template <
typename Iterator>
4599 const unsigned int f,
4600 std::integral_constant<int, 3>)
4602 return object->face(f)->center();
4629 template <
typename Iterator>
4633 const unsigned int structdim =
4634 Iterator::AccessorType::structure_dimension;
4636 double diameter =
object->diameter();
4637 for (
const unsigned int f : object->face_indices())
4638 for (
unsigned int e = f + 1;
e <
object->n_faces(); ++
e)
4643 std::integral_constant<int, structdim>())
4645 object,
e, std::integral_constant<int, structdim>())));
4656 template <
typename Iterator>
4660 const unsigned int structdim =
4661 Iterator::AccessorType::structure_dimension;
4662 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
4669 Assert(object->refinement_case() ==
4679 unsigned int iteration = 0;
4684 double initial_delta = 0;
4691 const double step_length =
diameter / 4 / (iteration + 1);
4696 for (
unsigned int d = 0;
d < spacedim; ++
d)
4698 const double eps = step_length / 10;
4712 if (gradient.norm() == 0)
4721 std::min(2 * current_value / (gradient * gradient),
4722 step_length / gradient.norm()) *
4727 const double previous_value = current_value;
4731 initial_delta = (previous_value - current_value);
4734 if ((iteration >= 1) &&
4735 ((previous_value - current_value < 0) ||
4736 (
std::fabs(previous_value - current_value) <
4737 0.001 * initial_delta)))
4742 while (iteration < 20);
4758 double old_min_product, new_min_product;
4763 parent_vertices[i] = object->vertex(i);
4765 Tensor<spacedim - structdim, spacedim>
4768 parent_vertices, parent_alternating_forms);
4774 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4775 for (
const unsigned int i : object->child(c)->vertex_indices())
4776 child_vertices[c][i] =
object->child(c)->vertex(i);
4778 Tensor<spacedim - structdim, spacedim> child_alternating_forms
4782 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4784 child_vertices[c], child_alternating_forms[c]);
4787 child_alternating_forms[0][0] * parent_alternating_forms[0];
4788 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4789 for (
const unsigned int i : object->child(c)->vertex_indices())
4790 for (
const unsigned int j : object->vertex_indices())
4791 old_min_product = std::min<double>(old_min_product,
4792 child_alternating_forms[c][i] *
4793 parent_alternating_forms[j]);
4801 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4803 1] = object_mid_point;
4805 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4807 child_vertices[c], child_alternating_forms[c]);
4810 child_alternating_forms[0][0] * parent_alternating_forms[0];
4811 for (
unsigned int c = 0; c <
object->n_children(); ++c)
4812 for (
const unsigned int i : object->child(c)->vertex_indices())
4813 for (
const unsigned int j : object->vertex_indices())
4814 new_min_product = std::min<double>(new_min_product,
4815 child_alternating_forms[c][i] *
4816 parent_alternating_forms[j]);
4824 if (new_min_product >= old_min_product)
4825 object->child(0)->vertex(
4832 return (
std::max(new_min_product, old_min_product) > 0);
4838 template <
int dim,
int spacedim>
4841 const typename ::Triangulation<dim, spacedim>::cell_iterator
4843 std::integral_constant<int, dim>,
4844 std::integral_constant<int, spacedim>)
4854 for (
auto f : cell->face_indices())
4857 Assert(cell->face(f)->refinement_case() ==
4861 bool subface_is_more_refined =
false;
4862 for (
unsigned int g = 0;
4863 g < GeometryInfo<dim>::max_children_per_face;
4865 if (cell->face(f)->child(g)->has_children())
4867 subface_is_more_refined =
true;
4871 if (subface_is_more_refined ==
true)
4882 template <
int dim,
int spacedim>
4890 dim != 1 && spacedim != 1,
4891 "This function is only valid when dim != 1 or spacedim != 1.");
4895 for (
typename std::list<
4897 cell_ptr = distorted_cells.distorted_cells.
begin();
4898 cell_ptr != distorted_cells.distorted_cells.
end();
4904 Assert(!cell->is_active(),
4906 "This function is only valid for a list of cells that "
4907 "have children (i.e., no cell in the list may be active)."));
4911 std::integral_constant<int, dim>(),
4912 std::integral_constant<int, spacedim>());
4916 unfixable_subset.distorted_cells.push_back(cell);
4919 return unfixable_subset;
4924 template <
int dim,
int spacedim>
4927 const bool reset_boundary_ids)
4930 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
4931 auto m_it = dst_manifold_ids.begin();
4932 for (
const auto b : src_boundary_ids)
4937 const std::vector<types::boundary_id> reset_boundary_id =
4938 reset_boundary_ids ?
4939 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
4949 template <
int dim,
int spacedim>
4952 const std::vector<types::boundary_id> &src_boundary_ids,
4953 const std::vector<types::manifold_id> &dst_manifold_ids,
4955 const std::vector<types::boundary_id> &reset_boundary_ids_)
4958 const auto reset_boundary_ids =
4959 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
4968 for (
auto f : cell->face_indices())
4969 if (cell->face(f)->at_boundary())
4970 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
4972 const auto bid = cell->face(f)->line(
e)->boundary_id();
4973 const unsigned int ind = std::find(src_boundary_ids.begin(),
4974 src_boundary_ids.end(),
4976 src_boundary_ids.begin();
4977 if (ind < src_boundary_ids.size())
4978 cell->face(f)->line(
e)->set_manifold_id(
4979 dst_manifold_ids[ind]);
4984 for (
auto f : cell->face_indices())
4985 if (cell->face(f)->at_boundary())
4987 const auto bid = cell->face(f)->boundary_id();
4988 const unsigned int ind =
4989 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
4990 src_boundary_ids.begin();
4992 if (ind < src_boundary_ids.size())
4995 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
4997 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
5001 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
5003 const auto bid = cell->face(f)->line(
e)->boundary_id();
5004 const unsigned int ind = std::find(src_boundary_ids.begin(),
5005 src_boundary_ids.end(),
5007 src_boundary_ids.begin();
5008 if (ind < src_boundary_ids.size())
5009 cell->face(f)->line(
e)->set_boundary_id(
5010 reset_boundary_ids[ind]);
5016 template <
int dim,
int spacedim>
5019 const bool compute_face_ids)
5025 for (; cell != endc; ++cell)
5027 cell->set_manifold_id(cell->material_id());
5028 if (compute_face_ids ==
true)
5030 for (
auto f : cell->face_indices())
5032 if (cell->at_boundary(f) ==
false)
5033 cell->face(f)->set_manifold_id(
5035 cell->neighbor(f)->material_id()));
5037 cell->face(f)->set_manifold_id(cell->material_id());
5044 template <
int dim,
int spacedim>
5049 const std::set<types::manifold_id> &)> &disambiguation_function,
5050 bool overwrite_only_flat_manifold_ids)
5055 const unsigned int n_subobjects =
5059 std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
5060 std::vector<unsigned int> backup;
5064 unsigned next_index = 1;
5068 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
5070 if (cell->line(
l)->user_index() == 0)
5073 manifold_ids[next_index].insert(cell->manifold_id());
5074 cell->line(
l)->set_user_index(next_index++);
5077 manifold_ids[cell->line(
l)->user_index()].insert(
5078 cell->manifold_id());
5081 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
5083 if (cell->quad(
l)->user_index() == 0)
5086 manifold_ids[next_index].insert(cell->manifold_id());
5087 cell->quad(
l)->set_user_index(next_index++);
5090 manifold_ids[cell->quad(
l)->user_index()].insert(
5091 cell->manifold_id());
5097 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
5099 const auto id = cell->line(
l)->user_index();
5103 if (cell->line(
l)->manifold_id() ==
5105 overwrite_only_flat_manifold_ids ==
false)
5106 cell->line(
l)->set_manifold_id(
5107 disambiguation_function(manifold_ids[
id]));
5108 cell->line(
l)->set_user_index(0);
5112 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
5114 const auto id = cell->quad(
l)->user_index();
5118 if (cell->quad(
l)->manifold_id() ==
5120 overwrite_only_flat_manifold_ids ==
false)
5121 cell->quad(
l)->set_manifold_id(
5122 disambiguation_function(manifold_ids[
id]));
5123 cell->quad(
l)->set_user_index(0);
5132 template <
int dim,
int spacedim>
5133 std::pair<unsigned int, double>
5137 double max_ratio = 1;
5138 unsigned int index = 0;
5140 for (
unsigned int i = 0; i < dim; ++i)
5141 for (
unsigned int j = i + 1; j < dim; ++j)
5143 unsigned int ax = i % dim;
5144 unsigned int next_ax = j % dim;
5147 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
5149 if (ratio > max_ratio)
5154 else if (1.0 / ratio > max_ratio)
5156 max_ratio = 1.0 / ratio;
5160 return std::make_pair(index, max_ratio);
5164 template <
int dim,
int spacedim>
5167 const bool isotropic,
5168 const unsigned int max_iterations)
5170 unsigned int iter = 0;
5171 bool continue_refinement =
true;
5173 while (continue_refinement && (iter < max_iterations))
5177 continue_refinement =
false;
5180 for (
const unsigned int j : cell->face_indices())
5181 if (cell->at_boundary(j) ==
false &&
5182 cell->neighbor(j)->has_children())
5186 cell->set_refine_flag();
5187 continue_refinement =
true;
5190 continue_refinement |= cell->flag_for_face_refinement(j);
5197 template <
int dim,
int spacedim>
5200 const double max_ratio,
5201 const unsigned int max_iterations)
5203 unsigned int iter = 0;
5204 bool continue_refinement =
true;
5206 while (continue_refinement && (iter < max_iterations))
5209 continue_refinement =
false;
5212 std::pair<unsigned int, double> info =
5213 GridTools::get_longest_direction<dim, spacedim>(cell);
5214 if (info.second > max_ratio)
5216 cell->set_refine_flag(
5218 continue_refinement =
true;
5226 template <
int dim,
int spacedim>
5229 const double limit_angle_fraction)
5237 "have hanging nodes."));
5241 bool has_cells_with_more_than_dim_faces_on_boundary =
true;
5242 bool has_cells_with_dim_faces_on_boundary =
false;
5244 unsigned int refinement_cycles = 0;
5246 while (has_cells_with_more_than_dim_faces_on_boundary)
5248 has_cells_with_more_than_dim_faces_on_boundary =
false;
5252 unsigned int boundary_face_counter = 0;
5253 for (
auto f : cell->face_indices())
5254 if (cell->face(f)->at_boundary())
5255 boundary_face_counter++;
5256 if (boundary_face_counter > dim)
5258 has_cells_with_more_than_dim_faces_on_boundary =
true;
5261 else if (boundary_face_counter == dim)
5262 has_cells_with_dim_faces_on_boundary =
true;
5264 if (has_cells_with_more_than_dim_faces_on_boundary)
5267 refinement_cycles++;
5271 if (has_cells_with_dim_faces_on_boundary)
5274 refinement_cycles++;
5278 while (refinement_cycles > 0)
5281 cell->set_coarsen_flag();
5283 refinement_cycles--;
5293 std::vector<CellData<dim>> cells_to_add;
5297 const unsigned int v0 = 0,
v1 = 1, v2 = (dim > 1 ? 2 : 0),
5298 v3 = (dim > 1 ? 3 : 0);
5302 double angle_fraction = 0;
5308 p0[spacedim > 1 ? 1 : 0] = 1;
5312 if (cell->face(
v0)->at_boundary() && cell->face(v3)->at_boundary())
5314 p0 = cell->vertex(
v0) - cell->vertex(v2);
5315 p1 = cell->vertex(v3) - cell->vertex(v2);
5316 vertex_at_corner = v2;
5318 else if (cell->face(v3)->at_boundary() &&
5319 cell->face(
v1)->at_boundary())
5321 p0 = cell->vertex(v2) - cell->vertex(v3);
5322 p1 = cell->vertex(
v1) - cell->vertex(v3);
5323 vertex_at_corner = v3;
5325 else if (cell->face(1)->at_boundary() &&
5326 cell->face(2)->at_boundary())
5328 p0 = cell->vertex(
v0) - cell->vertex(
v1);
5329 p1 = cell->vertex(v3) - cell->vertex(
v1);
5330 vertex_at_corner =
v1;
5332 else if (cell->face(2)->at_boundary() &&
5333 cell->face(0)->at_boundary())
5335 p0 = cell->vertex(v2) - cell->vertex(
v0);
5336 p1 = cell->vertex(
v1) - cell->vertex(
v0);
5337 vertex_at_corner =
v0;
5348 if (angle_fraction > limit_angle_fraction)
5350 auto flags_removal = [&](
unsigned int f1,
5353 unsigned int n2) ->
void {
5354 cells_to_remove[cell->active_cell_index()] =
true;
5355 cells_to_remove[cell->neighbor(n1)->active_cell_index()] =
true;
5356 cells_to_remove[cell->neighbor(n2)->active_cell_index()] =
true;
5358 faces_to_remove[cell->face(f1)->index()] =
true;
5359 faces_to_remove[cell->face(f2)->index()] =
true;
5361 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] =
true;
5362 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] =
true;
5365 auto cell_creation = [&](
const unsigned int vv0,
5366 const unsigned int vv1,
5367 const unsigned int f0,
5368 const unsigned int f1,
5370 const unsigned int n0,
5371 const unsigned int v0n0,
5372 const unsigned int v1n0,
5374 const unsigned int n1,
5375 const unsigned int v0n1,
5376 const unsigned int v1n1) {
5382 c1.
vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
5383 c1.
vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
5389 c2.
vertices[
v1] = cell->neighbor(n1)->vertex_index(v0n1);
5390 c2.
vertices[v2] = cell->vertex_index(vv1);
5391 c2.
vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
5396 l1.
vertices[0] = cell->vertex_index(vv0);
5397 l1.
vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
5403 l2.
vertices[0] = cell->vertex_index(vv0);
5404 l2.
vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
5410 cells_to_add.push_back(c1);
5411 cells_to_add.push_back(c2);
5416 switch (vertex_at_corner)
5419 flags_removal(0, 2, 3, 1);
5420 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
5423 flags_removal(1, 2, 3, 0);
5424 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
5427 flags_removal(3, 0, 1, 2);
5428 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
5431 flags_removal(3, 1, 0, 2);
5432 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
5445 if (cells_to_add.size() == 0)
5447 while (refinement_cycles > 0)
5450 cell->set_coarsen_flag();
5452 refinement_cycles--;
5460 if (cells_to_remove[cell->active_cell_index()] ==
false)
5463 for (
const unsigned int v : cell->vertex_indices())
5464 c.
vertices[v] = cell->vertex_index(v);
5467 cells_to_add.push_back(c);
5475 for (; face != endf; ++face)
5476 if ((face->at_boundary() ||
5478 faces_to_remove[face->index()] ==
false)
5480 for (
unsigned int l = 0;
l < face->
n_lines(); ++
l)
5485 for (
const unsigned int v : face->vertex_indices())
5486 line.
vertices[v] = face->vertex_index(v);
5492 for (
const unsigned int v : face->line(
l)->vertex_indices())
5493 line.
vertices[v] = face->line(
l)->vertex_index(v);
5502 for (
const unsigned int v : face->vertex_indices())
5503 quad.
vertices[v] = face->vertex_index(v);
5511 subcelldata_to_add);
5516 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
5535 template <
int dim,
int spacedim>
5538 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5539 std::vector<std::vector<Point<dim>>>,
5540 std::vector<std::vector<unsigned int>>>
5552 auto &cells = std::get<0>(cqmp);
5553 auto &qpoints = std::get<1>(cqmp);
5554 auto &maps = std::get<2>(cqmp);
5556 return std::make_tuple(std::move(cells),
5563 template <
int dim,
int spacedim>
5566 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5567 std::vector<std::vector<Point<dim>>>,
5568 std::vector<std::vector<unsigned int>>,
5569 std::vector<unsigned int>>
5579 Assert((dim == spacedim),
5580 ExcMessage(
"Only implemented for dim==spacedim."));
5589 const unsigned int np = points.size();
5591 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5593 std::vector<std::vector<Point<dim>>> qpoints_out;
5594 std::vector<std::vector<unsigned int>> maps_out;
5595 std::vector<unsigned int> missing_points_out;
5599 return std::make_tuple(std::move(cells_out),
5600 std::move(qpoints_out),
5601 std::move(maps_out),
5602 std::move(missing_points_out));
5610 std::vector<std::pair<Point<spacedim>,
unsigned int>> points_and_ids(np);
5611 for (
unsigned int i = 0; i < np; ++i)
5612 points_and_ids[i] = std::make_pair(points[i], i);
5613 const auto p_tree =
pack_rtree(points_and_ids);
5616 std::vector<bool> found_points(points.size(),
false);
5619 const auto already_found = [&found_points](
const auto &id) {
5621 return found_points[
id.second];
5627 const auto store_cell_point_and_id =
5631 const unsigned int &id) {
5632 const auto it = std::find(cells_out.rbegin(), cells_out.rend(), cell);
5633 if (it != cells_out.rend())
5635 const auto cell_id =
5636 (cells_out.size() - 1 - (it - cells_out.rbegin()));
5637 qpoints_out[cell_id].emplace_back(ref_point);
5638 maps_out[cell_id].emplace_back(
id);
5642 cells_out.emplace_back(cell);
5643 qpoints_out.emplace_back(std::vector<
Point<dim>>({ref_point}));
5644 maps_out.emplace_back(std::vector<unsigned int>({
id}));
5649 const auto check_all_points_within_box = [&](
const auto &leaf) {
5650 const auto &box = leaf.first;
5651 const auto &cell_hint = leaf.second;
5653 for (
const auto &point_and_id :
5654 p_tree | bgi::adaptors::queried(!bgi::satisfies(already_found) &&
5655 bgi::intersects(box)))
5657 const auto id = point_and_id.second;
5658 const auto cell_and_ref =
5662 const auto &cell = cell_and_ref.first;
5663 const auto &ref_point = cell_and_ref.second;
5666 store_cell_point_and_id(cell, ref_point,
id);
5668 missing_points_out.emplace_back(
id);
5671 found_points[id] =
true;
5677 check_all_points_within_box(
5678 std::make_pair(mapping.get_bounding_box(cell_hint), cell_hint));
5681 for (
unsigned int i = 0; i < np; ++i)
5682 if (found_points[i] ==
false)
5685 const auto leaf = b_tree.qbegin(bgi::nearest(points[i], 1));
5687 if (leaf != b_tree.qend())
5688 check_all_points_within_box(*leaf);
5696 for (
unsigned int i = 0; i < np; ++i)
5697 if (found_points[i] ==
false)
5698 missing_points_out.emplace_back(i);
5705 unsigned int c = cells_out.size();
5706 unsigned int qps = 0;
5712 for (
unsigned int n = 0; n < c; ++n)
5715 qps += qpoints_out[n].size();
5718 Assert(qps + missing_points_out.size() == np,
5722 return std::make_tuple(std::move(cells_out),
5723 std::move(qpoints_out),
5724 std::move(maps_out),
5725 std::move(missing_points_out));
5730 template <
int dim,
int spacedim>
5733 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5734 std::vector<std::vector<Point<dim>>>,
5735 std::vector<std::vector<unsigned int>>,
5736 std::vector<std::vector<Point<spacedim>>>,
5737 std::vector<std::vector<unsigned int>>>
5745 const double tolerance)
5749 cache, points, global_bboxes, {}, tolerance,
false,
true)
5754 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5755 std::vector<std::vector<Point<dim>>>,
5756 std::vector<std::vector<unsigned int>>,
5757 std::vector<std::vector<Point<spacedim>>>,
5758 std::vector<std::vector<unsigned int>>>
5761 std::pair<int, int> dummy{-1, -1};
5763 for (
unsigned int i = 0; i < all.size(); ++i)
5765 if (dummy != std::get<0>(all[i]))
5767 std::get<0>(result).push_back(
5770 std::get<0>(all[i]).
first,
5771 std::get<0>(all[i]).
second});
5773 const unsigned int new_size = std::get<0>(result).size();
5775 std::get<1>(result).resize(new_size);
5776 std::get<2>(result).resize(new_size);
5777 std::get<3>(result).resize(new_size);
5778 std::get<4>(result).resize(new_size);
5780 dummy = std::get<0>(all[i]);
5783 std::get<1>(result).back().push_back(
5784 std::get<3>(all[i]));
5785 std::get<2>(result).back().push_back(std::get<2>(all[i]));
5786 std::get<3>(result).back().push_back(std::get<4>(all[i]));
5787 std::get<4>(result).back().push_back(std::get<1>(all[i]));
5797 template <
int spacedim>
5798 std::tuple<std::vector<unsigned int>,
5799 std::vector<unsigned int>,
5800 std::vector<unsigned int>>
5804 const double tolerance)
5806 std::vector<std::pair<unsigned int, unsigned int>> ranks_and_indices;
5807 ranks_and_indices.reserve(points.size());
5809 for (
unsigned int i = 0; i < points.size(); ++i)
5811 const auto &
point = points[i];
5812 for (
unsigned rank = 0; rank < global_bboxes.size(); ++rank)
5813 for (
const auto &box : global_bboxes[rank])
5814 if (box.point_inside(
point, tolerance))
5816 ranks_and_indices.emplace_back(rank, i);
5822 std::sort(ranks_and_indices.begin(), ranks_and_indices.end());
5824 std::vector<unsigned int> ranks;
5825 std::vector<unsigned int> ptr;
5826 std::vector<unsigned int> indices;
5830 for (
const auto &i : ranks_and_indices)
5832 if (dummy_rank != i.first)
5834 dummy_rank = i.first;
5835 ranks.push_back(dummy_rank);
5836 ptr.push_back(indices.size());
5839 indices.push_back(i.second);
5841 ptr.push_back(indices.size());
5843 return std::make_tuple(std::move(ranks),
5845 std::move(indices));
5850 template <
int dim,
int spacedim>
5852 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5858 const std::vector<bool> &marked_vertices,
5859 const double tolerance,
5860 const bool enforce_unique_mapping)
5863 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5865 locally_owned_active_cells_around_point;
5882 cell_hint = first_cell.first;
5885 const auto active_cells_around_point =
5893 if (enforce_unique_mapping)
5900 for (
const auto &cell : active_cells_around_point)
5901 lowes_rank =
std::min(lowes_rank, cell.first->subdomain_id());
5903 if (lowes_rank != my_rank)
5907 locally_owned_active_cells_around_point.reserve(
5908 active_cells_around_point.size());
5910 for (
const auto &cell : active_cells_around_point)
5911 if (cell.first->is_locally_owned())
5912 locally_owned_active_cells_around_point.push_back(cell);
5915 std::sort(locally_owned_active_cells_around_point.begin(),
5916 locally_owned_active_cells_around_point.end(),
5917 [](
const auto &a,
const auto &
b) { return a.first < b.first; });
5919 if (enforce_unique_mapping &&
5920 locally_owned_active_cells_around_point.size() > 1)
5922 return {locally_owned_active_cells_around_point.front()};
5924 return locally_owned_active_cells_around_point;
5929 template <
int dim,
int spacedim>
5930 DistributedComputePointLocationsInternal<dim, spacedim>
5935 const std::vector<bool> & marked_vertices,
5936 const double tolerance,
5937 const bool perform_handshake,
5938 const bool enforce_unique_mapping)
5949 const auto potential_owners =
5952 const auto &potential_owners_ranks = std::get<0>(potential_owners);
5953 const auto &potential_owners_ptrs = std::get<1>(potential_owners);
5954 const auto &potential_owners_indices = std::get<2>(potential_owners);
5958 const auto translate = [&](
const unsigned int other_rank) {
5959 const auto ptr = std::find(potential_owners_ranks.begin(),
5960 potential_owners_ranks.end(),
5965 const auto other_rank_index =
5966 std::distance(potential_owners_ranks.begin(), ptr);
5968 return other_rank_index;
5972 (marked_vertices.size() == 0) ||
5975 "The marked_vertices vector has to be either empty or its size has "
5976 "to equal the number of vertices of the triangulation."));
5978 using RequestType = std::vector<std::pair<unsigned int, Point<spacedim>>>;
5979 using AnswerType = std::vector<unsigned int>;
5985 const bool has_relevant_vertices =
5986 (marked_vertices.size() == 0) ||
5987 (std::find(marked_vertices.begin(), marked_vertices.end(),
true) !=
5988 marked_vertices.end());
5990 const auto create_request = [&](
const unsigned int other_rank) {
5991 const auto other_rank_index = translate(other_rank);
5993 RequestType request;
5994 request.reserve(potential_owners_ptrs[other_rank_index + 1] -
5995 potential_owners_ptrs[other_rank_index]);
5997 for (
unsigned int i = potential_owners_ptrs[other_rank_index];
5998 i < potential_owners_ptrs[other_rank_index + 1];
6000 request.emplace_back(potential_owners_indices[i],
6001 points[potential_owners_indices[i]]);
6006 const auto answer_request =
6007 [&](
const unsigned int &other_rank,
6008 const RequestType & request) -> AnswerType {
6009 AnswerType answer(request.size(), 0);
6011 if (has_relevant_vertices)
6015 for (
unsigned int i = 0; i < request.size(); ++i)
6017 const auto &index_and_point = request[i];
6019 const auto cells_and_reference_positions =
6022 index_and_point.second,
6026 enforce_unique_mapping);
6028 for (
const auto &cell_and_reference_position :
6029 cells_and_reference_positions)
6031 send_components.emplace_back(
6032 std::pair<int, int>(
6033 cell_and_reference_position.first->level(),
6034 cell_and_reference_position.first->index()),
6036 index_and_point.first,
6037 cell_and_reference_position.second,
6038 index_and_point.second,
6042 answer[i] = cells_and_reference_positions.size();
6046 if (perform_handshake)
6052 const auto process_answer = [&](
const unsigned int other_rank,
6053 const AnswerType & answer) {
6054 if (perform_handshake)
6056 const auto other_rank_index = translate(other_rank);
6058 for (
unsigned int i = 0; i < answer.size(); ++i)
6059 for (
unsigned int j = 0; j < answer[i]; ++j)
6060 recv_components.emplace_back(
6062 potential_owners_indices
6063 [i + potential_owners_ptrs[other_rank_index]],
6068 Utilities::MPI::ConsensusAlgorithms::selector<RequestType, AnswerType>(
6069 potential_owners_ranks,
6079 std::sort(send_components.begin(),
6080 send_components.end(),
6081 [&](
const auto &a,
const auto &
b) {
6082 if (std::get<1>(a) != std::get<1>(b))
6083 return std::get<1>(a) < std::get<1>(b);
6085 if (std::get<2>(a) != std::get<2>(b))
6086 return std::get<2>(a) < std::get<2>(b);
6088 return std::get<0>(a) < std::get<0>(b);
6093 i < send_components.size();
6096 std::get<5>(send_components[i]) = i;
6098 if (dummy != std::get<1>(send_components[i]))
6100 dummy = std::get<1>(send_components[i]);
6101 send_ranks.push_back(dummy);
6102 send_ptrs.push_back(i);
6105 send_ptrs.push_back(send_components.size());
6109 std::sort(send_components.begin(),
6110 send_components.end(),
6111 [&](
const auto &a,
const auto &
b) {
6112 if (std::get<0>(a) != std::get<0>(b))
6113 return std::get<0>(a) < std::get<0>(b);
6115 if (std::get<1>(a) != std::get<1>(b))
6116 return std::get<1>(a) < std::get<1>(b);
6118 if (std::get<2>(a) != std::get<2>(b))
6119 return std::get<2>(a) < std::get<2>(b);
6121 return std::get<5>(a) < std::get<5>(b);
6125 if (perform_handshake)
6128 std::sort(recv_components.begin(),
6129 recv_components.end(),
6130 [&](
const auto &a,
const auto &
b) {
6131 if (std::get<0>(a) != std::get<0>(b))
6132 return std::get<0>(a) < std::get<0>(b);
6134 return std::get<1>(a) < std::get<1>(b);
6139 i < recv_components.size();
6142 std::get<2>(recv_components[i]) = i;
6144 if (dummy != std::get<0>(recv_components[i]))
6146 dummy = std::get<0>(recv_components[i]);
6147 recv_ranks.push_back(dummy);
6148 recv_ptrs.push_back(i);
6151 recv_ptrs.push_back(recv_components.size());
6155 std::sort(recv_components.begin(),
6156 recv_components.end(),
6157 [&](
const auto &a,
const auto &
b) {
6158 if (std::get<1>(a) != std::get<1>(b))
6159 return std::get<1>(a) < std::get<1>(b);
6161 if (std::get<0>(a) != std::get<0>(b))
6162 return std::get<0>(a) < std::get<0>(b);
6164 return std::get<2>(a) < std::get<2>(b);
6174 template <
int dim,
int spacedim>
6175 std::map<unsigned int, Point<spacedim>>
6179 std::map<unsigned int, Point<spacedim>> result;
6182 if (!cell->is_artificial())
6185 for (
unsigned int i = 0; i < vs.size(); ++i)
6186 result[cell->vertex_index(i)] = vs[i];
6193 template <
int spacedim>
6198 auto id_and_v = std::min_element(
6203 return p1.second.distance(p) < p2.second.distance(p);
6205 return id_and_v->first;
6209 template <
int dim,
int spacedim>
6210 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
6217 const std::vector<bool> &marked_vertices,
6218 const double tolerance)
6223 const auto &vertex_to_cell_centers =
6231 vertex_to_cell_centers,
6234 used_vertices_rtree,
6238 template <
int spacedim>
6239 std::vector<std::vector<BoundingBox<spacedim>>>
6242 const MPI_Comm & mpi_communicator)
6244 #ifndef DEAL_II_WITH_MPI
6246 (void)mpi_communicator;
6249 "GridTools::exchange_local_bounding_boxes() requires MPI."));
6253 unsigned int n_bboxes = local_bboxes.size();
6255 int n_local_data = 2 * spacedim * n_bboxes;
6257 std::vector<double> loc_data_array(n_local_data);
6258 for (
unsigned int i = 0; i < n_bboxes; ++i)
6259 for (
unsigned int d = 0;
d < spacedim; ++
d)
6262 loc_data_array[2 * i * spacedim +
d] =
6263 local_bboxes[i].get_boundary_points().first[
d];
6264 loc_data_array[2 * i * spacedim + spacedim +
d] =
6265 local_bboxes[i].get_boundary_points().second[
d];
6272 std::vector<int> size_all_data(n_procs);
6275 int ierr = MPI_Allgather(&n_local_data,
6278 size_all_data.data(),
6286 std::vector<int> rdispls(n_procs);
6288 for (
unsigned int i = 1; i < n_procs; ++i)
6289 rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
6293 std::vector<double> data_array(rdispls.back() + size_all_data.back());
6295 ierr = MPI_Allgatherv(loc_data_array.data(),
6299 size_all_data.data(),
6306 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
6307 unsigned int begin_idx = 0;
6308 for (
unsigned int i = 0; i < n_procs; ++i)
6311 unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
6312 global_bboxes[i].resize(n_bbox_i);
6313 for (
unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
6316 for (
unsigned int d = 0;
d < spacedim; ++
d)
6318 p1[
d] = data_array[begin_idx + 2 * bbox * spacedim +
d];
6320 data_array[begin_idx + 2 * bbox * spacedim + spacedim +
d];
6323 global_bboxes[i][bbox] = loc_bbox;
6326 begin_idx += size_all_data[i];
6328 return global_bboxes;
6334 template <
int spacedim>
6338 const MPI_Comm & mpi_communicator)
6340 #ifndef DEAL_II_WITH_MPI
6341 (void)mpi_communicator;
6343 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> boxes_index(
6344 local_description.size());
6346 for (
unsigned int i = 0; i < local_description.size(); ++i)
6347 boxes_index[i] = std::make_pair(local_description[i], 0u);
6351 const std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes =
6355 const unsigned int n_procs =