17 #include <deal.II/base/mpi.templates.h> 60 #include <boost/random/mersenne_twister.hpp> 61 #include <boost/random/uniform_real_distribution.hpp> 71 #include <unordered_map> 78 template <
int dim,
int spacedim>
88 #if defined(DEAL_II_WITH_P4EST) && defined(DEBUG) 103 std::vector<bool> boundary_vertices(vertices.size(),
false);
109 for (; cell != endc; ++cell)
110 for (
const unsigned int face : cell->face_indices())
111 if (cell->face(face)->at_boundary())
112 for (
unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
113 boundary_vertices[cell->face(face)->vertex_index(i)] =
true;
117 double max_distance_sqr = 0;
118 std::vector<bool>::const_iterator pi = boundary_vertices.begin();
119 const unsigned int N = boundary_vertices.size();
120 for (
unsigned int i = 0; i <
N; ++i, ++pi)
122 std::vector<bool>::const_iterator pj = pi + 1;
123 for (
unsigned int j = i + 1; j <
N; ++j, ++pj)
124 if ((*pi ==
true) && (*pj ==
true) &&
125 ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
126 max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
129 return std::sqrt(max_distance_sqr);
134 template <
int dim,
int spacedim>
140 unsigned int mapping_degree = 1;
143 mapping_degree = p->get_degree();
144 else if (
const auto *p =
146 mapping_degree = p->get_degree();
149 const QGauss<dim> quadrature_formula(mapping_degree + 1);
150 const unsigned int n_q_points = quadrature_formula.
size();
163 endc = triangulation.
end();
165 double local_volume = 0;
168 for (; cell != endc; ++cell)
169 if (cell->is_locally_owned())
171 fe_values.reinit(cell);
172 for (
unsigned int q = 0; q < n_q_points; ++q)
173 local_volume += fe_values.JxW(q);
176 double global_volume = 0;
178 #ifdef DEAL_II_WITH_MPI 186 global_volume = local_volume;
188 return global_volume;
210 struct TransformR2UAffine
225 [1] = {{-1.000000}, {1.000000}};
229 {1.000000, 0.000000};
240 [2] = {{-0.500000, -0.500000},
241 {0.500000, -0.500000},
242 {-0.500000, 0.500000},
243 {0.500000, 0.500000}};
253 {0.750000, 0.250000, 0.250000, -0.250000};
259 {-0.250000, -0.250000, -0.250000},
260 {0.250000, -0.250000, -0.250000},
261 {-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}
285 template <
int dim,
int spacedim>
294 for (
unsigned int d = 0;
d < spacedim; ++
d)
295 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
296 for (
unsigned int e = 0;
e < dim; ++
e)
301 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
304 return std::make_pair(A, b);
318 Vector<double> aspect_ratio_vector(triangulation.
n_active_cells());
323 if (cell->is_locally_owned())
325 double aspect_ratio_cell = 0.0;
330 for (
unsigned int q = 0; q < quadrature.
size(); ++q)
341 aspect_ratio_cell = std::numeric_limits<double>::infinity();
346 for (
unsigned int i = 0; i < dim; i++)
347 for (
unsigned int j = 0; j < dim; j++)
348 J(i, j) = jacobian[i][j];
354 double const ar = max_sv / min_sv;
361 aspect_ratio_cell =
std::max(aspect_ratio_cell, ar);
366 aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
370 return aspect_ratio_vector;
381 Vector<double> aspect_ratio_vector =
391 template <
int dim,
int spacedim>
397 const auto predicate = [](
const iterator &) {
return true; };
400 tria, std::function<
bool(
const iterator &)>(predicate));
426 template <
int structdim>
450 "Two CellData objects with equal vertices must " 451 "have the same material/boundary ids and manifold " 476 template <
class FaceIteratorType>
481 for (
unsigned int vertex_n = 0; vertex_n < face->n_vertices();
483 face_cell_data.
vertices[vertex_n] = face->vertex_index(vertex_n);
484 face_cell_data.boundary_id = face->boundary_id();
485 face_cell_data.manifold_id = face->manifold_id();
487 face_data.insert(face_cell_data);
515 template <
class FaceIteratorType>
530 template <
int dim,
int spacedim>
532 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>,
SubCellData>
536 ExcMessage(
"The input triangulation must be non-empty."));
538 std::vector<Point<spacedim>>
vertices;
539 std::vector<CellData<dim>> cells;
541 unsigned int max_level_0_vertex_n = 0;
543 for (
const unsigned int cell_vertex_n : cell->vertex_indices())
544 max_level_0_vertex_n =
545 std::max(cell->vertex_index(cell_vertex_n), max_level_0_vertex_n);
546 vertices.resize(max_level_0_vertex_n + 1);
556 for (
const unsigned int cell_vertex_n : cell->vertex_indices())
558 Assert(cell->vertex_index(cell_vertex_n) < vertices.size(),
560 vertices[cell->vertex_index(cell_vertex_n)] =
561 cell->vertex(cell_vertex_n);
562 cell_data.vertices[cell_vertex_n] =
563 cell->vertex_index(cell_vertex_n);
565 cell_data.material_id = cell->material_id();
566 cell_data.manifold_id = cell->manifold_id();
567 cells.push_back(cell_data);
572 for (
const unsigned int face_n : cell->face_indices())
578 for (
unsigned int line_n = 0; line_n < cell->n_lines(); ++line_n)
580 const auto line = cell->line(line_n);
582 for (
unsigned int vertex_n = 0; vertex_n < line->n_vertices();
585 line->vertex_index(vertex_n);
589 line_data.insert(line_cell_data);
597 std::vector<bool> used_vertices(vertices.size());
599 for (
const auto v : cell_data.vertices)
600 used_vertices[v] =
true;
601 Assert(std::find(used_vertices.begin(), used_vertices.end(),
false) ==
603 ExcMessage(
"The level zero vertices should form a contiguous " 611 for (
const CellData<1> &face_line_data : line_data)
614 return std::tuple<std::vector<Point<spacedim>>,
615 std::vector<CellData<dim>>,
618 std::move(subcell_data));
623 template <
int dim,
int spacedim>
632 "Invalid SubCellData supplied according to ::check_consistency(). " 633 "This is caused by data containing objects for the wrong dimension."));
636 std::vector<bool> vertex_used(
vertices.size(),
false);
637 for (
unsigned int c = 0; c < cells.size(); ++c)
638 for (
unsigned int v = 0; v < cells[c].vertices.size(); ++v)
641 ExcMessage(
"Invalid vertex index encountered! cells[" +
645 " is invalid, because only " +
647 " vertices were supplied."));
648 vertex_used[cells[c].vertices[v]] =
true;
655 std::vector<unsigned int> new_vertex_numbers(
vertices.size(),
657 unsigned int next_free_number = 0;
658 for (
unsigned int i = 0; i <
vertices.size(); ++i)
659 if (vertex_used[i] ==
true)
661 new_vertex_numbers[i] = next_free_number;
666 for (
unsigned int c = 0; c < cells.size(); ++c)
668 v = new_vertex_numbers[v];
673 for (
unsigned int v = 0;
678 new_vertex_numbers.size(),
680 "Invalid vertex index in subcelldata.boundary_lines. " 681 "subcelldata.boundary_lines[" +
686 " is invalid, because only " +
688 " vertices were supplied."));
695 for (
unsigned int v = 0;
700 new_vertex_numbers.size(),
702 "Invalid vertex index in subcelldata.boundary_quads. " 703 "subcelldata.boundary_quads[" +
708 " is invalid, because only " +
710 " vertices were supplied."));
718 std::vector<Point<spacedim>> tmp;
719 tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(),
true));
720 for (
unsigned int v = 0; v < vertices.size(); ++v)
721 if (vertex_used[v] ==
true)
722 tmp.push_back(vertices[v]);
728 template <
int dim,
int spacedim>
733 std::vector<unsigned int> & considered_vertices,
739 std::vector<unsigned int> new_vertex_numbers(
vertices.size());
740 std::iota(new_vertex_numbers.begin(), new_vertex_numbers.end(), 0);
743 if (considered_vertices.size() == 0)
744 considered_vertices = new_vertex_numbers;
758 const auto &
min = bbox.get_boundary_points().first;
759 const auto &
max = bbox.get_boundary_points().second;
761 unsigned int longest_coordinate_direction = 0;
762 double longest_coordinate_length =
max[0] -
min[0];
763 for (
unsigned int d = 1;
d < spacedim; ++
d)
765 const double coordinate_length =
max[
d] - min[
d];
766 if (longest_coordinate_length < coordinate_length)
768 longest_coordinate_length = coordinate_length;
769 longest_coordinate_direction =
d;
775 std::vector<std::pair<unsigned int, Point<spacedim>>> sorted_vertices;
776 sorted_vertices.reserve(
vertices.size());
777 for (
const unsigned int vertex_n : considered_vertices)
780 sorted_vertices.emplace_back(vertex_n,
vertices[vertex_n]);
782 std::sort(sorted_vertices.begin(),
783 sorted_vertices.end(),
784 [&](
const std::pair<unsigned int, Point<spacedim>> &a,
785 const std::pair<unsigned int, Point<spacedim>> &
b) {
786 return a.second[longest_coordinate_direction] <
787 b.second[longest_coordinate_direction];
792 for (
unsigned int d = 0;
d < spacedim; ++
d)
793 if (std::abs(a[
d] -
b[
d]) > tol)
800 auto range_start = sorted_vertices.begin();
801 while (range_start != sorted_vertices.end())
803 auto range_end = range_start + 1;
804 while (range_end != sorted_vertices.end() &&
805 std::abs(range_end->second[longest_coordinate_direction] -
806 range_start->second[longest_coordinate_direction]) <
812 std::sort(range_start,
816 return a.first <
b.first;
825 for (
auto reference = range_start; reference != range_end; ++reference)
828 for (
auto it = reference + 1; it != range_end; ++it)
830 if (within_tolerance(reference->second, it->second))
832 new_vertex_numbers[it->first] = reference->first;
838 range_start = range_end;
845 for (
auto &cell : cells)
846 for (
auto &vertex_index : cell.vertices)
847 vertex_index = new_vertex_numbers[vertex_index];
849 for (
auto &vertex_index : quad.vertices)
850 vertex_index = new_vertex_numbers[vertex_index];
852 for (
auto &vertex_index : line.vertices)
853 vertex_index = new_vertex_numbers[vertex_index];
863 template <
int spacedim>
895 std::cos(
angle) * p(1) - std::sin(
angle) * p(2),
896 std::sin(
angle) * p(1) + std::cos(
angle) * p(2)};
898 return {std::cos(
angle) * p(0) + std::sin(
angle) * p(2),
900 -std::sin(
angle) * p(0) + std::cos(
angle) * p(2)};
902 return {std::cos(
angle) * p(0) - std::sin(
angle) * p(1),
903 std::sin(
angle) * p(0) + std::cos(
angle) * p(1),
912 template <
int spacedim>
931 template <
int dim,
int spacedim>
943 const unsigned int axis,
951 template <
int dim,
int spacedim>
973 const unsigned int n_dofs = S.
n();
983 Vector<double> f(n_dofs);
985 const auto constrained_rhs =
987 solver.
solve(SF, u, constrained_rhs, prec);
1000 const bool solve_for_absolute_positions)
1028 std::array<AffineConstraints<double>, dim> constraints;
1029 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
1037 for (
const unsigned int vertex_no : cell->vertex_indices())
1039 const unsigned int vertex_index = cell->vertex_index(vertex_no);
1040 const Point<dim> & vertex_point = cell->vertex(vertex_no);
1042 const typename std::map<unsigned int, Point<dim>>::const_iterator
1043 map_iter = new_points.find(vertex_index);
1045 if (map_iter != map_end)
1046 for (
unsigned int i = 0; i < dim; ++i)
1048 constraints[i].add_line(cell->vertex_dof_index(vertex_no, 0));
1049 constraints[i].set_inhomogeneity(
1050 cell->vertex_dof_index(vertex_no, 0),
1051 (solve_for_absolute_positions ?
1052 map_iter->second(i) :
1053 map_iter->second(i) - vertex_point[i]));
1058 for (
unsigned int i = 0; i < dim; ++i)
1059 constraints[i].close();
1062 Vector<double> us[dim];
1063 for (
unsigned int i = 0; i < dim; ++i)
1068 for (
unsigned int i = 0; i < dim; ++i)
1075 std::vector<bool> vertex_touched(triangulation.
n_vertices(),
false);
1077 for (
const unsigned int vertex_no : cell->vertex_indices())
1078 if (vertex_touched[cell->vertex_index(vertex_no)] ==
false)
1083 cell->vertex_dof_index(vertex_no, 0);
1084 for (
unsigned int i = 0; i < dim; ++i)
1085 if (solve_for_absolute_positions)
1086 v(i) = us[i](dof_index);
1088 v(i) += us[i](dof_index);
1090 vertex_touched[cell->vertex_index(vertex_no)] =
true;
1094 template <
int dim,
int spacedim>
1095 std::map<unsigned int, Point<spacedim>>
1098 std::map<unsigned int, Point<spacedim>> vertex_map;
1102 for (; cell != endc; ++cell)
1104 for (
unsigned int i : cell->face_indices())
1108 if (face->at_boundary())
1110 for (
unsigned j = 0; j < face->n_vertices(); ++j)
1113 const unsigned int vertex_index = face->vertex_index(j);
1114 vertex_map[vertex_index] = vertex;
1126 template <
int dim,
int spacedim>
1130 const bool keep_boundary,
1131 const unsigned int seed)
1148 double almost_infinite_length = 0;
1150 triangulation.
begin(0);
1151 cell != triangulation.
end(0);
1153 almost_infinite_length += cell->diameter();
1155 std::vector<double> minimal_length(triangulation.
n_vertices(),
1156 almost_infinite_length);
1159 std::vector<bool> at_boundary(keep_boundary ? triangulation.
n_vertices() :
1164 const bool is_parallel_shared =
1168 if (is_parallel_shared || cell->is_locally_owned())
1172 for (
unsigned int i = 0; i < cell->n_lines(); ++i)
1175 line = cell->line(i);
1177 if (keep_boundary && line->at_boundary())
1179 at_boundary[line->vertex_index(0)] =
true;
1180 at_boundary[line->vertex_index(1)] =
true;
1183 minimal_length[line->vertex_index(0)] =
1185 minimal_length[line->vertex_index(0)]);
1186 minimal_length[line->vertex_index(1)] =
1188 minimal_length[line->vertex_index(1)]);
1194 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
1195 if (cell->at_boundary(vertex) ==
true)
1196 at_boundary[cell->vertex_index(vertex)] =
true;
1198 minimal_length[cell->vertex_index(0)] =
1200 minimal_length[cell->vertex_index(0)]);
1201 minimal_length[cell->vertex_index(1)] =
1203 minimal_length[cell->vertex_index(1)]);
1208 boost::random::mt19937 rng(seed);
1209 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
1213 if (
auto distributed_triangulation =
1217 const std::vector<bool> locally_owned_vertices =
1219 std::vector<bool> vertex_moved(triangulation.
n_vertices(),
false);
1223 if (cell->is_locally_owned())
1225 for (
const unsigned int vertex_no : cell->vertex_indices())
1227 const unsigned global_vertex_no =
1228 cell->vertex_index(vertex_no);
1233 if ((keep_boundary && at_boundary[global_vertex_no]) ||
1234 vertex_moved[global_vertex_no] ||
1235 !locally_owned_vertices[global_vertex_no])
1240 for (
unsigned int d = 0;
d < spacedim; ++
d)
1241 shift_vector(
d) = uniform_distribution(rng);
1243 shift_vector *= factor * minimal_length[global_vertex_no] /
1244 std::sqrt(shift_vector.
square());
1247 cell->vertex(vertex_no) += shift_vector;
1248 vertex_moved[global_vertex_no] =
true;
1252 distributed_triangulation->communicate_locally_moved_vertices(
1253 locally_owned_vertices);
1262 const unsigned int n_vertices = triangulation.
n_vertices();
1263 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
1264 const std::vector<Point<spacedim>> &old_vertex_locations =
1267 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1271 if (keep_boundary && at_boundary[vertex])
1272 new_vertex_locations[vertex] = old_vertex_locations[vertex];
1277 for (
unsigned int d = 0;
d < spacedim; ++
d)
1278 shift_vector(
d) = uniform_distribution(rng);
1280 shift_vector *= factor * minimal_length[vertex] /
1281 std::sqrt(shift_vector.
square());
1284 new_vertex_locations[vertex] =
1285 old_vertex_locations[vertex] + shift_vector;
1291 for (
const unsigned int vertex_no : cell->vertex_indices())
1292 cell->vertex(vertex_no) =
1293 new_vertex_locations[cell->vertex_index(vertex_no)];
1308 endc = triangulation.
end();
1309 for (; cell != endc; ++cell)
1310 if (!cell->is_artificial())
1311 for (
const unsigned int face : cell->face_indices())
1312 if (cell->face(face)->has_children() &&
1313 !cell->face(face)->at_boundary())
1317 cell->face(face)->child(0)->vertex(1) =
1318 (cell->face(face)->vertex(0) +
1319 cell->face(face)->vertex(1)) /
1323 cell->face(face)->child(0)->vertex(1) =
1324 .5 * (cell->face(face)->vertex(0) +
1325 cell->face(face)->vertex(1));
1326 cell->face(face)->child(0)->vertex(2) =
1327 .5 * (cell->face(face)->vertex(0) +
1328 cell->face(face)->vertex(2));
1329 cell->face(face)->child(1)->vertex(3) =
1330 .5 * (cell->face(face)->vertex(1) +
1331 cell->face(face)->vertex(3));
1332 cell->face(face)->child(2)->vertex(3) =
1333 .5 * (cell->face(face)->vertex(2) +
1334 cell->face(face)->vertex(3));
1337 cell->face(face)->child(0)->vertex(3) =
1338 .25 * (cell->face(face)->vertex(0) +
1339 cell->face(face)->vertex(1) +
1340 cell->face(face)->vertex(2) +
1341 cell->face(face)->vertex(3));
1349 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1353 const std::vector<bool> & marked_vertices)
1362 marked_vertices.size() == 0,
1364 marked_vertices.size()));
1371 marked_vertices.size() == 0 ||
1372 std::equal(marked_vertices.begin(),
1373 marked_vertices.end(),
1375 [](
bool p,
bool q) {
return !p || q; }),
1377 "marked_vertices should be a subset of used vertices in the triangulation " 1378 "but marked_vertices contains one or more vertices that are not used vertices!"));
1382 const std::vector<bool> &vertices_to_use = (marked_vertices.size() == 0) ?
1388 std::vector<bool>::const_iterator
first =
1389 std::find(vertices_to_use.begin(), vertices_to_use.end(),
true);
1394 unsigned int best_vertex = std::distance(vertices_to_use.begin(),
first);
1395 double best_dist = (p - vertices[best_vertex]).norm_square();
1399 for (
unsigned int j = best_vertex + 1; j < vertices.size(); j++)
1400 if (vertices_to_use[j])
1402 const double dist = (p - vertices[j]).norm_square();
1403 if (dist < best_dist)
1415 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1418 const MeshType<dim, spacedim> &mesh,
1420 const std::vector<bool> & marked_vertices)
1433 marked_vertices.size() == 0,
1435 marked_vertices.size()));
1443 marked_vertices.size() == 0 ||
1444 std::equal(marked_vertices.begin(),
1445 marked_vertices.end(),
1447 [](
bool p,
bool q) {
return !p || q; }),
1449 "marked_vertices should be a subset of used vertices in the triangulation " 1450 "but marked_vertices contains one or more vertices that are not used vertices!"));
1453 if (marked_vertices.size() != 0)
1456 if (marked_vertices[it->first] ==
false)
1471 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1473 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1476 typename ::internal::
1477 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1480 const unsigned int vertex)
1486 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
1493 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1529 for (
const auto &cell : mesh.active_cell_iterators())
1531 for (
const unsigned int v : cell->vertex_indices())
1532 if (cell->vertex_index(v) == vertex)
1537 adjacent_cells.insert(cell);
1544 for (
const auto face :
1545 cell->reference_cell().faces_for_given_vertex(v))
1546 if (!cell->at_boundary(face) &&
1547 cell->neighbor(face)->is_active())
1559 adjacent_cells.insert(cell->neighbor(face));
1569 for (
unsigned int e = 0;
e < cell->n_lines(); ++
e)
1570 if (cell->line(
e)->has_children())
1574 if (cell->line(
e)->child(0)->vertex_index(1) == vertex)
1576 adjacent_cells.insert(cell);
1598 typename ::internal::
1599 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
1600 adjacent_cells.begin(), adjacent_cells.end());
1605 template <
int dim,
int spacedim>
1606 std::vector<std::vector<Tensor<1, spacedim>>>
1614 const unsigned int n_vertices = vertex_to_cells.size();
1619 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
1621 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1624 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
1625 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
1627 typename std::set<typename Triangulation<dim, spacedim>::
1628 active_cell_iterator>::iterator it =
1629 vertex_to_cells[vertex].begin();
1630 for (
unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
1632 vertex_to_cell_centers[vertex][cell] =
1633 (*it)->center() - vertices[vertex];
1634 vertex_to_cell_centers[vertex][cell] /=
1635 vertex_to_cell_centers[vertex][cell].norm();
1638 return vertex_to_cell_centers;
1644 template <
int spacedim>
1647 const unsigned int a,
1648 const unsigned int b,
1652 const double scalar_product_a = center_directions[a] * point_direction;
1653 const double scalar_product_b = center_directions[
b] * point_direction;
1658 return (scalar_product_a > scalar_product_b);
1662 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1664 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1666 std::pair<typename ::internal::
1667 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1672 const MeshType<dim, spacedim> &mesh,
1675 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
1678 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
1679 const std::vector<bool> & marked_vertices,
1681 const double tolerance)
1683 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1688 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1690 cell_and_position_approx;
1692 bool found_cell =
false;
1693 bool approx_cell =
false;
1695 unsigned int closest_vertex_index = 0;
1697 auto current_cell = cell_hint;
1699 while (found_cell ==
false)
1705 const auto cell_vertices = mapping.
get_vertices(current_cell);
1706 const unsigned int closest_vertex =
1707 find_closest_vertex_of_cell<dim, spacedim>(current_cell,
1710 vertex_to_point = p - cell_vertices[closest_vertex];
1711 closest_vertex_index = current_cell->vertex_index(closest_vertex);
1715 if (!used_vertices_rtree.empty())
1718 using ValueType = std::pair<Point<spacedim>,
unsigned int>;
1719 std::function<bool(const ValueType &)> marked;
1720 if (marked_vertices.size() == mesh.n_vertices())
1721 marked = [&marked_vertices](
const ValueType &value) ->
bool {
1722 return marked_vertices[value.second];
1725 marked = [](
const ValueType &) ->
bool {
return true; };
1727 std::vector<std::pair<Point<spacedim>,
unsigned int>> res;
1728 used_vertices_rtree.query(
1729 boost::geometry::index::nearest(p, 1) &&
1730 boost::geometry::index::satisfies(marked),
1731 std::back_inserter(res));
1735 closest_vertex_index = res[0].second;
1740 mapping, mesh, p, marked_vertices);
1742 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
1745 const double vertex_point_norm = vertex_to_point.
norm();
1746 if (vertex_point_norm > 0)
1747 vertex_to_point /= vertex_point_norm;
1749 const unsigned int n_neighbor_cells =
1750 vertex_to_cells[closest_vertex_index].size();
1753 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
1755 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
1756 neighbor_permutation[i] = i;
1758 auto comp = [&](
const unsigned int a,
const unsigned int b) ->
bool {
1759 return internal::compare_point_association<spacedim>(
1763 vertex_to_cell_centers[closest_vertex_index]);
1766 std::sort(neighbor_permutation.begin(),
1767 neighbor_permutation.end(),
1772 double best_distance = tolerance;
1776 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
1780 auto cell = vertex_to_cells[closest_vertex_index].begin();
1783 if (!(*cell)->is_artificial())
1790 cell_and_position.first = *cell;
1791 cell_and_position.second = p_unit;
1793 approx_cell =
false;
1802 if (dist < best_distance)
1804 best_distance = dist;
1805 cell_and_position_approx.first = *cell;
1806 cell_and_position_approx.second = p_unit;
1815 if (found_cell ==
true)
1816 return cell_and_position;
1817 else if (approx_cell ==
true)
1818 return cell_and_position_approx;
1834 mapping, mesh, p, marked_vertices, tolerance);
1836 current_cell =
typename MeshType<dim, spacedim>::active_cell_iterator();
1838 return cell_and_position;
1843 template <
int dim,
int spacedim>
1852 unsigned int closest_vertex = 0;
1854 for (
unsigned int v = 1; v < cell->n_vertices(); ++v)
1857 if (vertex_distance < minimum_distance)
1860 minimum_distance = vertex_distance;
1863 return closest_vertex;
1870 namespace BoundingBoxPredicate
1872 template <
class MeshType>
1873 std::tuple<BoundingBox<MeshType::space_dimension>,
bool>
1875 const typename MeshType::cell_iterator &parent_cell,
1876 const std::function<
1877 bool(
const typename MeshType::active_cell_iterator &)> &predicate)
1879 bool has_predicate =
1881 std::vector<typename MeshType::active_cell_iterator> active_cells;
1882 if (parent_cell->is_active())
1883 active_cells = {parent_cell};
1887 active_cells = get_active_child_cells<MeshType>(parent_cell);
1889 const unsigned int spacedim = MeshType::space_dimension;
1893 while (i < active_cells.size() && !predicate(active_cells[i]))
1897 if (active_cells.size() == 0 || i == active_cells.size())
1900 return std::make_tuple(bbox, has_predicate);
1907 for (; i < active_cells.size(); ++i)
1908 if (predicate(active_cells[i]))
1909 for (
const unsigned int v : active_cells[i]->vertex_indices())
1910 for (
unsigned int d = 0;
d < spacedim; ++
d)
1912 minp[
d] =
std::min(minp[
d], active_cells[i]->vertex(v)[d]);
1913 maxp[
d] =
std::max(maxp[d], active_cells[i]->vertex(v)[d]);
1916 has_predicate =
true;
1918 return std::make_tuple(bbox, has_predicate);
1925 template <
class MeshType>
1926 std::vector<BoundingBox<MeshType::space_dimension>>
1928 const MeshType &mesh,
1929 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1931 const unsigned int refinement_level,
1932 const bool allow_merge,
1933 const unsigned int max_boxes)
1940 refinement_level <= mesh.n_levels(),
1942 "Error: refinement level is higher then total levels in the triangulation!"));
1944 const unsigned int spacedim = MeshType::space_dimension;
1945 std::vector<BoundingBox<spacedim>> bounding_boxes;
1949 for (
unsigned int i = 0; i < refinement_level; ++i)
1950 for (
const typename MeshType::cell_iterator &cell :
1951 mesh.active_cell_iterators_on_level(i))
1953 bool has_predicate =
false;
1955 std::tie(bbox, has_predicate) =
1957 MeshType>(cell, predicate);
1959 bounding_boxes.push_back(bbox);
1963 for (
const typename MeshType::cell_iterator &cell :
1964 mesh.cell_iterators_on_level(refinement_level))
1966 bool has_predicate =
false;
1968 std::tie(bbox, has_predicate) =
1970 MeshType>(cell, predicate);
1972 bounding_boxes.push_back(bbox);
1977 return bounding_boxes;
1983 std::vector<unsigned int> merged_boxes_idx;
1984 bool found_neighbors =
true;
1989 while (found_neighbors)
1991 found_neighbors =
false;
1992 for (
unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
1994 if (std::find(merged_boxes_idx.begin(),
1995 merged_boxes_idx.end(),
1996 i) == merged_boxes_idx.end())
1997 for (
unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
1998 if (std::find(merged_boxes_idx.begin(),
1999 merged_boxes_idx.end(),
2000 j) == merged_boxes_idx.end() &&
2001 bounding_boxes[i].get_neighbor_type(
2002 bounding_boxes[j]) ==
2005 bounding_boxes[i].merge_with(bounding_boxes[j]);
2006 merged_boxes_idx.push_back(j);
2007 found_neighbors =
true;
2013 std::vector<BoundingBox<spacedim>> merged_b_boxes;
2014 for (
unsigned int i = 0; i < bounding_boxes.size(); ++i)
2015 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
2016 merged_boxes_idx.end())
2017 merged_b_boxes.push_back(bounding_boxes[i]);
2022 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
2024 std::vector<double> volumes;
2025 for (
unsigned int i = 0; i < merged_b_boxes.size(); ++i)
2026 volumes.push_back(merged_b_boxes[i].volume());
2028 while (merged_b_boxes.size() > max_boxes)
2030 unsigned int min_idx =
2031 std::min_element(volumes.begin(), volumes.end()) -
2033 volumes.erase(volumes.begin() + min_idx);
2035 bool not_removed =
true;
2036 for (
unsigned int i = 0;
2037 i < merged_b_boxes.size() && not_removed;
2042 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
2043 merged_b_boxes[min_idx]) ==
2045 merged_b_boxes[i].get_neighbor_type(
2046 merged_b_boxes[min_idx]) ==
2049 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
2050 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
2051 not_removed =
false;
2054 ExcMessage(
"Error: couldn't merge bounding boxes!"));
2057 Assert(merged_b_boxes.size() <= max_boxes,
2059 "Error: couldn't reach target number of bounding boxes!"));
2060 return merged_b_boxes;
2066 template <
int spacedim>
2068 std::tuple<std::vector<std::vector<unsigned int>>,
2069 std::map<unsigned int, unsigned int>,
2070 std::map<unsigned int, std::vector<unsigned int>>>
2078 unsigned int n_procs = global_bboxes.size();
2079 std::vector<std::vector<unsigned int>> point_owners(n_procs);
2080 std::map<unsigned int, unsigned int> map_owners_found;
2081 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
2083 unsigned int n_points = points.size();
2084 for (
unsigned int pt = 0; pt < n_points; ++pt)
2087 std::vector<unsigned int> owners_found;
2089 for (
unsigned int rk = 0; rk < n_procs; ++rk)
2092 if (bbox.point_inside(points[pt]))
2094 point_owners[rk].emplace_back(pt);
2095 owners_found.emplace_back(rk);
2099 Assert(owners_found.size() > 0,
2100 ExcMessage(
"No owners found for the point " +
2102 if (owners_found.size() == 1)
2103 map_owners_found[pt] = owners_found[0];
2106 map_owners_guessed[pt] = owners_found;
2109 return std::make_tuple(std::move(point_owners),
2110 std::move(map_owners_found),
2111 std::move(map_owners_guessed));
2114 template <
int spacedim>
2116 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2117 std::map<unsigned int, unsigned int>,
2118 std::map<unsigned int, std::vector<unsigned int>>>
2126 std::map<unsigned int, std::vector<unsigned int>> point_owners;
2127 std::map<unsigned int, unsigned int> map_owners_found;
2128 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
2129 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> search_result;
2131 unsigned int n_points = points.size();
2132 for (
unsigned int pt_n = 0; pt_n < n_points; ++pt_n)
2134 search_result.clear();
2137 covering_rtree.query(boost::geometry::index::intersects(points[pt_n]),
2138 std::back_inserter(search_result));
2141 std::set<unsigned int> owners_found;
2143 for (
const auto &rank_bbox : search_result)
2147 const bool pt_inserted = owners_found.insert(pt_n).second;
2149 point_owners[rank_bbox.second].emplace_back(pt_n);
2151 Assert(owners_found.size() > 0,
2152 ExcMessage(
"No owners found for the point " +
2154 if (owners_found.size() == 1)
2155 map_owners_found[pt_n] = *owners_found.begin();
2160 std::back_inserter(map_owners_guessed[pt_n]));
2163 return std::make_tuple(std::move(point_owners),
2164 std::move(map_owners_found),
2165 std::move(map_owners_guessed));
2169 template <
int dim,
int spacedim>
2171 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2175 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2179 endc = triangulation.
end();
2180 for (; cell != endc; ++cell)
2181 for (
const unsigned int i : cell->vertex_indices())
2186 for (; cell != endc; ++cell)
2188 for (
unsigned int i : cell->face_indices())
2190 if ((cell->at_boundary(i) ==
false) &&
2191 (cell->neighbor(i)->is_active()))
2194 adjacent_cell = cell->neighbor(i);
2195 for (
unsigned int j = 0; j < cell->face(i)->n_vertices(); ++j)
2204 for (
unsigned int i = 0; i < cell->n_lines(); ++i)
2205 if (cell->line(i)->has_children())
2219 template <
int dim,
int spacedim>
2220 std::map<unsigned int, types::global_vertex_index>
2224 std::map<unsigned int, types::global_vertex_index>
2225 local_to_global_vertex_index;
2227 #ifndef DEAL_II_WITH_MPI 2232 (void)triangulation;
2234 ExcMessage(
"This function does not make any sense " 2235 "for parallel::distributed::Triangulation " 2236 "objects if you do not have MPI enabled."));
2240 using active_cell_iterator =
2242 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
2247 unsigned int max_cellid_size = 0;
2248 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
2250 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
2256 active_cell_iterator cell = triangulation.
begin_active(),
2257 endc = triangulation.
end();
2258 std::set<active_cell_iterator> missing_vert_cells;
2259 std::set<unsigned int> used_vertex_index;
2260 for (; cell != endc; ++cell)
2262 if (cell->is_locally_owned())
2264 for (
const unsigned int i : cell->vertex_indices())
2266 types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
2267 typename std::set<active_cell_iterator>::iterator
2268 adjacent_cell = vertex_to_cell[cell->vertex_index(i)].begin(),
2269 end_adj_cell = vertex_to_cell[cell->vertex_index(i)].end();
2270 for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
2271 lowest_subdomain_id =
2273 (*adjacent_cell)->subdomain_id());
2276 if (lowest_subdomain_id == cell->subdomain_id())
2280 if (used_vertex_index.find(cell->vertex_index(i)) ==
2281 used_vertex_index.end())
2284 local_to_global_vertex_index[cell->vertex_index(i)] =
2290 vertex_to_cell[cell->vertex_index(i)].begin();
2291 for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
2292 if ((*adjacent_cell)->subdomain_id() !=
2293 cell->subdomain_id())
2296 types::global_vertex_index>
2297 tmp((*adjacent_cell)->subdomain_id(),
2298 cell->vertex_index(i));
2299 if (vertices_added.find(tmp) ==
2300 vertices_added.end())
2302 vertices_to_send[(*adjacent_cell)
2305 cell->vertex_index(i),
2306 cell->id().to_string());
2307 if (cell->id().to_string().size() >
2310 cell->id().to_string().size();
2311 vertices_added.insert(tmp);
2314 used_vertex_index.insert(cell->vertex_index(i));
2321 vertices_to_recv[lowest_subdomain_id].insert(
2322 cell->vertex_index(i));
2323 missing_vert_cells.insert(cell);
2330 if (cell->is_ghost())
2332 for (
unsigned int i : cell->face_indices())
2334 if (cell->at_boundary(i) ==
false)
2336 if (cell->neighbor(i)->is_active())
2339 spacedim>::active_cell_iterator
2340 adjacent_cell = cell->neighbor(i);
2341 if ((adjacent_cell->is_locally_owned()))
2343 types::subdomain_id adj_subdomain_id =
2344 adjacent_cell->subdomain_id();
2345 if (cell->subdomain_id() < adj_subdomain_id)
2346 for (
unsigned int j = 0;
2347 j < cell->face(i)->n_vertices();
2350 vertices_to_recv[cell->subdomain_id()].insert(
2351 cell->face(i)->vertex_index(j));
2352 missing_vert_cells.insert(cell);
2367 types::global_vertex_index
shift = 0;
2368 int ierr = MPI_Exscan(&next_index,
2376 std::map<unsigned int, types::global_vertex_index>::iterator
2377 global_index_it = local_to_global_vertex_index.begin(),
2378 global_index_end = local_to_global_vertex_index.end();
2379 for (; global_index_it != global_index_end; ++global_index_it)
2380 global_index_it->second += shift;
2394 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
2395 vertices_to_send.size());
2396 std::vector<MPI_Request> first_requests(vertices_to_send.size());
2400 std::string>>>::iterator
2401 vert_to_send_it = vertices_to_send.begin(),
2402 vert_to_send_end = vertices_to_send.end();
2403 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
2404 ++vert_to_send_it, ++i)
2406 int destination = vert_to_send_it->first;
2407 const unsigned int n_vertices = vert_to_send_it->second.size();
2408 const int buffer_size = 2 * n_vertices;
2409 vertices_send_buffers[i].resize(buffer_size);
2412 for (
unsigned int j = 0; j < n_vertices; ++j)
2414 vertices_send_buffers[i][2 * j] =
2415 std::get<0>(vert_to_send_it->second[j]);
2416 vertices_send_buffers[i][2 * j + 1] =
2417 local_to_global_vertex_index[std::get<1>(
2418 vert_to_send_it->second[j])];
2422 ierr = MPI_Isend(vertices_send_buffers[i].data(),
2428 &first_requests[i]);
2433 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
2434 vertices_to_recv.size());
2435 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
2436 vert_to_recv_it = vertices_to_recv.begin(),
2437 vert_to_recv_end = vertices_to_recv.end();
2438 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2439 ++vert_to_recv_it, ++i)
2441 int source = vert_to_recv_it->first;
2442 const unsigned int n_vertices = vert_to_recv_it->second.size();
2443 const int buffer_size = 2 * n_vertices;
2444 vertices_recv_buffers[i].resize(buffer_size);
2447 ierr = MPI_Recv(vertices_recv_buffers[i].data(),
2459 std::vector<std::vector<char>> cellids_send_buffers(
2460 vertices_to_send.size());
2461 std::vector<MPI_Request> second_requests(vertices_to_send.size());
2462 vert_to_send_it = vertices_to_send.begin();
2463 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
2464 ++vert_to_send_it, ++i)
2466 int destination = vert_to_send_it->first;
2467 const unsigned int n_vertices = vert_to_send_it->second.size();
2468 const int buffer_size = max_cellid_size * n_vertices;
2469 cellids_send_buffers[i].resize(buffer_size);
2472 unsigned int pos = 0;
2473 for (
unsigned int j = 0; j < n_vertices; ++j)
2475 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
2476 for (
unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
2478 if (k < cell_id.size())
2479 cellids_send_buffers[i][pos] = cell_id[k];
2483 cellids_send_buffers[i][pos] =
'-';
2488 ierr = MPI_Isend(cellids_send_buffers[i].data(),
2494 &second_requests[i]);
2499 std::vector<std::vector<char>> cellids_recv_buffers(
2500 vertices_to_recv.size());
2501 vert_to_recv_it = vertices_to_recv.begin();
2502 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2503 ++vert_to_recv_it, ++i)
2505 int source = vert_to_recv_it->first;
2506 const unsigned int n_vertices = vert_to_recv_it->second.size();
2507 const int buffer_size = max_cellid_size * n_vertices;
2508 cellids_recv_buffers[i].resize(buffer_size);
2511 ierr = MPI_Recv(cellids_recv_buffers[i].data(),
2523 vert_to_recv_it = vertices_to_recv.begin();
2524 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2525 ++i, ++vert_to_recv_it)
2527 for (
unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
2529 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
2530 const types::global_vertex_index global_id_recv =
2531 vertices_recv_buffers[i][2 * j + 1];
2532 const std::string cellid_recv(
2533 &cellids_recv_buffers[i][max_cellid_size * j],
2534 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
2536 typename std::set<active_cell_iterator>::iterator
2537 cell_set_it = missing_vert_cells.begin(),
2538 end_cell_set = missing_vert_cells.end();
2539 for (; (found ==
false) && (cell_set_it != end_cell_set);
2542 typename std::set<active_cell_iterator>::iterator
2544 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
2546 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
2547 for (; candidate_cell != end_cell; ++candidate_cell)
2549 std::string current_cellid =
2550 (*candidate_cell)->id().to_string();
2551 current_cellid.resize(max_cellid_size,
'-');
2552 if (current_cellid.compare(cellid_recv) == 0)
2554 local_to_global_vertex_index
2555 [(*candidate_cell)->vertex_index(local_pos_recv)] =
2567 return local_to_global_vertex_index;
2572 template <
int dim,
int spacedim>
2590 const unsigned int index = cell->active_cell_index();
2591 cell_connectivity.
add(index, index);
2592 for (
auto f : cell->face_indices())
2593 if ((cell->at_boundary(f) ==
false) &&
2594 (cell->neighbor(f)->has_children() ==
false))
2596 const unsigned int other_index =
2597 cell->neighbor(f)->active_cell_index();
2598 cell_connectivity.
add(index, other_index);
2599 cell_connectivity.
add(other_index, index);
2606 template <
int dim,
int spacedim>
2612 std::vector<std::vector<unsigned int>> vertex_to_cell(
2616 for (
const unsigned int v : cell->vertex_indices())
2617 vertex_to_cell[cell->vertex_index(v)].push_back(
2618 cell->active_cell_index());
2625 for (
const unsigned int v : cell->vertex_indices())
2626 for (
unsigned int n = 0;
2627 n < vertex_to_cell[cell->vertex_index(v)].size();
2629 cell_connectivity.
add(cell->active_cell_index(),
2630 vertex_to_cell[cell->vertex_index(v)][n]);
2635 template <
int dim,
int spacedim>
2639 const unsigned int level,
2642 std::vector<std::vector<unsigned int>> vertex_to_cell(
2645 triangulation.
begin(level);
2646 cell != triangulation.
end(level);
2649 for (
const unsigned int v : cell->vertex_indices())
2650 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
2654 triangulation.
n_cells(level));
2656 triangulation.
begin(level);
2657 cell != triangulation.
end(level);
2660 for (
const unsigned int v : cell->vertex_indices())
2661 for (
unsigned int n = 0;
2662 n < vertex_to_cell[cell->vertex_index(v)].size();
2664 cell_connectivity.
add(cell->index(),
2665 vertex_to_cell[cell->vertex_index(v)][n]);
2671 template <
int dim,
int spacedim>
2678 &triangulation) ==
nullptr),
2679 ExcMessage(
"Objects of type parallel::distributed::Triangulation " 2680 "are already partitioned implicitly and can not be " 2681 "partitioned again explicitly."));
2683 std::vector<unsigned int> cell_weights;
2686 if (!triangulation.
signals.cell_weight.empty())
2694 if (cell->is_locally_owned())
2695 cell_weights[cell->active_cell_index()] =
2696 triangulation.
signals.cell_weight(
2705 if (
const auto shared_tria =
2709 shared_tria->get_communicator(),
2722 template <
int dim,
int spacedim>
2725 const std::vector<unsigned int> &cell_weights,
2730 &triangulation) ==
nullptr),
2731 ExcMessage(
"Objects of type parallel::distributed::Triangulation " 2732 "are already partitioned implicitly and can not be " 2733 "partitioned again explicitly."));
2737 if (n_partitions == 1)
2740 cell->set_subdomain_id(0);
2754 sp_cell_connectivity.
copy_from(cell_connectivity);
2757 sp_cell_connectivity,
2764 template <
int dim,
int spacedim>
2772 &triangulation) ==
nullptr),
2773 ExcMessage(
"Objects of type parallel::distributed::Triangulation " 2774 "are already partitioned implicitly and can not be " 2775 "partitioned again explicitly."));
2777 std::vector<unsigned int> cell_weights;
2780 if (!triangulation.
signals.cell_weight.empty())
2788 if (cell->is_locally_owned())
2789 cell_weights[cell->active_cell_index()] =
2790 triangulation.
signals.cell_weight(
2799 if (
const auto shared_tria =
2803 shared_tria->get_communicator(),
2810 cell_connection_graph,
2817 template <
int dim,
int spacedim>
2820 const std::vector<unsigned int> &cell_weights,
2826 &triangulation) ==
nullptr),
2827 ExcMessage(
"Objects of type parallel::distributed::Triangulation " 2828 "are already partitioned implicitly and can not be " 2829 "partitioned again explicitly."));
2832 ExcMessage(
"Connectivity graph has wrong size"));
2834 ExcMessage(
"Connectivity graph has wrong size"));
2837 triangulation.
signals.pre_partition();
2840 if (n_partitions == 1)
2843 cell->set_subdomain_id(0);
2851 std::vector<unsigned int> partition_indices(triangulation.
n_active_cells());
2860 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
2872 unsigned int & current_proc_idx,
2873 unsigned int & current_cell_idx,
2875 const unsigned int n_partitions)
2877 if (cell->is_active())
2879 while (current_cell_idx >=
2880 std::floor(static_cast<uint_least64_t>(n_active_cells) *
2881 (current_proc_idx + 1) / n_partitions))
2883 cell->set_subdomain_id(current_proc_idx);
2888 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2898 template <
int dim,
int spacedim>
2902 const bool group_siblings)
2905 &triangulation) ==
nullptr),
2906 ExcMessage(
"Objects of type parallel::distributed::Triangulation " 2907 "are already partitioned implicitly and can not be " 2908 "partitioned again explicitly."));
2912 triangulation.
signals.pre_partition();
2915 if (n_partitions == 1)
2918 cell->set_subdomain_id(0);
2924 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
2925 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
2931 coarse_cell_to_p4est_tree_permutation.resize(triangulation.
n_cells(0));
2933 coarse_cell_to_p4est_tree_permutation);
2935 p4est_tree_to_coarse_cell_permutation =
2938 unsigned int current_proc_idx = 0;
2939 unsigned int current_cell_idx = 0;
2944 for (
unsigned int idx = 0; idx < triangulation.
n_cells(0); ++idx)
2946 const unsigned int coarse_cell_idx =
2947 p4est_tree_to_coarse_cell_permutation[idx];
2949 &triangulation, 0, coarse_cell_idx);
2968 cell = triangulation.
begin(),
2969 endc = triangulation.
end();
2970 for (; cell != endc; ++cell)
2972 if (cell->is_active())
2974 bool all_children_active =
true;
2975 std::map<unsigned int, unsigned int> map_cpu_n_cells;
2976 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2977 if (!cell->child(n)->is_active())
2979 all_children_active =
false;
2983 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
2985 if (!all_children_active)
2988 unsigned int new_owner = cell->child(0)->subdomain_id();
2989 for (std::map<unsigned int, unsigned int>::iterator it =
2990 map_cpu_n_cells.begin();
2991 it != map_cpu_n_cells.end();
2993 if (it->second > map_cpu_n_cells[new_owner])
2994 new_owner = it->first;
2996 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2997 cell->child(n)->set_subdomain_id(new_owner);
3003 template <
int dim,
int spacedim>
3007 unsigned int n_levels = triangulation.
n_levels();
3008 for (
int lvl = n_levels - 1; lvl >= 0; --lvl)
3011 cell = triangulation.
begin(lvl),
3012 endc = triangulation.
end(lvl);
3013 for (; cell != endc; ++cell)
3015 if (cell->is_active())
3016 cell->set_level_subdomain_id(cell->subdomain_id());
3019 Assert(cell->child(0)->level_subdomain_id() !=
3022 cell->set_level_subdomain_id(
3023 cell->child(0)->level_subdomain_id());
3031 template <
int dim,
int spacedim>
3032 std::vector<types::subdomain_id>
3034 const std::vector<CellId> & cell_ids)
3036 std::vector<types::subdomain_id> subdomain_ids;
3037 subdomain_ids.reserve(cell_ids.size());
3041 &triangulation) !=
nullptr)
3046 *parallel_tria =
dynamic_cast< 3050 #ifndef DEAL_II_WITH_P4EST 3054 "You are attempting to use a functionality that is only available " 3055 "if deal.II was configured to use p4est, but cmake did not find a " 3056 "valid p4est library."));
3061 for (
const auto &cell_id : cell_ids)
3064 typename ::internal::p4est::types<dim>::quadrant p4est_cell,
3067 ::internal::p4est::init_coarse_quadrant<dim>(p4est_cell);
3068 for (
const auto &child_index : cell_id.get_child_indices())
3070 ::internal::p4est::init_quadrant_children<dim>(
3071 p4est_cell, p4est_children);
3073 p4est_children[
static_cast<unsigned int>(child_index)];
3079 const_cast<typename ::internal::p4est::types<dim>::forest
3080 *
>(parallel_tria->get_p4est()),
3081 cell_id.get_coarse_cell_id(),
3084 parallel_tria->get_communicator()));
3088 subdomain_ids.push_back(owner);
3098 const std::vector<types::subdomain_id> &true_subdomain_ids_of_cells =
3099 shared_tria->get_true_subdomain_ids_of_cells();
3101 for (
const auto &cell_id : cell_ids)
3103 const unsigned int active_cell_index =
3104 shared_tria->create_cell_iterator(cell_id)->active_cell_index();
3105 subdomain_ids.push_back(
3106 true_subdomain_ids_of_cells[active_cell_index]);
3113 for (
const auto &cell_id : cell_ids)
3115 subdomain_ids.push_back(
3120 return subdomain_ids;
3125 template <
int dim,
int spacedim>
3128 std::vector<types::subdomain_id> & subdomain)
3134 subdomain[cell->active_cell_index()] = cell->subdomain_id();
3139 template <
int dim,
int spacedim>
3145 unsigned int count = 0;
3147 if (cell->subdomain_id() == subdomain)
3155 template <
int dim,
int spacedim>
3160 std::vector<bool> locally_owned_vertices =
3167 if (
const auto *tr =
dynamic_cast< 3171 if (cell->is_artificial() ||
3172 (cell->is_ghost() &&
3173 (cell->subdomain_id() < tr->locally_owned_subdomain())))
3174 for (
const unsigned int v : cell->vertex_indices())
3175 locally_owned_vertices[cell->vertex_index(v)] =
false;
3177 return locally_owned_vertices;
3182 template <
int dim,
int spacedim>
3189 if (!cell->is_artificial())
3190 min_diameter =
std::min(min_diameter, cell->diameter(mapping));
3192 double global_min_diameter = 0;
3194 #ifdef DEAL_II_WITH_MPI 3198 global_min_diameter =
3202 global_min_diameter = min_diameter;
3204 return global_min_diameter;
3209 template <
int dim,
int spacedim>
3214 double max_diameter = 0.;
3216 if (!cell->is_artificial())
3217 max_diameter =
std::max(max_diameter, cell->diameter(mapping));
3219 double global_max_diameter = 0;
3221 #ifdef DEAL_II_WITH_MPI 3225 global_max_diameter =
3229 global_max_diameter = max_diameter;
3231 return global_max_diameter;
3238 namespace FixUpDistortedChildCells
3261 template <
typename Iterator,
int spacedim>
3266 const unsigned int structdim =
3267 Iterator::AccessorType::structure_dimension;
3268 Assert(spacedim == Iterator::AccessorType::dimension,
3274 Assert(object->refinement_case() ==
3282 Tensor<spacedim - structdim, spacedim>
3285 for (
const unsigned int i : object->vertex_indices())
3286 parent_vertices[i] = object->vertex(i);
3289 parent_vertices, parent_alternating_forms);
3291 const Tensor<spacedim - structdim, spacedim>
3292 average_parent_alternating_form =
3293 std::accumulate(parent_alternating_forms,
3294 parent_alternating_forms +
3308 Tensor<spacedim - structdim, spacedim> child_alternating_forms
3312 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3313 for (
const unsigned int i : object->child(c)->vertex_indices())
3314 child_vertices[c][i] = object->child(c)->vertex(i);
3322 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3324 1] = object_mid_point;
3326 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3328 child_vertices[c], child_alternating_forms[c]);
3340 double objective = 0;
3341 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3342 for (
const unsigned int i : object->child(c)->vertex_indices())
3344 (child_alternating_forms[c][i] -
3345 average_parent_alternating_form / std::pow(2., 1. * structdim))
3357 template <
typename Iterator>
3360 const unsigned int f,
3361 std::integral_constant<int, 1>)
3363 return object->vertex(f);
3373 template <
typename Iterator>
3376 const unsigned int f,
3377 std::integral_constant<int, 2>)
3379 return object->line(f)->center();
3389 template <
typename Iterator>
3392 const unsigned int f,
3393 std::integral_constant<int, 3>)
3395 return object->face(f)->center();
3422 template <
typename Iterator>
3426 const unsigned int structdim =
3427 Iterator::AccessorType::structure_dimension;
3429 double diameter =
object->diameter();
3430 for (
const unsigned int f : object->face_indices())
3431 for (
unsigned int e = f + 1;
e <
object->n_faces(); ++
e)
3436 std::integral_constant<int, structdim>())
3438 object,
e, std::integral_constant<int, structdim>())));
3449 template <
typename Iterator>
3453 const unsigned int structdim =
3454 Iterator::AccessorType::structure_dimension;
3455 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
3462 Assert(object->refinement_case() ==
3472 unsigned int iteration = 0;
3477 double initial_delta = 0;
3484 const double step_length = diameter / 4 / (iteration + 1);
3489 for (
unsigned int d = 0;
d < spacedim; ++
d)
3491 const double eps = step_length / 10;
3505 if (gradient.
norm() == 0)
3514 std::min(2 * current_value / (gradient * gradient),
3515 step_length / gradient.norm()) *
3520 const double previous_value = current_value;
3524 initial_delta = (previous_value - current_value);
3527 if ((iteration >= 1) &&
3528 ((previous_value - current_value < 0) ||
3529 (
std::fabs(previous_value - current_value) <
3530 0.001 * initial_delta)))
3535 while (iteration < 20);
3551 double old_min_product, new_min_product;
3556 parent_vertices[i] = object->vertex(i);
3558 Tensor<spacedim - structdim, spacedim>
3561 parent_vertices, parent_alternating_forms);
3567 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3568 for (
const unsigned int i : object->child(c)->vertex_indices())
3569 child_vertices[c][i] = object->child(c)->vertex(i);
3571 Tensor<spacedim - structdim, spacedim> child_alternating_forms
3575 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3577 child_vertices[c], child_alternating_forms[c]);
3580 child_alternating_forms[0][0] * parent_alternating_forms[0];
3581 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3582 for (
const unsigned int i : object->child(c)->vertex_indices())
3583 for (
const unsigned int j : object->vertex_indices())
3584 old_min_product = std::min<double>(old_min_product,
3585 child_alternating_forms[c][i] *
3586 parent_alternating_forms[j]);
3594 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3596 1] = object_mid_point;
3598 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3600 child_vertices[c], child_alternating_forms[c]);
3603 child_alternating_forms[0][0] * parent_alternating_forms[0];
3604 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3605 for (
const unsigned int i : object->child(c)->vertex_indices())
3606 for (
const unsigned int j : object->vertex_indices())
3607 new_min_product = std::min<double>(new_min_product,
3608 child_alternating_forms[c][i] *
3609 parent_alternating_forms[j]);
3617 if (new_min_product >= old_min_product)
3618 object->child(0)->vertex(
3625 return (
std::max(new_min_product, old_min_product) > 0);
3631 template <
int dim,
int spacedim>
3634 const typename ::Triangulation<dim, spacedim>::cell_iterator
3636 std::integral_constant<int, dim>,
3637 std::integral_constant<int, spacedim>)
3647 for (
auto f : cell->face_indices())
3650 Assert(cell->face(f)->refinement_case() ==
3654 bool subface_is_more_refined =
false;
3655 for (
unsigned int g = 0;
3656 g < GeometryInfo<dim>::max_children_per_face;
3658 if (cell->face(f)->child(g)->has_children())
3660 subface_is_more_refined =
true;
3664 if (subface_is_more_refined ==
true)
3675 template <
int dim,
int spacedim>
3683 dim != 1 && spacedim != 1,
3684 "This function is only valid when dim != 1 or spacedim != 1.");
3688 for (
typename std::list<
3697 Assert(!cell->is_active(),
3699 "This function is only valid for a list of cells that " 3700 "have children (i.e., no cell in the list may be active)."));
3704 std::integral_constant<int, dim>(),
3705 std::integral_constant<int, spacedim>());
3712 return unfixable_subset;
3717 template <
int dim,
int spacedim>
3720 const bool reset_boundary_ids)
3723 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
3724 auto m_it = dst_manifold_ids.begin();
3725 for (
const auto b : src_boundary_ids)
3730 const std::vector<types::boundary_id> reset_boundary_id =
3731 reset_boundary_ids ?
3732 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
3742 template <
int dim,
int spacedim>
3745 const std::vector<types::boundary_id> &src_boundary_ids,
3746 const std::vector<types::manifold_id> &dst_manifold_ids,
3748 const std::vector<types::boundary_id> &reset_boundary_ids_)
3751 const auto reset_boundary_ids =
3752 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
3761 for (
auto f : cell->face_indices())
3762 if (cell->face(f)->at_boundary())
3763 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
3765 const auto bid = cell->face(f)->line(
e)->boundary_id();
3766 const unsigned int ind = std::find(src_boundary_ids.begin(),
3767 src_boundary_ids.end(),
3769 src_boundary_ids.begin();
3770 if (ind < src_boundary_ids.size())
3771 cell->face(f)->line(
e)->set_manifold_id(
3772 dst_manifold_ids[ind]);
3777 for (
auto f : cell->face_indices())
3778 if (cell->face(f)->at_boundary())
3780 const auto bid = cell->face(f)->boundary_id();
3781 const unsigned int ind =
3782 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
3783 src_boundary_ids.begin();
3785 if (ind < src_boundary_ids.size())
3788 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
3790 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
3794 for (
unsigned int e = 0;
e < cell->face(f)->n_lines(); ++
e)
3796 const auto bid = cell->face(f)->line(
e)->boundary_id();
3797 const unsigned int ind = std::find(src_boundary_ids.begin(),
3798 src_boundary_ids.end(),
3800 src_boundary_ids.begin();
3801 if (ind < src_boundary_ids.size())
3802 cell->face(f)->line(
e)->set_boundary_id(
3803 reset_boundary_ids[ind]);
3809 template <
int dim,
int spacedim>
3812 const bool compute_face_ids)
3818 for (; cell != endc; ++cell)
3820 cell->set_manifold_id(cell->material_id());
3821 if (compute_face_ids ==
true)
3823 for (
auto f : cell->face_indices())
3825 if (cell->at_boundary(f) ==
false)
3826 cell->face(f)->set_manifold_id(
3828 cell->neighbor(f)->material_id()));
3830 cell->face(f)->set_manifold_id(cell->material_id());
3837 template <
int dim,
int spacedim>
3842 const std::set<types::manifold_id> &)> &disambiguation_function,
3843 bool overwrite_only_flat_manifold_ids)
3848 const unsigned int n_subobjects =
3852 std::vector<std::set<types::manifold_id>> manifold_ids(n_subobjects + 1);
3853 std::vector<unsigned int> backup;
3857 unsigned next_index = 1;
3861 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
3863 if (cell->line(
l)->user_index() == 0)
3866 manifold_ids[next_index].insert(cell->manifold_id());
3867 cell->line(
l)->set_user_index(next_index++);
3870 manifold_ids[cell->line(
l)->user_index()].insert(
3871 cell->manifold_id());
3874 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
3876 if (cell->quad(
l)->user_index() == 0)
3879 manifold_ids[next_index].insert(cell->manifold_id());
3880 cell->quad(
l)->set_user_index(next_index++);
3883 manifold_ids[cell->quad(
l)->user_index()].insert(
3884 cell->manifold_id());
3890 for (
unsigned int l = 0;
l < cell->n_lines(); ++
l)
3892 const auto id = cell->line(
l)->user_index();
3896 if (cell->line(
l)->manifold_id() ==
3898 overwrite_only_flat_manifold_ids ==
false)
3899 cell->line(
l)->set_manifold_id(
3900 disambiguation_function(manifold_ids[
id]));
3901 cell->line(
l)->set_user_index(0);
3905 for (
unsigned int l = 0;
l < cell->n_faces(); ++
l)
3907 const auto id = cell->quad(
l)->user_index();
3911 if (cell->quad(
l)->manifold_id() ==
3913 overwrite_only_flat_manifold_ids ==
false)
3914 cell->quad(
l)->set_manifold_id(
3915 disambiguation_function(manifold_ids[
id]));
3916 cell->quad(
l)->set_user_index(0);
3925 template <
int dim,
int spacedim>
3926 std::pair<unsigned int, double>
3930 double max_ratio = 1;
3931 unsigned int index = 0;
3933 for (
unsigned int i = 0; i < dim; ++i)
3934 for (
unsigned int j = i + 1; j < dim; ++j)
3936 unsigned int ax = i % dim;
3937 unsigned int next_ax = j % dim;
3940 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
3942 if (ratio > max_ratio)
3947 else if (1.0 / ratio > max_ratio)
3949 max_ratio = 1.0 / ratio;
3953 return std::make_pair(index, max_ratio);
3957 template <
int dim,
int spacedim>
3960 const bool isotropic,
3961 const unsigned int max_iterations)
3963 unsigned int iter = 0;
3964 bool continue_refinement =
true;
3966 while (continue_refinement && (iter < max_iterations))
3970 continue_refinement =
false;
3973 for (
const unsigned int j : cell->face_indices())
3974 if (cell->at_boundary(j) ==
false &&
3975 cell->neighbor(j)->has_children())
3979 cell->set_refine_flag();
3980 continue_refinement =
true;
3983 continue_refinement |= cell->flag_for_face_refinement(j);
3990 template <
int dim,
int spacedim>
3993 const double max_ratio,
3994 const unsigned int max_iterations)
3996 unsigned int iter = 0;
3997 bool continue_refinement =
true;
3999 while (continue_refinement && (iter < max_iterations))
4002 continue_refinement =
false;
4005 std::pair<unsigned int, double> info =
4006 GridTools::get_longest_direction<dim, spacedim>(cell);
4007 if (info.second > max_ratio)
4009 cell->set_refine_flag(
4011 continue_refinement =
true;
4019 template <
int dim,
int spacedim>
4022 const double limit_angle_fraction)
4030 "have hanging nodes."));
4033 bool has_cells_with_more_than_dim_faces_on_boundary =
true;
4034 bool has_cells_with_dim_faces_on_boundary =
false;
4036 unsigned int refinement_cycles = 0;
4038 while (has_cells_with_more_than_dim_faces_on_boundary)
4040 has_cells_with_more_than_dim_faces_on_boundary =
false;
4044 unsigned int boundary_face_counter = 0;
4045 for (
auto f : cell->face_indices())
4046 if (cell->face(f)->at_boundary())
4047 boundary_face_counter++;
4048 if (boundary_face_counter > dim)
4050 has_cells_with_more_than_dim_faces_on_boundary =
true;
4053 else if (boundary_face_counter == dim)
4054 has_cells_with_dim_faces_on_boundary =
true;
4056 if (has_cells_with_more_than_dim_faces_on_boundary)
4059 refinement_cycles++;
4063 if (has_cells_with_dim_faces_on_boundary)
4066 refinement_cycles++;
4070 while (refinement_cycles > 0)
4073 cell->set_coarsen_flag();
4075 refinement_cycles--;
4083 std::vector<bool> faces_to_remove(tria.
n_raw_faces(),
false);
4085 std::vector<CellData<dim>> cells_to_add;
4089 const unsigned int v0 = 0, v1 = 1, v2 = (dim > 1 ? 2 : 0),
4090 v3 = (dim > 1 ? 3 : 0);
4094 double angle_fraction = 0;
4100 p0[spacedim > 1 ? 1 : 0] = 1;
4104 if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
4106 p0 = cell->vertex(v0) - cell->vertex(v2);
4107 p1 = cell->vertex(v3) - cell->vertex(v2);
4108 vertex_at_corner = v2;
4110 else if (cell->face(v3)->at_boundary() &&
4111 cell->face(v1)->at_boundary())
4113 p0 = cell->vertex(v2) - cell->vertex(v3);
4114 p1 = cell->vertex(v1) - cell->vertex(v3);
4115 vertex_at_corner = v3;
4117 else if (cell->face(1)->at_boundary() &&
4118 cell->face(2)->at_boundary())
4120 p0 = cell->vertex(v0) - cell->vertex(v1);
4121 p1 = cell->vertex(v3) - cell->vertex(v1);
4122 vertex_at_corner = v1;
4124 else if (cell->face(2)->at_boundary() &&
4125 cell->face(0)->at_boundary())
4127 p0 = cell->vertex(v2) - cell->vertex(v0);
4128 p1 = cell->vertex(v1) - cell->vertex(v0);
4129 vertex_at_corner = v0;
4140 if (angle_fraction > limit_angle_fraction)
4142 auto flags_removal = [&](
unsigned int f1,
4145 unsigned int n2) ->
void {
4146 cells_to_remove[cell->active_cell_index()] =
true;
4147 cells_to_remove[cell->neighbor(n1)->active_cell_index()] =
true;
4148 cells_to_remove[cell->neighbor(n2)->active_cell_index()] =
true;
4150 faces_to_remove[cell->face(f1)->index()] =
true;
4151 faces_to_remove[cell->face(f2)->index()] =
true;
4153 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] =
true;
4154 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] =
true;
4157 auto cell_creation = [&](
const unsigned int vv0,
4158 const unsigned int vv1,
4159 const unsigned int f0,
4160 const unsigned int f1,
4162 const unsigned int n0,
4163 const unsigned int v0n0,
4164 const unsigned int v1n0,
4166 const unsigned int n1,
4167 const unsigned int v0n1,
4168 const unsigned int v1n1) {
4172 c1.
vertices[v0] = cell->vertex_index(vv0);
4173 c1.
vertices[v1] = cell->vertex_index(vv1);
4174 c1.
vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
4175 c1.
vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
4180 c2.
vertices[v0] = cell->vertex_index(vv0);
4181 c2.
vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
4182 c2.
vertices[v2] = cell->vertex_index(vv1);
4183 c2.
vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
4188 l1.
vertices[0] = cell->vertex_index(vv0);
4189 l1.
vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
4195 l2.
vertices[0] = cell->vertex_index(vv0);
4196 l2.
vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
4202 cells_to_add.push_back(c1);
4203 cells_to_add.push_back(c2);
4208 switch (vertex_at_corner)
4211 flags_removal(0, 2, 3, 1);
4212 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
4215 flags_removal(1, 2, 3, 0);
4216 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
4219 flags_removal(3, 0, 1, 2);
4220 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
4223 flags_removal(3, 1, 0, 2);
4224 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
4237 if (cells_to_add.size() == 0)
4239 while (refinement_cycles > 0)
4242 cell->set_coarsen_flag();
4244 refinement_cycles--;
4252 if (cells_to_remove[cell->active_cell_index()] ==
false)
4255 for (
const unsigned int v : cell->vertex_indices())
4256 c.
vertices[v] = cell->vertex_index(v);
4259 cells_to_add.push_back(c);
4267 for (; face != endf; ++face)
4268 if ((face->at_boundary() ||
4270 faces_to_remove[face->index()] ==
false)
4272 for (
unsigned int l = 0;
l < face->n_lines(); ++
l)
4277 for (
const unsigned int v : face->vertex_indices())
4278 line.
vertices[v] = face->vertex_index(v);
4284 for (
const unsigned int v : face->line(
l)->vertex_indices())
4285 line.
vertices[v] = face->line(
l)->vertex_index(v);
4294 for (
const unsigned int v : face->vertex_indices())
4295 quad.
vertices[v] = face->vertex_index(v);
4303 subcelldata_to_add);
4308 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
4327 template <
int dim,
int spacedim>
4330 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
4331 std::vector<std::vector<Point<dim>>>,
4332 std::vector<std::vector<unsigned int>>>
4344 auto &cells = std::get<0>(cqmp);
4345 auto &qpoints = std::get<1>(cqmp);
4346 auto &maps = std::get<2>(cqmp);
4347 auto &missing_points = std::get<3>(cqmp);
4351 ExcPointNotFound<spacedim>(points[missing_points[0]]));
4353 (void)missing_points;
4355 return std::make_tuple(std::move(cells),
4362 template <
int dim,
int spacedim>
4365 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
4366 std::vector<std::vector<Point<dim>>>,
4367 std::vector<std::vector<unsigned int>>,
4368 std::vector<unsigned int>>
4379 const unsigned int np = points.size();
4381 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
4383 std::vector<std::vector<Point<dim>>> qpoints_out;
4384 std::vector<std::vector<unsigned int>> maps_out;
4385 std::vector<unsigned int> missing_points_out;
4389 return std::make_tuple(std::move(cells_out),
4390 std::move(qpoints_out),
4391 std::move(maps_out),
4392 std::move(missing_points_out));
4398 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4403 unsigned int points_checked = 0;
4415 catch (
const GridTools::ExcPointNotFound<dim> &)
4417 missing_points_out.emplace_back(0);
4426 std::pair<BoundingBox<spacedim>,
4431 int cell_candidate_idx = -1;
4435 bool use_try =
false;
4437 while (!found && points_checked < np)
4440 b_tree.query(boost::geometry::index::intersects(points[points_checked]),
4441 std::back_inserter(box_cell));
4444 cell_candidate_idx = -1;
4445 for (
unsigned int i = 0; i < box_cell.size(); ++i)
4448 if (!box_cell[i].
second->is_artificial())
4449 cell_candidate_idx = i;
4453 if (cell_candidate_idx != -1 &&
4454 (!box_cell[i].
second->is_locally_owned() ||
4455 box_cell[i].second->at_boundary()))
4459 if (cell_candidate_idx != -1)
4464 if (cell_candidate_idx != -1)
4472 points[points_checked],
4473 box_cell[cell_candidate_idx].
second);
4476 catch (
const GridTools::ExcPointNotFound<dim> &)
4478 missing_points_out.emplace_back(points_checked);
4485 points[points_checked],
4486 box_cell[cell_candidate_idx].
second);
4495 cache, points[points_checked]);
4500 cell_candidate_idx = box_cell.size();
4502 std::make_pair(my_pair.first->bounding_box(), my_pair.first));
4504 catch (
const GridTools::ExcPointNotFound<dim> &)
4506 missing_points_out.emplace_back(points_checked);
4517 cells_out.emplace_back(my_pair.first);
4518 qpoints_out.emplace_back(1, my_pair.second);
4519 maps_out.emplace_back(1, points_checked - 1);
4523 if (np == qpoints_out.size())
4524 return std::make_tuple(std::move(cells_out),
4525 std::move(qpoints_out),
4526 std::move(maps_out),
4527 std::move(missing_points_out));
4530 for (
unsigned int p = points_checked; p < np; ++p)
4533 if (cell_candidate_idx != -1)
4534 if (!box_cell[cell_candidate_idx].
first.point_inside(points[p]))
4536 cell_candidate_idx = -1;
4539 if (cell_candidate_idx == -1)
4543 b_tree.query(boost::geometry::index::intersects(points[p]),
4544 std::back_inserter(box_cell));
4547 cell_candidate_idx = -1;
4548 for (
unsigned int i = 0; i < box_cell.size(); ++i)
4551 if (!box_cell[i].
second->is_artificial())
4552 cell_candidate_idx = i;
4556 if (cell_candidate_idx != -1 &&
4557 (!box_cell[i].
second->is_locally_owned() ||
4558 box_cell[i].second->at_boundary()))
4562 if (cell_candidate_idx != -1)
4567 if (cell_candidate_idx == -1)
4578 cell_candidate_idx = box_cell.size();
4580 std::make_pair(my_pair.first->bounding_box(), my_pair.first));
4582 catch (
const GridTools::ExcPointNotFound<dim> &)
4584 missing_points_out.emplace_back(p);
4596 cache, points[p], box_cell[cell_candidate_idx].
second);
4598 catch (
const GridTools::ExcPointNotFound<dim> &)
4600 missing_points_out.push_back(p);
4607 cache, points[p], box_cell[cell_candidate_idx].
second);
4612 if (my_pair.first != box_cell[cell_candidate_idx].second)
4614 for (
unsigned int i = 0; i < box_cell.size(); ++i)
4616 if (my_pair.first == box_cell[i].second)
4618 cell_candidate_idx = i;
4623 if (my_pair.first != box_cell[cell_candidate_idx].second)
4627 cell_candidate_idx = box_cell.size();
4629 std::make_pair(my_pair.first->bounding_box(),
4638 if (my_pair.first == cells_out.back())
4641 qpoints_out.back().emplace_back(my_pair.second);
4642 maps_out.back().emplace_back(p);
4647 typename std::vector<typename Triangulation<dim, spacedim>::
4648 active_cell_iterator>::iterator cells_it =
4649 std::find(cells_out.begin(), cells_out.end() - 1, my_pair.first);
4651 if (cells_it == cells_out.end() - 1)
4654 cells_out.emplace_back(my_pair.first);
4655 qpoints_out.emplace_back(1, my_pair.second);
4656 maps_out.emplace_back(1, p);
4662 unsigned int current_cell = cells_it - cells_out.begin();
4663 qpoints_out[current_cell].emplace_back(my_pair.second);
4664 maps_out[current_cell].emplace_back(p);
4670 Assert(cells_out.size() == maps_out.size(),
4673 Assert(cells_out.size() == qpoints_out.size(),
4677 unsigned int c = cells_out.size();
4678 unsigned int qps = 0;
4684 for (
unsigned int n = 0; n < c; ++n)
4686 Assert(qpoints_out[n].size() == maps_out[n].size(),
4688 qps += qpoints_out[n].size();
4691 Assert(qps + missing_points_out.size() == np,
4695 return std::make_tuple(std::move(cells_out),
4696 std::move(qpoints_out),
4697 std::move(maps_out),
4698 std::move(missing_points_out));
4706 namespace DistributedComputePointLocations
4709 template <
int dim,
int spacedim>
4718 return k->active_cell_index();
4727 template <
int dim,
int spacedim>
4730 std::pair<std::vector<Point<dim>>, std::vector<unsigned int>>,
4735 const unsigned int n_points = points.size();
4738 typename Triangulation<dim, spacedim>::active_cell_iterator,
4739 std::pair<std::vector<Point<dim>>, std::vector<unsigned int>>,
4745 return cell_qpoint_map;
4748 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4750 point_and_reference_location;
4752 unsigned int counter = 0;
4754 while (counter < n_points)
4757 unsigned int i = counter;
4760 point_and_reference_location =
4766 if (counter == n_points)
4767 return cell_qpoint_map;
4770 auto last_cell = cell_qpoint_map.emplace(std::make_pair(
4771 point_and_reference_location.first,
4773 std::vector<
Point<dim>>{point_and_reference_location.second},
4774 std::vector<unsigned int>{counter - 1})));
4778 return cell_qpoint_map;
4781 point_and_reference_location.first->center();
4782 double cell_diameter = point_and_reference_location.first->diameter() *
4786 for (
unsigned int p = counter; p < n_points; ++p)
4790 if (cell_center.
distance(points[p]) < cell_diameter)
4793 point_and_reference_location =
4795 cache, points[p], last_cell.first->first);
4804 point_and_reference_location =
4812 if (last_cell.first->first == point_and_reference_location.first)
4814 last_cell.first->second.first.emplace_back(
4815 point_and_reference_location.second);
4816 last_cell.first->second.second.emplace_back(p);
4821 last_cell = cell_qpoint_map.emplace(
4822 std::make_pair(point_and_reference_location.first,
4825 point_and_reference_location.second},
4826 std::vector<unsigned int>{p})));
4828 if (last_cell.second ==
false)
4831 last_cell.first->second.first.emplace_back(
4832 point_and_reference_location.second);
4833 last_cell.first->second.second.emplace_back(p);
4838 cell_center = point_and_reference_location.first->center();
4840 point_and_reference_location.first->diameter() *
4847 unsigned int inserted_points = 0;
4852 for (
const auto &map_entry : cell_qpoint_map)
4854 Assert(map_entry.second.second.size() ==
4855 map_entry.second.first.size(),
4857 map_entry.second.first.size()));
4858 inserted_points += map_entry.second.second.size();
4861 return cell_qpoint_map;
4873 template <
int dim,
int spacedim>
4877 typename Triangulation<dim, spacedim>::active_cell_iterator> &cells,
4878 const std::vector<std::vector<
Point<dim>>> & qpoints,
4879 const std::vector<std::vector<unsigned int>> & maps,
4881 const unsigned int rank,
4883 typename Triangulation<dim, spacedim>::active_cell_iterator,
4885 std::vector<unsigned int>,
4887 std::vector<unsigned int>>,
4891 for (
unsigned int c = 0; c < cells.size(); ++c)
4894 auto current_c = point_locations.emplace(
4895 std::make_pair(cells[c],
4896 std::make_tuple(qpoints[c],
4899 std::vector<unsigned int>(
4900 points[c].size(), rank))));
4903 if (current_c.second ==
false)
4906 auto &cell_qpts = std::get<0>(current_c.first->second);
4907 auto &cell_maps = std::get<1>(current_c.first->second);
4908 auto &cell_pts = std::get<2>(current_c.first->second);
4909 auto &cell_ranks = std::get<3>(current_c.first->second);
4911 cell_qpts.insert(cell_qpts.end(),
4914 cell_maps.insert(cell_maps.end(),
4917 cell_pts.insert(cell_pts.end(),
4920 std::vector<unsigned int> ranks_tmp(points[c].size(), rank);
4921 cell_ranks.insert(cell_ranks.end(),
4942 template <
int dim,
int spacedim>
4947 const std::vector<unsigned int> & local_points_idx,
4949 typename Triangulation<dim, spacedim>::active_cell_iterator,
4951 std::vector<unsigned int>,
4953 std::vector<unsigned int>>,
4955 std::map<
unsigned int,
4956 std::tuple<std::vector<CellId>,
4958 std::vector<std::vector<unsigned int>>,
4960 & ghost_cell_locations,
4961 std::vector<unsigned int> &found_location_indices)
4963 auto point_location_data =
4965 cache, local_points);
4969 for (
const auto &cell_tuples : point_location_data)
4971 auto &cell = cell_tuples.first;
4972 auto &q_loc = std::get<0>(cell_tuples.second);
4973 auto &indices_loc = std::get<1>(cell_tuples.second);
4976 if (cell->is_locally_owned())
4978 std::vector<Point<spacedim>> cell_points(indices_loc.size());
4979 std::vector<unsigned int> cell_points_idx(indices_loc.size());
4980 for (
unsigned int i = 0; i < indices_loc.size(); ++i)
4983 cell_points[i] = local_points[indices_loc[i]];
4988 cell_points_idx[i] = local_points_idx[indices_loc[i]];
4989 found_location_indices.emplace_back(
4990 local_points_idx[indices_loc[i]]);
4992 locally_owned_locations.emplace(
4993 std::make_pair(cell,
4994 std::make_tuple(q_loc,
4997 std::vector<unsigned int>(
4999 cell->subdomain_id()))));
5002 else if (cell->is_ghost())
5004 std::vector<Point<spacedim>> cell_points(indices_loc.size());
5005 std::vector<unsigned int> cell_points_idx(indices_loc.size());
5006 for (
unsigned int i = 0; i < indices_loc.size(); ++i)
5008 cell_points[i] = local_points[indices_loc[i]];
5009 cell_points_idx[i] = local_points_idx[indices_loc[i]];
5010 found_location_indices.emplace_back(
5011 local_points_idx[indices_loc[i]]);
5017 auto &map_tuple_owner =
5018 ghost_cell_locations[cell->subdomain_id()];
5020 std::get<0>(map_tuple_owner).emplace_back(cell->id());
5021 std::get<1>(map_tuple_owner).emplace_back(q_loc);
5022 std::get<2>(map_tuple_owner).emplace_back(cell_points_idx);
5023 std::get<3>(map_tuple_owner).emplace_back(cell_points);
5037 template <
int dim,
int spacedim>
5044 &received_point_locations,
5046 typename Triangulation<dim, spacedim>::active_cell_iterator,
5048 std::vector<unsigned int>,
5050 std::vector<unsigned int>>,
5052 const bool check_owned)
5056 for (
const auto &rank_and_points : received_point_locations)
5059 const auto &received_process = rank_and_points.first;
5060 const auto &received_points = rank_and_points.second.first;
5061 const auto &received_map = rank_and_points.second.second;
5066 typename Triangulation<dim, spacedim>::active_cell_iterator>
5068 std::vector<std::vector<Point<dim>>> in_qpoints;
5069 std::vector<std::vector<unsigned int>> in_maps;
5070 std::vector<std::vector<Point<spacedim>>> in_points;
5072 const auto computed_point_locations =
5075 for (
const auto &map_c_pt_idx : computed_point_locations)
5078 const auto &proc_cell = map_c_pt_idx.first;
5079 const auto &proc_qpoints = map_c_pt_idx.second.first;
5080 const auto &proc_maps = map_c_pt_idx.second.second;
5084 if (check_owned ==
false || proc_cell->is_locally_owned())
5086 in_cell.emplace_back(proc_cell);
5087 in_qpoints.emplace_back(proc_qpoints);
5089 unsigned int loc_size = proc_qpoints.size();
5090 std::vector<unsigned int> cell_maps(loc_size);
5091 std::vector<Point<spacedim>> cell_points(loc_size);
5092 for (
unsigned int pt = 0; pt < loc_size; ++pt)
5094 cell_maps[pt] = received_map[proc_maps[pt]];
5095 cell_points[pt] = received_points[proc_maps[pt]];
5097 in_maps.emplace_back(cell_maps);
5098 in_points.emplace_back(cell_points);
5117 template <
int dim,
int spacedim>
5120 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5121 std::vector<std::vector<Point<dim>>>,
5122 std::vector<std::vector<unsigned int>>,
5123 std::vector<std::vector<Point<spacedim>>>,
5124 std::vector<std::vector<unsigned int>>>
5133 #ifndef DEAL_II_WITH_MPI 5136 (void)global_bboxes;
5139 "GridTools::distributed_compute_point_locations() requires MPI."));
5141 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5142 std::vector<std::vector<Point<dim>>>,
5143 std::vector<std::vector<unsigned int>>,
5144 std::vector<std::vector<Point<spacedim>>>,
5145 std::vector<std::vector<unsigned int>>>
5150 const auto &tria_mpi =
5158 "GridTools::distributed_compute_point_locations() requires a parallel triangulation."));
5159 auto mpi_communicator = tria_mpi->get_communicator();
5162 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
5163 std::vector<std::vector<Point<dim>>>,
5164 std::vector<std::vector<unsigned int>>,
5165 std::vector<std::vector<Point<spacedim>>>,
5166 std::vector<std::vector<unsigned int>>>
5172 std::tuple<std::vector<Point<dim>>,
5173 std::vector<unsigned int>,
5174 std::vector<Point<spacedim>>,
5175 std::vector<unsigned int>>,
5181 const unsigned int my_rank =
5185 std::tuple<std::vector<std::vector<unsigned int>>,
5186 std::map<unsigned int, unsigned int>,
5187 std::map<unsigned int, std::vector<unsigned int>>>
5192 const auto & guess_loc_idx = std::get<0>(guessed_points)[my_rank];
5193 const unsigned int n_local_guess = guess_loc_idx.size();
5196 std::vector<Point<spacedim>> guess_local_points(n_local_guess);
5197 for (
unsigned int i = 0; i < n_local_guess; ++i)
5198 guess_local_points[i] = local_points[guess_loc_idx[i]];
5201 std::map<
unsigned int,
5202 std::tuple<std::vector<CellId>,
5203 std::vector<std::vector<Point<dim>>>,
5204 std::vector<std::vector<unsigned int>>,
5205 std::vector<std::vector<Point<spacedim>>>>>
5210 std::vector<unsigned int> found_point_indices;
5215 compute_and_classify_points<dim, spacedim>,
5221 found_point_indices);
5225 const auto ¬_locally_owned_idx = std::get<1>(guessed_points);
5226 std::map<
unsigned int,
5227 std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>>
5228 not_locally_owned_points;
5230 for (
const auto &indices : not_locally_owned_idx)
5231 if (indices.second != my_rank)
5234 auto &points_to_send = not_locally_owned_points[indices.second];
5236 points_to_send.first.emplace_back(local_points[indices.first]);
5237 points_to_send.second.emplace_back(indices.first);
5241 auto received_points =
5244 compute_locations_task.
join();
5249 merge_received_point_locations<dim, spacedim>,
5256 auto received_ghost_points =
5263 std::map<
unsigned int,
5264 std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>>
5269 const std::map<unsigned int, std::vector<unsigned int>>
5270 &points_to_probable_owners = std::get<2>(guessed_points);
5277 std::sort(found_point_indices.begin(), found_point_indices.end());
5279 for (
const auto &probable_owners : points_to_probable_owners)
5281 const auto &point_idx = probable_owners.first;
5282 const auto &probable_owner_ranks = probable_owners.second;
5283 if (!std::binary_search(found_point_indices.begin(),
5284 found_point_indices.end(),
5287 for (
const unsigned int probable_owner_rank : probable_owner_ranks)
5288 if (probable_owner_rank != my_rank)
5291 auto &points_to_send = uncertain_points[probable_owner_rank];
5292 points_to_send.first.emplace_back(local_points[point_idx]);
5293 points_to_send.second.emplace_back(point_idx);
5298 const auto received_uncertain_points =
5301 merge_locally_owned_points_task.
join();
5304 for (
const auto &received_ghost_point : received_ghost_points)
5307 const auto & cell_ids = std::get<0>(received_ghost_point.second);
5308 const unsigned int n_cells = cell_ids.size();
5309 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
5311 for (
unsigned int c = 0; c <
n_cells; ++c)
5317 std::get<1>(received_ghost_point.second),
5318 std::get<2>(received_ghost_point.second),
5319 std::get<3>(received_ghost_point.second),
5320 received_ghost_point.first,
5327 cache, received_uncertain_points, found_points,
true);
5331 const unsigned int size_output = found_points.size();
5332 auto & out_cells = std::get<0>(output_tuple);
5333 auto & out_qpoints = std::get<1>(output_tuple);
5334 auto & out_maps = std::get<2>(output_tuple);
5335 auto & out_points = std::get<3>(output_tuple);
5336 auto & out_ranks = std::get<4>(output_tuple);
5338 out_cells.resize(size_output);
5339 out_qpoints.resize(size_output);
5340 out_maps.resize(size_output);
5341 out_points.resize(size_output);
5342 out_ranks.resize(size_output);
5345 for (
const auto &cell_and_data : found_points)
5347 out_cells[c] = cell_and_data.first;
5348 out_qpoints[c] = std::get<0>(cell_and_data.second);
5349 out_maps[c] = std::get<1>(cell_and_data.second);
5350 out_points[c] = std::get<2>(cell_and_data.second);
5351 out_ranks[c] = std::get<3>(cell_and_data.second);
5355 return output_tuple;
5360 template <
int dim,
int spacedim>
5361 std::map<unsigned int, Point<spacedim>>
5365 std::map<unsigned int, Point<spacedim>> result;
5368 if (!cell->is_artificial())
5371 for (
unsigned int i = 0; i < vs.size(); ++i)
5372 result[cell->vertex_index(i)] = vs[i];
5379 template <
int spacedim>
5384 auto id_and_v = std::min_element(
5387 [&](
const std::pair<const unsigned int, Point<spacedim>> &p1,
5388 const std::pair<const unsigned int, Point<spacedim>> &p2) ->
bool {
5389 return p1.second.
distance(p) < p2.second.distance(p);
5391 return id_and_v->first;
5395 template <
int dim,
int spacedim>
5396 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
5403 const std::vector<bool> &marked_vertices,
5404 const double tolerance)
5409 const auto &vertex_to_cell_centers =
5417 vertex_to_cell_centers,
5420 used_vertices_rtree,
5424 template <
int spacedim>
5425 std::vector<std::vector<BoundingBox<spacedim>>>
5430 #ifndef DEAL_II_WITH_MPI 5432 (void)mpi_communicator;
5435 "GridTools::exchange_local_bounding_boxes() requires MPI."));
5439 unsigned int n_bboxes = local_bboxes.size();
5441 int n_local_data = 2 * spacedim * n_bboxes;
5443 std::vector<double> loc_data_array(n_local_data);
5444 for (
unsigned int i = 0; i < n_bboxes; ++i)
5445 for (
unsigned int d = 0;
d < spacedim; ++
d)
5448 loc_data_array[2 * i * spacedim +
d] =
5449 local_bboxes[i].get_boundary_points().first[
d];
5450 loc_data_array[2 * i * spacedim + spacedim +
d] =
5451 local_bboxes[i].get_boundary_points().second[
d];
5458 std::vector<int> size_all_data(n_procs);
5461 int ierr = MPI_Allgather(&n_local_data,
5464 size_all_data.data(),
5472 std::vector<int> rdispls(n_procs);
5474 for (
unsigned int i = 1; i < n_procs; ++i)
5475 rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
5479 std::vector<double> data_array(rdispls.back() + size_all_data.back());
5481 ierr = MPI_Allgatherv(loc_data_array.data(),
5485 size_all_data.data(),
5492 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
5493 unsigned int begin_idx = 0;
5494 for (
unsigned int i = 0; i < n_procs; ++i)
5497 unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
5498 global_bboxes[i].resize(n_bbox_i);
5499 for (
unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
5502 for (
unsigned int d = 0;
d < spacedim; ++
d)
5504 p1[
d] = data_array[begin_idx + 2 * bbox * spacedim +
d];
5506 data_array[begin_idx + 2 * bbox * spacedim + spacedim +
d];
5509 global_bboxes[i][bbox] = loc_bbox;
5512 begin_idx += size_all_data[i];
5514 return global_bboxes;
5515 #endif // DEAL_II_WITH_MPI 5520 template <
int spacedim>
5526 #ifndef DEAL_II_WITH_MPI 5527 (void)mpi_communicator;
5529 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>> boxes_index(
5530 local_description.size());
5532 for (
unsigned int i = 0; i < local_description.size(); ++i)
5533 boxes_index[i] = std::make_pair(local_description[i], 0u);
5537 const std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes =
5541 const unsigned int n_procs =
5545 std::vector<unsigned int> bboxes_position(n_procs);
5547 unsigned int tot_bboxes = 0;
5548 for (
const auto &process_bboxes : global_bboxes)
5549 tot_bboxes += process_bboxes.size();
5552 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>>
5554 flat_global_bboxes.reserve(tot_bboxes);
5555 unsigned int process_index = 0;
5556 for (
const auto &process_bboxes : global_bboxes)
5559 std::vector<std::pair<BoundingBox<spacedim>,
unsigned int>>
5560 boxes_and_indices(process_bboxes.size());
5563 for (
unsigned int i = 0; i < process_bboxes.size(); ++i)
5564 boxes_and_indices[i] =
5565 std::make_pair(process_bboxes[i], process_index);
5567 flat_global_bboxes.insert(flat_global_bboxes.end(),
5568 boxes_and_indices.begin(),
5569 boxes_and_indices.end());
5576 return RTree<std::pair<BoundingBox<spacedim>,
unsigned int>>(
5577 flat_global_bboxes.begin(), flat_global_bboxes.end());
5578 #endif // DEAL_II_WITH_MPI 5583 template <
int dim,
int spacedim>
5587 std::map<
unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
5588 std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group)
5593 static const int lookup_table_2d[2][2] =
5600 static const int lookup_table_3d[2][2][2][4] =
5622 if (pair.first.first->level() != pair.second.first.first->level())
5625 const auto face_a = pair.first.first->face(pair.first.second);
5627 pair.second.first.first->face(pair.second.first.second);
5628 const auto mask = pair.second.second;
5633 for (
unsigned int i = 0; i < face_a->n_vertices(); ++i)
5635 const bool face_orientation = mask[0];
5636 const bool face_flip = mask[1];
5637 const bool face_rotation = mask[2];
5647 j = lookup_table_2d[face_flip][i];
5650 j = lookup_table_3d[face_orientation][face_flip]
5658 const auto vertex_a = face_a->vertex_index(i);
5659 const auto vertex_b = face_b->vertex_index(j);
5660 unsigned int temp =
std::min(vertex_a, vertex_b);
5662 auto it_a = vertex_to_coinciding_vertex_group.find(vertex_a);
5663 if (it_a != vertex_to_coinciding_vertex_group.end())
5664 temp =
std::min(temp, it_a->second);
5666 auto it_b = vertex_to_coinciding_vertex_group.find(vertex_b);
5667 if (it_b != vertex_to_coinciding_vertex_group.end())
5668 temp =
std::min(temp, it_b->second);
5670 if (it_a != vertex_to_coinciding_vertex_group.end())
5671 it_a->second = temp;
5673 vertex_to_coinciding_vertex_group[vertex_a] = temp;
5675 if (it_b != vertex_to_coinciding_vertex_group.end())
5676 it_b->second = temp;
5678 vertex_to_coinciding_vertex_group[vertex_b] = temp;
5684 for (
auto &p : vertex_to_coinciding_vertex_group)
5686 if (p.first == p.second)
5688 unsigned int temp = p.second;
5689 while (temp != vertex_to_coinciding_vertex_group[temp])
5690 temp = vertex_to_coinciding_vertex_group[temp];
5696 for (
auto p : vertex_to_coinciding_vertex_group)
5697 coinciding_vertex_groups[p.second] = {};
5699 for (
auto p : vertex_to_coinciding_vertex_group)
5700 coinciding_vertex_groups[p.second].push_back(p.first);
5706 template <
int dim,
int spacedim>
5707 std::map<unsigned int, std::set<::types::subdomain_id>>
5717 std::map<unsigned int, std::vector<unsigned int>> coinciding_vertex_groups;
5718 std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
5721 coinciding_vertex_groups,
5722 vertex_to_coinciding_vertex_group);
5725 std::vector<bool> vertex_of_own_cell(tria.
n_vertices(),
false);
5727 if (cell->is_locally_owned())
5728 for (
const unsigned int v : cell->vertex_indices())
5729 vertex_of_own_cell[cell->vertex_index(v)] =
true;
5733 std::map<unsigned int, std::set<types::subdomain_id>> result;
5737 if (cell->is_ghost())
5742 for (
const unsigned int v : cell->vertex_indices())
5745 if (vertex_of_own_cell[cell->vertex_index(v)])
5746 result[cell->vertex_index(v)].insert(owner);
5749 auto coinciding_vertex_group =
5750 vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
5751 if (coinciding_vertex_group !=
5752 vertex_to_coinciding_vertex_group.end())
5753 for (
auto coinciding_vertex :
5754 coinciding_vertex_groups[coinciding_vertex_group->second])
5755 if (vertex_of_own_cell[coinciding_vertex])
5756 result[coinciding_vertex].insert(owner);
5767 #include "grid_tools.inst"
void map_boundary_to_manifold_ids(const std::vector< types::boundary_id > &src_boundary_ids, const std::vector< types::manifold_id > &dst_manifold_ids, Triangulation< dim, spacedim > &tria, const std::vector< types::boundary_id > &reset_boundary_ids={})
std::vector< CellData< 1 > > boundary_lines
Transformed quadrature weights.
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
unsigned int n_active_cells() const
unsigned int n_vertices() const
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
static void reorder_cells(std::vector< CellData< dim >> &original_cells, const bool use_new_style_ordering=false)
const types::manifold_id flat_manifold_id
static const unsigned int invalid_unsigned_int
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
#define AssertDimension(dim1, dim2)
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
active_face_iterator begin_active_face() const
void create_laplace_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, SparseMatrix< double > &matrix, const Function< spacedim > *const a=nullptr, const AffineConstraints< double > &constraints=AffineConstraints< double >())
typename IteratorSelector::line_iterator line_iterator
virtual bool has_hanging_nodes() const
unsigned int n_cells() const
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Task< RT > new_task(const std::function< RT()> &function)
void add(const size_type i, const size_type j)
IteratorRange< active_cell_iterator > active_cell_iterators() const
#define AssertIndexRange(index, range)
std::vector< unsigned int > vertices
LinearOperator< Range, Domain, Payload > constrained_linear_operator(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop)
active_cell_iterator begin_active(const unsigned int level=0) const
#define AssertThrow(cond, exc)
types::boundary_id boundary_id
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
const DerivativeForm< 1, dim, spacedim > & jacobian(const unsigned int quadrature_point) const
cell_iterator begin(const unsigned int level=0) const
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
unsigned int n_levels() const
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
static double distance_to_unit_cell(const Point< dim > &p)
cell_iterator end() const
virtual void execute_coarsening_and_refinement()
RTree< typename LeafTypeIterator::value_type, IndexType, IndexableGetter > pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
bool check_consistency(const unsigned int dim) const
unsigned int subdomain_id
T sum(const T &t, const MPI_Comm &mpi_communicator)
Expression acos(const Expression &x)
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
#define Assert(cond, exc)
IteratorRange< active_cell_iterator > active_cell_iterators() const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
types::global_dof_index n_dofs() const
Abstract base class for mapping classes.
std::list< typename Triangulation< dim, spacedim >::cell_iterator > distorted_cells
unsigned int n_quads() const
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim - dim, spacedim >(&forms)[vertices_per_cell])
void save_user_indices(std::vector< unsigned int > &v) const
types::material_id material_id
const std::vector< Point< spacedim > > & get_vertices() const
#define DEAL_II_NAMESPACE_CLOSE
void load_user_indices(const std::vector< unsigned int > &v)
unsigned int n_lines() const
VectorType::value_type * end(VectorType &V)
std::vector< Integer > invert_permutation(const std::vector< Integer > &permutation)
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
uint64_t global_vertex_index
Expression fabs(const Expression &x)
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
cell_iterator create_cell_iterator(const CellId &cell_id) const
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
PackagedOperation< Range > constrained_right_hand_side(const AffineConstraints< typename Range::value_type > &constraints, const LinearOperator< Range, Domain, Payload > &linop, const Range &right_hand_side)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
unsigned int size() const
virtual const MPI_Comm & get_communicator() const
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
number singular_value(const size_type i) const
types::manifold_id manifold_id
unsigned int n_raw_faces() const
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
const types::subdomain_id artificial_subdomain_id
__global__ void set(Number *val, const Number s, const size_type N)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
void advance(std::tuple< I1, I2 > &t, const unsigned int n)
void distribute(VectorType &vec) const
#define AssertThrowMPI(error_code)
static constexpr double PI
const std::vector< bool > & get_used_vertices() const
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NAMESPACE_OPEN
VectorType::value_type * begin(VectorType &V)
T min(const T &t, const MPI_Comm &mpi_communicator)
std::vector< CellData< 2 > > boundary_quads
numbers::NumberTraits< Number >::real_type square() const
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void refine_global(const unsigned int times=1)