16 #include <deal.II/base/mpi.h> 17 #include <deal.II/base/mpi.templates.h> 18 #include <deal.II/base/quadrature_lib.h> 19 #include <deal.II/base/thread_management.h> 21 #include <deal.II/distributed/shared_tria.h> 22 #include <deal.II/distributed/tria.h> 24 #include <deal.II/dofs/dof_accessor.h> 25 #include <deal.II/dofs/dof_handler.h> 26 #include <deal.II/dofs/dof_tools.h> 28 #include <deal.II/fe/fe_nothing.h> 29 #include <deal.II/fe/fe_q.h> 30 #include <deal.II/fe/fe_values.h> 31 #include <deal.II/fe/mapping_q.h> 32 #include <deal.II/fe/mapping_q1.h> 33 #include <deal.II/fe/mapping_q_generic.h> 35 #include <deal.II/grid/filtered_iterator.h> 36 #include <deal.II/grid/grid_reordering.h> 37 #include <deal.II/grid/grid_tools.h> 38 #include <deal.II/grid/grid_tools_cache.h> 39 #include <deal.II/grid/manifold.h> 40 #include <deal.II/grid/tria.h> 41 #include <deal.II/grid/tria_accessor.h> 42 #include <deal.II/grid/tria_iterator.h> 44 #include <deal.II/lac/dynamic_sparsity_pattern.h> 45 #include <deal.II/lac/filtered_matrix.h> 46 #include <deal.II/lac/precondition.h> 47 #include <deal.II/lac/solver_cg.h> 48 #include <deal.II/lac/sparse_matrix.h> 49 #include <deal.II/lac/sparsity_pattern.h> 50 #include <deal.II/lac/sparsity_tools.h> 51 #include <deal.II/lac/vector.h> 52 #include <deal.II/lac/vector_memory.h> 54 #include <deal.II/numerics/matrix_tools.h> 56 #include <boost/random/mersenne_twister.hpp> 57 #include <boost/random/uniform_real_distribution.hpp> 66 #include <unordered_map> 68 DEAL_II_NAMESPACE_OPEN
73 template <
int dim,
int spacedim>
83 #if defined(DEAL_II_WITH_P4EST) && defined(DEBUG) 97 const std::vector<Point<spacedim>> &vertices = tria.
get_vertices();
98 std::vector<bool> boundary_vertices(vertices.size(),
false);
104 for (; cell != endc; ++cell)
105 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
107 if (cell->face(face)->at_boundary())
108 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face;
110 boundary_vertices[cell->face(face)->vertex_index(i)] =
true;
114 double max_distance_sqr = 0;
115 std::vector<bool>::const_iterator pi = boundary_vertices.begin();
116 const unsigned int N = boundary_vertices.size();
117 for (
unsigned int i = 0; i < N; ++i, ++pi)
119 std::vector<bool>::const_iterator pj = pi + 1;
120 for (
unsigned int j = i + 1; j < N; ++j, ++pj)
121 if ((*pi ==
true) && (*pj ==
true) &&
122 ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
123 max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
126 return std::sqrt(max_distance_sqr);
131 template <
int dim,
int spacedim>
137 unsigned int mapping_degree = 1;
140 mapping_degree = p->get_degree();
141 else if (
const auto *p =
143 mapping_degree = p->get_degree();
146 const QGauss<dim> quadrature_formula(mapping_degree + 1);
147 const unsigned int n_q_points = quadrature_formula.
size();
160 endc = triangulation.
end();
162 double local_volume = 0;
165 for (; cell != endc; ++cell)
166 if (cell->is_locally_owned())
169 for (
unsigned int q = 0; q < n_q_points; ++q)
170 local_volume += fe_values.
JxW(q);
173 double global_volume = 0;
175 #ifdef DEAL_II_WITH_MPI 183 global_volume = local_volume;
185 return global_volume;
193 const std::vector<Point<1>> &all_vertices,
196 return all_vertices[vertex_indices[1]][0] -
197 all_vertices[vertex_indices[0]][0];
205 const std::vector<Point<3>> &all_vertices,
214 const double x[8] = {all_vertices[vertex_indices[0]](0),
215 all_vertices[vertex_indices[1]](0),
216 all_vertices[vertex_indices[2]](0),
217 all_vertices[vertex_indices[3]](0),
218 all_vertices[vertex_indices[4]](0),
219 all_vertices[vertex_indices[5]](0),
220 all_vertices[vertex_indices[6]](0),
221 all_vertices[vertex_indices[7]](0)};
222 const double y[8] = {all_vertices[vertex_indices[0]](1),
223 all_vertices[vertex_indices[1]](1),
224 all_vertices[vertex_indices[2]](1),
225 all_vertices[vertex_indices[3]](1),
226 all_vertices[vertex_indices[4]](1),
227 all_vertices[vertex_indices[5]](1),
228 all_vertices[vertex_indices[6]](1),
229 all_vertices[vertex_indices[7]](1)};
230 const double z[8] = {all_vertices[vertex_indices[0]](2),
231 all_vertices[vertex_indices[1]](2),
232 all_vertices[vertex_indices[2]](2),
233 all_vertices[vertex_indices[3]](2),
234 all_vertices[vertex_indices[4]](2),
235 all_vertices[vertex_indices[5]](2),
236 all_vertices[vertex_indices[6]](2),
237 all_vertices[vertex_indices[7]](2)};
276 const double t3 = y[3] * x[2];
277 const double t5 = z[1] * x[5];
278 const double t9 = z[3] * x[2];
279 const double t11 = x[1] * y[0];
280 const double t14 = x[4] * y[0];
281 const double t18 = x[5] * y[7];
282 const double t20 = y[1] * x[3];
283 const double t22 = y[5] * x[4];
284 const double t26 = z[7] * x[6];
285 const double t28 = x[0] * y[4];
287 z[3] * x[1] * y[2] + t3 * z[1] - t5 * y[7] + y[7] * x[4] * z[6] +
288 t9 * y[6] - t11 * z[4] - t5 * y[3] - t14 * z[2] + z[1] * x[4] * y[0] -
289 t18 * z[3] + t20 * z[0] - t22 * z[0] - y[0] * x[5] * z[4] - t26 * y[3] +
290 t28 * z[2] - t9 * y[1] - y[1] * x[4] * z[0] - t11 * z[5];
291 const double t37 = y[1] * x[0];
292 const double t44 = x[1] * y[5];
293 const double t46 = z[1] * x[0];
294 const double t49 = x[0] * y[2];
295 const double t52 = y[5] * x[7];
296 const double t54 = x[3] * y[7];
297 const double t56 = x[2] * z[0];
298 const double t58 = x[3] * y[2];
299 const double t64 = -x[6] * y[4] * z[2] - t37 * z[2] + t18 * z[6] -
300 x[3] * y[6] * z[2] + t11 * z[2] + t5 * y[0] +
301 t44 * z[4] - t46 * y[4] - t20 * z[7] - t49 * z[6] -
302 t22 * z[1] + t52 * z[3] - t54 * z[2] - t56 * y[4] -
303 t58 * z[0] + y[1] * x[2] * z[0] + t9 * y[7] + t37 * z[4];
304 const double t66 = x[1] * y[7];
305 const double t68 = y[0] * x[6];
306 const double t70 = x[7] * y[6];
307 const double t73 = z[5] * x[4];
308 const double t76 = x[6] * y[7];
309 const double t90 = x[4] * z[0];
310 const double t92 = x[1] * y[3];
311 const double t95 = -t66 * z[3] - t68 * z[2] - t70 * z[2] + t26 * y[5] -
312 t73 * y[6] - t14 * z[6] + t76 * z[2] - t3 * z[6] +
313 x[6] * y[2] * z[4] - z[3] * x[6] * y[2] + t26 * y[4] -
314 t44 * z[3] - x[1] * y[2] * z[0] + x[5] * y[6] * z[4] +
315 t54 * z[5] + t90 * y[2] - t92 * z[2] + t46 * y[2];
316 const double t102 = x[2] * y[0];
317 const double t107 = y[3] * x[7];
318 const double t114 = x[0] * y[6];
320 y[0] * x[3] * z[2] - z[7] * x[5] * y[6] - x[2] * y[6] * z[4] +
321 t102 * z[6] - t52 * z[6] + x[2] * y[4] * z[6] - t107 * z[5] - t54 * z[6] +
322 t58 * z[6] - x[7] * y[4] * z[6] + t37 * z[5] - t114 * z[4] + t102 * z[4] -
323 z[1] * x[2] * y[0] + t28 * z[6] - y[5] * x[6] * z[4] -
324 z[5] * x[1] * y[4] - t73 * y[7];
325 const double t129 = z[0] * x[6];
326 const double t133 = y[1] * x[7];
327 const double t145 = y[1] * x[5];
328 const double t156 = t90 * y[6] - t129 * y[4] + z[7] * x[2] * y[6] -
329 t133 * z[5] + x[5] * y[3] * z[7] - t26 * y[2] -
330 t70 * z[3] + t46 * y[3] + z[5] * x[7] * y[4] +
331 z[7] * x[3] * y[6] - t49 * z[4] + t145 * z[7] -
332 x[2] * y[7] * z[6] + t70 * z[5] + t66 * z[5] -
333 z[7] * x[4] * y[6] + t18 * z[4] + x[1] * y[4] * z[0];
334 const double t160 = x[5] * y[4];
335 const double t165 = z[1] * x[7];
336 const double t178 = z[1] * x[3];
338 t107 * z[6] + t22 * z[7] + t76 * z[3] + t160 * z[1] - x[4] * y[2] * z[6] +
339 t70 * z[4] + t165 * y[5] + x[7] * y[2] * z[6] - t76 * z[5] - t76 * z[4] +
340 t133 * z[3] - t58 * z[1] + y[5] * x[0] * z[4] + t114 * z[2] - t3 * z[7] +
341 t20 * z[2] + t178 * y[7] + t129 * y[2];
342 const double t207 = t92 * z[7] + t22 * z[6] + z[3] * x[0] * y[2] -
343 x[0] * y[3] * z[2] - z[3] * x[7] * y[2] - t165 * y[3] -
344 t9 * y[0] + t58 * z[7] + y[3] * x[6] * z[2] +
345 t107 * z[2] + t73 * y[0] - x[3] * y[5] * z[7] +
346 t3 * z[0] - t56 * y[6] - z[5] * x[0] * y[4] +
347 t73 * y[1] - t160 * z[6] + t160 * z[0];
348 const double t228 = -t44 * z[7] + z[5] * x[6] * y[4] - t52 * z[4] -
349 t145 * z[4] + t68 * z[4] + t92 * z[5] - t92 * z[0] +
350 t11 * z[3] + t44 * z[0] + t178 * y[5] - t46 * y[5] -
351 t178 * y[0] - t145 * z[0] - t20 * z[5] - t37 * z[3] -
352 t160 * z[7] + t145 * z[3] + x[4] * y[6] * z[2];
354 return (t34 + t64 + t95 + t125 + t156 + t181 + t207 + t228) / 12.;
362 const std::vector<Point<2>> &all_vertices,
401 const double x[4] = {all_vertices[vertex_indices[0]](0),
402 all_vertices[vertex_indices[1]](0),
403 all_vertices[vertex_indices[2]](0),
404 all_vertices[vertex_indices[3]](0)};
406 const double y[4] = {all_vertices[vertex_indices[0]](1),
407 all_vertices[vertex_indices[1]](1),
408 all_vertices[vertex_indices[2]](1),
409 all_vertices[vertex_indices[3]](1)};
411 return (-x[1] * y[0] + x[1] * y[3] + y[0] * x[2] + x[0] * y[1] -
412 x[0] * y[2] - y[1] * x[3] - x[2] * y[3] + x[3] * y[2]) /
418 template <
int dim,
int spacedim>
424 const auto predicate = [](
const iterator &) {
return true; };
427 tria, std::function<
bool(
const iterator &)>(predicate));
432 template <
int dim,
int spacedim>
441 "Invalid SubCellData supplied according to ::check_consistency(). " 442 "This is caused by data containing objects for the wrong dimension."));
445 std::vector<bool> vertex_used(vertices.size(),
false);
446 for (
unsigned int c = 0; c < cells.size(); ++c)
447 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
449 Assert(cells[c].vertices[v] < vertices.size(),
450 ExcMessage(
"Invalid vertex index encountered! cells[" +
454 " is invalid, because only " +
456 " vertices were supplied."));
457 vertex_used[cells[c].vertices[v]] =
true;
464 std::vector<unsigned int> new_vertex_numbers(vertices.size(),
466 unsigned int next_free_number = 0;
467 for (
unsigned int i = 0; i < vertices.size(); ++i)
468 if (vertex_used[i] ==
true)
470 new_vertex_numbers[i] = next_free_number;
475 for (
unsigned int c = 0; c < cells.size(); ++c)
476 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
477 cells[c].vertices[v] = new_vertex_numbers[cells[c].vertices[v]];
480 for (
unsigned int c = 0; c < subcelldata.
boundary_lines.size(); ++c)
481 for (
unsigned int v = 0; v < GeometryInfo<1>::vertices_per_cell; ++v)
484 new_vertex_numbers.size(),
486 "Invalid vertex index in subcelldata.boundary_lines. " 487 "subcelldata.boundary_lines[" +
492 " is invalid, because only " +
494 " vertices were supplied."));
499 for (
unsigned int c = 0; c < subcelldata.
boundary_quads.size(); ++c)
500 for (
unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
503 new_vertex_numbers.size(),
505 "Invalid vertex index in subcelldata.boundary_quads. " 506 "subcelldata.boundary_quads[" +
511 " is invalid, because only " +
513 " vertices were supplied."));
521 std::vector<Point<spacedim>> tmp;
522 tmp.reserve(std::count(vertex_used.begin(), vertex_used.end(),
true));
523 for (
unsigned int v = 0; v < vertices.size(); ++v)
524 if (vertex_used[v] ==
true)
525 tmp.push_back(vertices[v]);
531 template <
int dim,
int spacedim>
536 std::vector<unsigned int> & considered_vertices,
542 std::vector<unsigned int> new_vertex_numbers(vertices.size());
543 for (
unsigned int i = 0; i < vertices.size(); ++i)
544 new_vertex_numbers[i] = i;
548 if (considered_vertices.size() == 0)
549 considered_vertices = new_vertex_numbers;
557 for (
unsigned int i = 0; i < considered_vertices.size(); ++i)
560 if (new_vertex_numbers[considered_vertices[i]] !=
561 considered_vertices[i])
572 for (
unsigned int j = i + 1; j < considered_vertices.size(); ++j)
575 for (
unsigned int d = 0; d < spacedim; ++d)
576 equal &= (fabs(vertices[considered_vertices[j]](d) -
577 vertices[considered_vertices[i]](d)) < tol);
580 new_vertex_numbers[considered_vertices[j]] =
581 considered_vertices[i];
590 for (
auto &cell : cells)
591 for (
auto &vertex_index : cell.vertices)
592 vertex_index = new_vertex_numbers[vertex_index];
594 for (
auto &vertex_index : quad.vertices)
595 vertex_index = new_vertex_numbers[vertex_index];
597 for (
auto &vertex_index : line.vertices)
598 vertex_index = new_vertex_numbers[vertex_index];
608 template <
int spacedim>
632 explicit Rotate2d(
const double angle)
638 return Point<2>(std::cos(angle) * p(0) - std::sin(angle) * p(1),
639 std::sin(angle) * p(0) + std::cos(angle) * p(1));
650 Rotate3d(
const double angle,
const unsigned int axis)
660 std::cos(angle) * p(1) - std::sin(angle) * p(2),
661 std::sin(angle) * p(1) + std::cos(angle) * p(2));
663 return Point<3>(std::cos(angle) * p(0) + std::sin(angle) * p(2),
665 -std::sin(angle) * p(0) + std::cos(angle) * p(2));
667 return Point<3>(std::cos(angle) * p(0) - std::sin(angle) * p(1),
668 std::sin(angle) * p(0) + std::cos(angle) * p(1),
674 const unsigned int axis;
677 template <
int spacedim>
681 explicit Scale(
const double factor)
696 template <
int dim,
int spacedim>
701 transform(Shift<spacedim>(shift_vector), triangulation);
709 transform(Rotate2d(angle), triangulation);
715 const unsigned int axis,
720 transform(Rotate3d(angle, axis), triangulation);
723 template <
int dim,
int spacedim>
729 transform(Scale<spacedim>(scaling_factor), triangulation);
742 const std::map<types::global_dof_index, double> &fixed_dofs,
745 const unsigned int n_dofs = S.
n();
759 solver.
solve(SF, u, f, PF);
783 const bool solve_for_absolute_positions)
808 std::map<types::global_dof_index, double> fixed_dofs[dim];
809 typename std::map<unsigned int, Point<dim>>::const_iterator map_end =
814 endc = dof_handler.
end();
815 for (; cell != endc; ++cell)
819 for (
unsigned int vertex_no = 0;
820 vertex_no < GeometryInfo<dim>::vertices_per_cell;
823 const unsigned int vertex_index = cell->vertex_index(vertex_no);
824 const Point<dim> & vertex_point = cell->vertex(vertex_no);
826 const typename std::map<unsigned int, Point<dim>>::const_iterator
827 map_iter = new_points.find(vertex_index);
829 if (map_iter != map_end)
830 for (
unsigned int i = 0; i < dim; ++i)
831 fixed_dofs[i].insert(std::pair<types::global_dof_index, double>(
832 cell->vertex_dof_index(vertex_no, 0),
833 (solve_for_absolute_positions ?
834 map_iter->second(i) :
835 map_iter->second(i) - vertex_point[i])));
841 for (
unsigned int i = 0; i < dim; ++i)
842 us[i].reinit(dof_handler.
n_dofs());
846 for (
unsigned int i = 0; i < dim; ++i)
852 std::vector<bool> vertex_touched(triangulation.
n_vertices(),
false);
853 for (cell = dof_handler.
begin_active(); cell != endc; ++cell)
854 for (
unsigned int vertex_no = 0;
855 vertex_no < GeometryInfo<dim>::vertices_per_cell;
857 if (vertex_touched[cell->vertex_index(vertex_no)] ==
false)
862 cell->vertex_dof_index(vertex_no, 0);
863 for (
unsigned int i = 0; i < dim; ++i)
864 if (solve_for_absolute_positions)
865 v(i) = us[i](dof_index);
867 v(i) += us[i](dof_index);
869 vertex_touched[cell->vertex_index(vertex_no)] =
true;
873 template <
int dim,
int spacedim>
874 std::map<unsigned int, Point<spacedim>>
877 std::map<unsigned int, Point<spacedim>> vertex_map;
881 for (; cell != endc; ++cell)
883 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
887 if (face->at_boundary())
889 for (
unsigned j = 0; j < GeometryInfo<dim>::vertices_per_face;
893 const unsigned int vertex_index = face->vertex_index(j);
894 vertex_map[vertex_index] = vertex;
906 template <
int dim,
int spacedim>
910 const bool keep_boundary)
927 double almost_infinite_length = 0;
929 triangulation.
begin(0);
930 cell != triangulation.
end(0);
932 almost_infinite_length += cell->diameter();
934 std::vector<double> minimal_length(triangulation.
n_vertices(),
935 almost_infinite_length);
938 std::vector<bool> at_boundary(keep_boundary ? triangulation.
n_vertices() :
943 const bool is_parallel_shared =
945 &triangulation) !=
nullptr);
948 cell != triangulation.
end();
950 if (is_parallel_shared || cell->is_locally_owned())
954 for (
unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell;
957 const typename Triangulation<dim, spacedim>::line_iterator
958 line = cell->line(i);
960 if (keep_boundary && line->at_boundary())
962 at_boundary[line->vertex_index(0)] =
true;
963 at_boundary[line->vertex_index(1)] =
true;
966 minimal_length[line->vertex_index(0)] =
967 std::min(line->diameter(),
968 minimal_length[line->vertex_index(0)]);
969 minimal_length[line->vertex_index(1)] =
970 std::min(line->diameter(),
971 minimal_length[line->vertex_index(1)]);
977 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
978 if (cell->at_boundary(vertex) ==
true)
979 at_boundary[cell->vertex_index(vertex)] =
true;
981 minimal_length[cell->vertex_index(0)] =
982 std::min(cell->diameter(),
983 minimal_length[cell->vertex_index(0)]);
984 minimal_length[cell->vertex_index(1)] =
985 std::min(cell->diameter(),
986 minimal_length[cell->vertex_index(1)]);
995 boost::random::mt19937 rng;
996 boost::random::uniform_real_distribution<> uniform_distribution(-1, 1);
1001 *distributed_triangulation =
1005 const std::vector<bool> locally_owned_vertices =
1007 std::vector<bool> vertex_moved(triangulation.
n_vertices(),
false);
1012 cell != triangulation.
end();
1014 if (cell->is_locally_owned())
1016 for (
unsigned int vertex_no = 0;
1017 vertex_no < GeometryInfo<dim>::vertices_per_cell;
1020 const unsigned global_vertex_no =
1021 cell->vertex_index(vertex_no);
1026 if ((keep_boundary && at_boundary[global_vertex_no]) ||
1027 vertex_moved[global_vertex_no] ||
1028 !locally_owned_vertices[global_vertex_no])
1033 for (
unsigned int d = 0; d < spacedim; ++d)
1034 shift_vector(d) = uniform_distribution(rng);
1036 shift_vector *= factor * minimal_length[global_vertex_no] /
1037 std::sqrt(shift_vector.
square());
1040 cell->vertex(vertex_no) += shift_vector;
1041 vertex_moved[global_vertex_no] =
true;
1045 #ifdef DEAL_II_WITH_P4EST 1046 distributed_triangulation->communicate_locally_moved_vertices(
1047 locally_owned_vertices);
1049 (void)distributed_triangulation;
1060 const unsigned int n_vertices = triangulation.
n_vertices();
1061 std::vector<Point<spacedim>> new_vertex_locations(n_vertices);
1062 const std::vector<Point<spacedim>> &old_vertex_locations =
1065 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1069 if (keep_boundary && at_boundary[vertex])
1070 new_vertex_locations[vertex] = old_vertex_locations[vertex];
1075 for (
unsigned int d = 0; d < spacedim; ++d)
1076 shift_vector(d) = uniform_distribution(rng);
1078 shift_vector *= factor * minimal_length[vertex] /
1079 std::sqrt(shift_vector.
square());
1082 new_vertex_locations[vertex] =
1083 old_vertex_locations[vertex] + shift_vector;
1090 cell != triangulation.
end();
1092 for (
unsigned int vertex_no = 0;
1093 vertex_no < GeometryInfo<dim>::vertices_per_cell;
1095 cell->vertex(vertex_no) =
1096 new_vertex_locations[cell->vertex_index(vertex_no)];
1111 endc = triangulation.
end();
1112 for (; cell != endc; ++cell)
1113 if (!cell->is_artificial())
1114 for (
unsigned int face = 0;
1115 face < GeometryInfo<dim>::faces_per_cell;
1117 if (cell->face(face)->has_children() &&
1118 !cell->face(face)->at_boundary())
1122 cell->face(face)->child(0)->vertex(1) =
1123 (cell->face(face)->vertex(0) +
1124 cell->face(face)->vertex(1)) /
1128 cell->face(face)->child(0)->vertex(1) =
1129 .5 * (cell->face(face)->vertex(0) +
1130 cell->face(face)->vertex(1));
1131 cell->face(face)->child(0)->vertex(2) =
1132 .5 * (cell->face(face)->vertex(0) +
1133 cell->face(face)->vertex(2));
1134 cell->face(face)->child(1)->vertex(3) =
1135 .5 * (cell->face(face)->vertex(1) +
1136 cell->face(face)->vertex(3));
1137 cell->face(face)->child(2)->vertex(3) =
1138 .5 * (cell->face(face)->vertex(2) +
1139 cell->face(face)->vertex(3));
1142 cell->face(face)->child(0)->vertex(3) =
1143 .25 * (cell->face(face)->vertex(0) +
1144 cell->face(face)->vertex(1) +
1145 cell->face(face)->vertex(2) +
1146 cell->face(face)->vertex(3));
1154 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1158 const std::vector<bool> & marked_vertices)
1166 const std::vector<Point<spacedim>> &vertices = tria.
get_vertices();
1169 marked_vertices.size() == 0,
1171 marked_vertices.size()));
1181 marked_vertices.size() == 0 ||
1182 std::equal(marked_vertices.begin(),
1183 marked_vertices.end(),
1185 [](
bool p,
bool q) {
return !p || q; }),
1187 "marked_vertices should be a subset of used vertices in the triangulation " 1188 "but marked_vertices contains one or more vertices that are not used vertices!"));
1196 const std::vector<bool> &used = (marked_vertices.size() == 0) ?
1202 std::vector<bool>::const_iterator first =
1203 std::find(used.begin(), used.end(),
true);
1209 unsigned int best_vertex = std::distance(used.begin(), first);
1210 double best_dist = (p - vertices[best_vertex]).norm_square();
1214 for (
unsigned int j = best_vertex + 1; j < vertices.size(); j++)
1217 double dist = (p - vertices[j]).norm_square();
1218 if (dist < best_dist)
1230 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1233 const MeshType<dim, spacedim> &mesh,
1235 const std::vector<bool> & marked_vertices)
1250 marked_vertices.size() == 0,
1252 marked_vertices.size()));
1262 marked_vertices.size() == 0 ||
1263 std::equal(marked_vertices.begin(),
1264 marked_vertices.end(),
1266 [](
bool p,
bool q) {
return !p || q; }),
1268 "marked_vertices should be a subset of used vertices in the triangulation " 1269 "but marked_vertices contains one or more vertices that are not used vertices!"));
1272 if (marked_vertices.size())
1273 for (
auto it = vertices.begin(); it != vertices.end();)
1275 if (marked_vertices[it->first] ==
false)
1277 vertices.erase(it++);
1290 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1292 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1295 typename ::internal::
1296 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1299 const unsigned int vertex)
1304 Assert(vertex < mesh.get_triangulation().n_vertices(),
1305 ExcIndexRange(0, mesh.get_triangulation().n_vertices(), vertex));
1306 Assert(mesh.get_triangulation().get_used_vertices()[vertex],
1312 std::set<typename ::internal::
1313 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1316 typename ::internal::
1317 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1318 cell = mesh.begin_active(),
1354 for (; cell != endc; ++cell)
1356 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; v++)
1357 if (cell->vertex_index(v) == vertex)
1362 adjacent_cells.insert(cell);
1369 for (
unsigned int vface = 0; vface < dim; vface++)
1371 const unsigned int face =
1374 if (!cell->at_boundary(face) &&
1375 cell->neighbor(face)->active())
1387 adjacent_cells.insert(cell->neighbor(face));
1398 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_cell; ++e)
1399 if (cell->line(e)->has_children())
1403 if (cell->line(e)->child(0)->vertex_index(1) == vertex)
1405 adjacent_cells.insert(cell);
1427 typename ::internal::
1428 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
1429 adjacent_cells.begin(), adjacent_cells.end());
1434 template <
int dim,
int spacedim>
1435 std::vector<std::vector<Tensor<1, spacedim>>>
1442 const std::vector<Point<spacedim>> &vertices = mesh.
get_vertices();
1443 const unsigned int n_vertices = vertex_to_cells.size();
1448 std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers(
1450 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
1453 const unsigned int n_neighbor_cells = vertex_to_cells[vertex].size();
1454 vertex_to_cell_centers[vertex].resize(n_neighbor_cells);
1456 typename std::set<typename Triangulation<dim, spacedim>::
1457 active_cell_iterator>::iterator it =
1458 vertex_to_cells[vertex].begin();
1459 for (
unsigned int cell = 0; cell < n_neighbor_cells; ++cell, ++it)
1461 vertex_to_cell_centers[vertex][cell] =
1462 (*it)->center() - vertices[vertex];
1463 vertex_to_cell_centers[vertex][cell] /=
1464 vertex_to_cell_centers[vertex][cell].norm();
1467 return vertex_to_cell_centers;
1473 template <
int spacedim>
1475 compare_point_association(
1476 const unsigned int a,
1477 const unsigned int b,
1481 const double scalar_product_a = center_directions[a] * point_direction;
1482 const double scalar_product_b = center_directions[b] * point_direction;
1487 return (scalar_product_a > scalar_product_b);
1491 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1493 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1495 std::pair<typename ::internal::
1496 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1501 const MeshType<dim, spacedim> &mesh,
1504 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
1507 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint,
1508 const std::vector<bool> &marked_vertices)
1510 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1515 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1517 cell_and_position_approx;
1519 bool found_cell =
false;
1520 bool approx_cell =
false;
1522 unsigned int closest_vertex_index = 0;
1524 auto current_cell = cell_hint;
1526 while (found_cell ==
false)
1532 const unsigned int closest_vertex =
1533 find_closest_vertex_of_cell<dim, spacedim>(current_cell, p);
1534 vertex_to_point = p - current_cell->vertex(closest_vertex);
1535 closest_vertex_index = current_cell->vertex_index(closest_vertex);
1539 closest_vertex_index =
1541 vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
1544 const double vertex_point_norm = vertex_to_point.
norm();
1545 if (vertex_point_norm > 0)
1546 vertex_to_point /= vertex_point_norm;
1548 const unsigned int n_neighbor_cells =
1549 vertex_to_cells[closest_vertex_index].size();
1552 std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
1554 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
1555 neighbor_permutation[i] = i;
1557 auto comp = [&](
const unsigned int a,
const unsigned int b) ->
bool {
1558 return compare_point_association<spacedim>(
1562 vertex_to_cell_centers[closest_vertex_index]);
1565 std::sort(neighbor_permutation.begin(),
1566 neighbor_permutation.end(),
1572 double best_distance = 1e-10;
1576 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
1580 auto cell = vertex_to_cells[closest_vertex_index].begin();
1581 std::advance(cell, neighbor_permutation[i]);
1586 cell_and_position.first = *cell;
1587 cell_and_position.second = p_unit;
1589 approx_cell =
false;
1597 if (dist < best_distance)
1599 best_distance = dist;
1600 cell_and_position_approx.first = *cell;
1601 cell_and_position_approx.second = p_unit;
1609 if (found_cell ==
true)
1610 return cell_and_position;
1611 else if (approx_cell ==
true)
1612 return cell_and_position_approx;
1626 ExcPointNotFound<spacedim>(p));
1628 current_cell =
typename MeshType<dim, spacedim>::active_cell_iterator();
1630 return cell_and_position;
1635 template <
int dim,
int spacedim>
1642 unsigned int closest_vertex = 0;
1644 for (
unsigned int v = 1; v < GeometryInfo<dim>::vertices_per_cell; ++v)
1646 const double vertex_distance =
1648 if (vertex_distance < minimum_distance)
1651 minimum_distance = vertex_distance;
1654 return closest_vertex;
1661 namespace BoundingBoxPredicate
1663 template <
class MeshType>
1664 std::tuple<BoundingBox<MeshType::space_dimension>,
bool>
1665 compute_cell_predicate_bounding_box(
1666 const typename MeshType::cell_iterator &parent_cell,
1667 const std::function<
1668 bool(
const typename MeshType::active_cell_iterator &)> &predicate)
1670 bool has_predicate =
1672 std::vector<typename MeshType::active_cell_iterator> active_cells;
1673 if (parent_cell->active())
1674 active_cells = {parent_cell};
1678 active_cells = get_active_child_cells<MeshType>(parent_cell);
1680 const unsigned int spacedim = MeshType::space_dimension;
1684 while (i < active_cells.size() && !predicate(active_cells[i]))
1688 if (active_cells.size() == 0 || i == active_cells.size())
1691 return std::make_tuple(bbox, has_predicate);
1698 for (; i < active_cells.size(); ++i)
1699 if (predicate(active_cells[i]))
1700 for (
unsigned int v = 0;
1701 v < GeometryInfo<spacedim>::vertices_per_cell;
1703 for (
unsigned int d = 0; d < spacedim; ++d)
1705 minp[d] = std::min(minp[d], active_cells[i]->vertex(v)[d]);
1706 maxp[d] = std::max(maxp[d], active_cells[i]->vertex(v)[d]);
1709 has_predicate =
true;
1711 return std::make_tuple(bbox, has_predicate);
1718 template <
class MeshType>
1719 std::vector<BoundingBox<MeshType::space_dimension>>
1721 const MeshType &mesh,
1722 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1724 const unsigned int refinement_level,
1725 const bool allow_merge,
1726 const unsigned int max_boxes)
1733 refinement_level <= mesh.n_levels(),
1735 "Error: refinement level is higher then total levels in the triangulation!"));
1737 const unsigned int spacedim = MeshType::space_dimension;
1738 std::vector<BoundingBox<spacedim>> bounding_boxes;
1742 for (
unsigned int i = 0; i < refinement_level; ++i)
1743 for (
const typename MeshType::cell_iterator &cell :
1744 mesh.active_cell_iterators_on_level(i))
1746 bool has_predicate =
false;
1748 std::tie(bbox, has_predicate) =
1749 internal::BoundingBoxPredicate::compute_cell_predicate_bounding_box<
1750 MeshType>(cell, predicate);
1752 bounding_boxes.push_back(bbox);
1756 for (
const typename MeshType::cell_iterator &cell :
1757 mesh.cell_iterators_on_level(refinement_level))
1759 bool has_predicate =
false;
1761 std::tie(bbox, has_predicate) =
1762 internal::BoundingBoxPredicate::compute_cell_predicate_bounding_box<
1763 MeshType>(cell, predicate);
1765 bounding_boxes.push_back(bbox);
1770 return bounding_boxes;
1776 std::vector<unsigned int> merged_boxes_idx;
1777 bool found_neighbors =
true;
1782 while (found_neighbors)
1784 found_neighbors =
false;
1785 for (
unsigned int i = 0; i < bounding_boxes.size() - 1; ++i)
1787 if (std::find(merged_boxes_idx.begin(),
1788 merged_boxes_idx.end(),
1789 i) == merged_boxes_idx.end())
1790 for (
unsigned int j = i + 1; j < bounding_boxes.size(); ++j)
1791 if (std::find(merged_boxes_idx.begin(),
1792 merged_boxes_idx.end(),
1793 j) == merged_boxes_idx.end() &&
1794 bounding_boxes[i].get_neighbor_type(
1795 bounding_boxes[j]) ==
1796 NeighborType::mergeable_neighbors)
1798 bounding_boxes[i].merge_with(bounding_boxes[j]);
1799 merged_boxes_idx.push_back(j);
1800 found_neighbors =
true;
1806 std::vector<BoundingBox<spacedim>> merged_b_boxes;
1807 for (
unsigned int i = 0; i < bounding_boxes.size(); ++i)
1808 if (std::find(merged_boxes_idx.begin(), merged_boxes_idx.end(), i) ==
1809 merged_boxes_idx.end())
1810 merged_b_boxes.push_back(bounding_boxes[i]);
1815 if ((merged_b_boxes.size() > max_boxes) && (spacedim > 1))
1817 std::vector<double> volumes;
1818 for (
unsigned int i = 0; i < merged_b_boxes.size(); ++i)
1819 volumes.push_back(merged_b_boxes[i].volume());
1821 while (merged_b_boxes.size() > max_boxes)
1823 unsigned int min_idx =
1824 std::min_element(volumes.begin(), volumes.end()) -
1826 volumes.erase(volumes.begin() + min_idx);
1828 bool not_removed =
true;
1829 for (
unsigned int i = 0;
1830 i < merged_b_boxes.size() && not_removed;
1835 if (i != min_idx && (merged_b_boxes[i].get_neighbor_type(
1836 merged_b_boxes[min_idx]) ==
1837 NeighborType::attached_neighbors ||
1838 merged_b_boxes[i].get_neighbor_type(
1839 merged_b_boxes[min_idx]) ==
1840 NeighborType::mergeable_neighbors))
1842 merged_b_boxes[i].merge_with(merged_b_boxes[min_idx]);
1843 merged_b_boxes.erase(merged_b_boxes.begin() + min_idx);
1844 not_removed =
false;
1847 ExcMessage(
"Error: couldn't merge bounding boxes!"));
1850 Assert(merged_b_boxes.size() <= max_boxes,
1852 "Error: couldn't reach target number of bounding boxes!"));
1853 return merged_b_boxes;
1859 template <
int spacedim>
1861 std::tuple<std::vector<std::vector<unsigned int>>,
1862 std::map<unsigned int, unsigned int>,
1863 std::map<unsigned int, std::vector<unsigned int>>>
1871 unsigned int n_procs = global_bboxes.size();
1872 std::vector<std::vector<unsigned int>> point_owners(n_procs);
1873 std::map<unsigned int, unsigned int> map_owners_found;
1874 std::map<unsigned int, std::vector<unsigned int>> map_owners_guessed;
1876 unsigned int n_points = points.size();
1877 for (
unsigned int pt = 0; pt < n_points; ++pt)
1880 std::vector<unsigned int> owners_found;
1882 for (
unsigned int rk = 0; rk < n_procs; ++rk)
1885 if (bbox.point_inside(points[pt]))
1887 point_owners[rk].emplace_back(pt);
1888 owners_found.emplace_back(rk);
1892 Assert(owners_found.size() > 0,
1893 ExcMessage(
"No owners found for the point " +
1894 std::to_string(pt)));
1895 if (owners_found.size() == 1)
1896 map_owners_found[pt] = owners_found[0];
1899 map_owners_guessed[pt] = owners_found;
1902 std::tuple<std::vector<std::vector<unsigned int>>,
1903 std::map<unsigned int, unsigned int>,
1904 std::map<unsigned int, std::vector<unsigned int>>>
1907 std::get<0>(output_tuple) = point_owners;
1908 std::get<1>(output_tuple) = map_owners_found;
1909 std::get<2>(output_tuple) = map_owners_guessed;
1911 return output_tuple;
1916 template <
int dim,
int spacedim>
1918 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1922 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1926 endc = triangulation.
end();
1927 for (; cell != endc; ++cell)
1928 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
1933 for (; cell != endc; ++cell)
1935 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
1937 if ((cell->at_boundary(i) ==
false) &&
1938 (cell->neighbor(i)->active()))
1941 adjacent_cell = cell->neighbor(i);
1942 for (
unsigned int j = 0;
1943 j < GeometryInfo<dim>::vertices_per_face;
1953 for (
unsigned int i = 0; i < GeometryInfo<dim>::lines_per_cell; ++i)
1954 if (cell->line(i)->has_children())
1968 template <
int dim,
int spacedim>
1969 std::map<unsigned int, types::global_vertex_index>
1973 std::map<unsigned int, types::global_vertex_index>
1974 local_to_global_vertex_index;
1976 #ifndef DEAL_II_WITH_MPI 1981 (void)triangulation;
1983 ExcMessage(
"This function does not make any sense " 1984 "for parallel::distributed::Triangulation " 1985 "objects if you do not have MPI enabled."));
1989 using active_cell_iterator =
1991 const std::vector<std::set<active_cell_iterator>> vertex_to_cell =
1996 unsigned int max_cellid_size = 0;
1997 std::set<std::pair<types::subdomain_id, types::global_vertex_index>>
1999 std::map<types::subdomain_id, std::set<unsigned int>> vertices_to_recv;
2005 active_cell_iterator cell = triangulation.
begin_active(),
2006 endc = triangulation.
end();
2007 std::set<active_cell_iterator> missing_vert_cells;
2008 std::set<unsigned int> used_vertex_index;
2009 for (; cell != endc; ++cell)
2011 if (cell->is_locally_owned())
2013 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
2016 types::subdomain_id lowest_subdomain_id = cell->subdomain_id();
2017 typename std::set<active_cell_iterator>::iterator
2018 adjacent_cell = vertex_to_cell[cell->vertex_index(i)].begin(),
2019 end_adj_cell = vertex_to_cell[cell->vertex_index(i)].end();
2020 for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
2021 lowest_subdomain_id =
2022 std::min(lowest_subdomain_id,
2023 (*adjacent_cell)->subdomain_id());
2026 if (lowest_subdomain_id == cell->subdomain_id())
2030 if (used_vertex_index.find(cell->vertex_index(i)) ==
2031 used_vertex_index.end())
2034 local_to_global_vertex_index[cell->vertex_index(i)] =
2040 vertex_to_cell[cell->vertex_index(i)].begin();
2041 for (; adjacent_cell != end_adj_cell; ++adjacent_cell)
2042 if ((*adjacent_cell)->subdomain_id() !=
2043 cell->subdomain_id())
2046 types::global_vertex_index>
2047 tmp((*adjacent_cell)->subdomain_id(),
2048 cell->vertex_index(i));
2049 if (vertices_added.find(tmp) ==
2050 vertices_added.end())
2052 vertices_to_send[(*adjacent_cell)
2055 cell->vertex_index(i),
2056 cell->id().to_string());
2057 if (cell->id().to_string().size() >
2060 cell->id().to_string().size();
2061 vertices_added.insert(tmp);
2064 used_vertex_index.insert(cell->vertex_index(i));
2071 vertices_to_recv[lowest_subdomain_id].insert(
2072 cell->vertex_index(i));
2073 missing_vert_cells.insert(cell);
2080 if (cell->is_ghost())
2082 for (
unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
2084 if (cell->at_boundary(i) ==
false)
2086 if (cell->neighbor(i)->active())
2089 spacedim>::active_cell_iterator
2090 adjacent_cell = cell->neighbor(i);
2091 if ((adjacent_cell->is_locally_owned()))
2093 types::subdomain_id adj_subdomain_id =
2094 adjacent_cell->subdomain_id();
2095 if (cell->subdomain_id() < adj_subdomain_id)
2096 for (
unsigned int j = 0;
2097 j < GeometryInfo<dim>::vertices_per_face;
2100 vertices_to_recv[cell->subdomain_id()].insert(
2101 cell->face(i)->vertex_index(j));
2102 missing_vert_cells.insert(cell);
2117 const unsigned int n_cpu =
2119 std::vector<types::global_vertex_index> indices(n_cpu);
2120 int ierr = MPI_Allgather(&next_index,
2122 DEAL_II_VERTEX_INDEX_MPI_TYPE,
2125 DEAL_II_VERTEX_INDEX_MPI_TYPE,
2131 const types::global_vertex_index
shift =
2132 std::accumulate(indices.begin(),
2136 std::map<unsigned int, types::global_vertex_index>::iterator
2137 global_index_it = local_to_global_vertex_index.begin(),
2138 global_index_end = local_to_global_vertex_index.end();
2139 for (; global_index_it != global_index_end; ++global_index_it)
2140 global_index_it->second += shift;
2147 std::vector<std::vector<types::global_vertex_index>> vertices_send_buffers(
2148 vertices_to_send.size());
2149 std::vector<MPI_Request> first_requests(vertices_to_send.size());
2153 std::string>>>::iterator
2154 vert_to_send_it = vertices_to_send.begin(),
2155 vert_to_send_end = vertices_to_send.end();
2156 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
2157 ++vert_to_send_it, ++i)
2159 int destination = vert_to_send_it->first;
2160 const unsigned int n_vertices = vert_to_send_it->second.size();
2161 const int buffer_size = 2 * n_vertices;
2162 vertices_send_buffers[i].resize(buffer_size);
2165 for (
unsigned int j = 0; j < n_vertices; ++j)
2167 vertices_send_buffers[i][2 * j] =
2168 std::get<0>(vert_to_send_it->second[j]);
2169 vertices_send_buffers[i][2 * j + 1] =
2170 local_to_global_vertex_index[std::get<1>(
2171 vert_to_send_it->second[j])];
2175 ierr = MPI_Isend(&vertices_send_buffers[i][0],
2177 DEAL_II_VERTEX_INDEX_MPI_TYPE,
2181 &first_requests[i]);
2186 std::vector<std::vector<types::global_vertex_index>> vertices_recv_buffers(
2187 vertices_to_recv.size());
2188 typename std::map<types::subdomain_id, std::set<unsigned int>>::iterator
2189 vert_to_recv_it = vertices_to_recv.begin(),
2190 vert_to_recv_end = vertices_to_recv.end();
2191 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2192 ++vert_to_recv_it, ++i)
2194 int source = vert_to_recv_it->first;
2195 const unsigned int n_vertices = vert_to_recv_it->second.size();
2196 const int buffer_size = 2 * n_vertices;
2197 vertices_recv_buffers[i].resize(buffer_size);
2200 ierr = MPI_Recv(&vertices_recv_buffers[i][0],
2202 DEAL_II_VERTEX_INDEX_MPI_TYPE,
2212 std::vector<std::vector<char>> cellids_send_buffers(
2213 vertices_to_send.size());
2214 std::vector<MPI_Request> second_requests(vertices_to_send.size());
2215 vert_to_send_it = vertices_to_send.begin();
2216 for (
unsigned int i = 0; vert_to_send_it != vert_to_send_end;
2217 ++vert_to_send_it, ++i)
2219 int destination = vert_to_send_it->first;
2220 const unsigned int n_vertices = vert_to_send_it->second.size();
2221 const int buffer_size = max_cellid_size * n_vertices;
2222 cellids_send_buffers[i].resize(buffer_size);
2225 unsigned int pos = 0;
2226 for (
unsigned int j = 0; j < n_vertices; ++j)
2228 std::string cell_id = std::get<2>(vert_to_send_it->second[j]);
2229 for (
unsigned int k = 0; k < max_cellid_size; ++k, ++pos)
2231 if (k < cell_id.size())
2232 cellids_send_buffers[i][pos] = cell_id[k];
2236 cellids_send_buffers[i][pos] =
'-';
2241 ierr = MPI_Isend(&cellids_send_buffers[i][0],
2247 &second_requests[i]);
2252 std::vector<std::vector<char>> cellids_recv_buffers(
2253 vertices_to_recv.size());
2254 vert_to_recv_it = vertices_to_recv.begin();
2255 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2256 ++vert_to_recv_it, ++i)
2258 int source = vert_to_recv_it->first;
2259 const unsigned int n_vertices = vert_to_recv_it->second.size();
2260 const int buffer_size = max_cellid_size * n_vertices;
2261 cellids_recv_buffers[i].resize(buffer_size);
2264 ierr = MPI_Recv(&cellids_recv_buffers[i][0],
2276 vert_to_recv_it = vertices_to_recv.begin();
2277 for (
unsigned int i = 0; vert_to_recv_it != vert_to_recv_end;
2278 ++i, ++vert_to_recv_it)
2280 for (
unsigned int j = 0; j < vert_to_recv_it->second.size(); ++j)
2282 const unsigned int local_pos_recv = vertices_recv_buffers[i][2 * j];
2283 const types::global_vertex_index global_id_recv =
2284 vertices_recv_buffers[i][2 * j + 1];
2285 const std::string cellid_recv(
2286 &cellids_recv_buffers[i][max_cellid_size * j],
2287 &cellids_recv_buffers[i][max_cellid_size * j] + max_cellid_size);
2289 typename std::set<active_cell_iterator>::iterator
2290 cell_set_it = missing_vert_cells.begin(),
2291 end_cell_set = missing_vert_cells.end();
2292 for (; (found ==
false) && (cell_set_it != end_cell_set);
2295 typename std::set<active_cell_iterator>::iterator
2297 vertex_to_cell[(*cell_set_it)->vertex_index(i)].begin(),
2299 vertex_to_cell[(*cell_set_it)->vertex_index(i)].end();
2300 for (; candidate_cell != end_cell; ++candidate_cell)
2302 std::string current_cellid =
2303 (*candidate_cell)->id().to_string();
2304 current_cellid.resize(max_cellid_size,
'-');
2305 if (current_cellid.compare(cellid_recv) == 0)
2307 local_to_global_vertex_index
2308 [(*candidate_cell)->vertex_index(local_pos_recv)] =
2320 return local_to_global_vertex_index;
2325 template <
int dim,
int spacedim>
2339 std::map<std::pair<unsigned int, unsigned int>,
unsigned int> indexmap;
2340 for (typename ::internal::
2343 cell != triangulation.
end();
2345 indexmap[std::pair<unsigned int, unsigned int>(cell->level(),
2347 cell->active_cell_index();
2356 for (typename ::internal::
2359 cell != triangulation.
end();
2362 const unsigned int index = cell->active_cell_index();
2363 cell_connectivity.
add(index, index);
2364 for (
unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
2365 if ((cell->at_boundary(f) ==
false) &&
2366 (cell->neighbor(f)->has_children() ==
false))
2368 unsigned int other_index =
2370 .find(std::pair<unsigned int, unsigned int>(
2371 cell->neighbor(f)->level(), cell->neighbor(f)->index()))
2373 cell_connectivity.
add(index, other_index);
2374 cell_connectivity.
add(other_index, index);
2381 template <
int dim,
int spacedim>
2387 std::vector<std::vector<unsigned int>> vertex_to_cell(
2391 cell != triangulation.
end();
2394 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2395 vertex_to_cell[cell->vertex_index(v)].push_back(
2396 cell->active_cell_index());
2403 cell != triangulation.
end();
2406 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2407 for (
unsigned int n = 0;
2408 n < vertex_to_cell[cell->vertex_index(v)].size();
2410 cell_connectivity.
add(cell->active_cell_index(),
2411 vertex_to_cell[cell->vertex_index(v)][n]);
2416 template <
int dim,
int spacedim>
2420 const unsigned int level,
2423 std::vector<std::vector<unsigned int>> vertex_to_cell(
2426 triangulation.
begin(level);
2427 cell != triangulation.
end(level);
2430 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2431 vertex_to_cell[cell->vertex_index(v)].push_back(cell->index());
2435 triangulation.
n_cells(level));
2437 triangulation.
begin(level);
2438 cell != triangulation.
end(level);
2441 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2442 for (
unsigned int n = 0;
2443 n < vertex_to_cell[cell->vertex_index(v)].size();
2445 cell_connectivity.
add(cell->index(),
2446 vertex_to_cell[cell->vertex_index(v)][n]);
2452 template <
int dim,
int spacedim>
2458 std::vector<unsigned int> cell_weights;
2461 if (!triangulation.
signals.cell_weight.empty())
2464 std::numeric_limits<unsigned int>::max());
2469 endc = triangulation.
end();
2470 for (; cell != endc; ++cell, ++c)
2471 cell_weights[c] = triangulation.
signals.cell_weight(
2484 template <
int dim,
int spacedim>
2487 const std::vector<unsigned int> &cell_weights,
2492 &triangulation) ==
nullptr),
2493 ExcMessage(
"Objects of type parallel::distributed::Triangulation " 2494 "are already partitioned implicitly and can not be " 2495 "partitioned again explicitly."));
2499 if (n_partitions == 1)
2501 for (typename ::internal::ActiveCellIterator<
2506 cell != triangulation.
end();
2508 cell->set_subdomain_id(0);
2522 sp_cell_connectivity.
copy_from(cell_connectivity);
2525 sp_cell_connectivity,
2532 template <
int dim,
int spacedim>
2539 std::vector<unsigned int> cell_weights;
2542 if (!triangulation.
signals.cell_weight.empty())
2545 std::numeric_limits<unsigned int>::max());
2550 endc = triangulation.
end();
2551 for (; cell != endc; ++cell, ++c)
2552 cell_weights[c] = triangulation.
signals.cell_weight(
2559 cell_connection_graph,
2566 template <
int dim,
int spacedim>
2569 const std::vector<unsigned int> &cell_weights,
2575 &triangulation) ==
nullptr),
2576 ExcMessage(
"Objects of type parallel::distributed::Triangulation " 2577 "are already partitioned implicitly and can not be " 2578 "partitioned again explicitly."));
2581 ExcMessage(
"Connectivity graph has wrong size"));
2583 ExcMessage(
"Connectivity graph has wrong size"));
2586 if (n_partitions == 1)
2588 for (typename ::internal::ActiveCellIterator<
2593 cell != triangulation.
end();
2595 cell->set_subdomain_id(0);
2603 std::vector<unsigned int> partition_indices(triangulation.
n_active_cells());
2612 for (typename ::internal::
2615 cell != triangulation.
end();
2617 cell->set_subdomain_id(partition_indices[cell->active_cell_index()]);
2628 set_subdomain_id_in_zorder_recursively(IT cell,
2629 unsigned int & current_proc_idx,
2630 unsigned int & current_cell_idx,
2631 const unsigned int n_active_cells,
2632 const unsigned int n_partitions)
2638 floor((
long)n_active_cells * (current_proc_idx + 1) / n_partitions))
2640 cell->set_subdomain_id(current_proc_idx);
2645 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2646 set_subdomain_id_in_zorder_recursively(cell->child(n),
2655 template <
int dim,
int spacedim>
2661 &triangulation) ==
nullptr),
2662 ExcMessage(
"Objects of type parallel::distributed::Triangulation " 2663 "are already partitioned implicitly and can not be " 2664 "partitioned again explicitly."));
2668 if (n_partitions == 1)
2670 for (typename ::internal::ActiveCellIterator<
2675 cell != triangulation.
end();
2677 cell->set_subdomain_id(0);
2683 std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
2684 std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
2690 coarse_cell_to_p4est_tree_permutation.resize(triangulation.
n_cells(0));
2692 coarse_cell_to_p4est_tree_permutation);
2694 p4est_tree_to_coarse_cell_permutation =
2697 unsigned int current_proc_idx = 0;
2698 unsigned int current_cell_idx = 0;
2699 const unsigned int n_active_cells = triangulation.
n_active_cells();
2703 for (
unsigned int idx = 0; idx < triangulation.
n_cells(0); ++idx)
2705 const unsigned int coarse_cell_idx =
2706 p4est_tree_to_coarse_cell_permutation[idx];
2708 &triangulation, 0, coarse_cell_idx);
2710 set_subdomain_id_in_zorder_recursively(coarse_cell,
2726 cell = triangulation.
begin(),
2727 endc = triangulation.
end();
2728 for (; cell != endc; ++cell)
2732 bool all_children_active =
true;
2733 std::map<unsigned int, unsigned int> map_cpu_n_cells;
2734 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2735 if (!cell->child(n)->active())
2737 all_children_active =
false;
2741 ++map_cpu_n_cells[cell->child(n)->subdomain_id()];
2743 if (!all_children_active)
2746 unsigned int new_owner = cell->child(0)->subdomain_id();
2747 for (std::map<unsigned int, unsigned int>::iterator it =
2748 map_cpu_n_cells.begin();
2749 it != map_cpu_n_cells.end();
2751 if (it->second > map_cpu_n_cells[new_owner])
2752 new_owner = it->first;
2754 for (
unsigned int n = 0; n < cell->n_children(); ++n)
2755 cell->child(n)->set_subdomain_id(new_owner);
2761 template <
int dim,
int spacedim>
2765 unsigned int n_levels = triangulation.
n_levels();
2766 for (
int lvl = n_levels - 1; lvl >= 0; --lvl)
2769 cell = triangulation.
begin(lvl),
2770 endc = triangulation.
end(lvl);
2771 for (; cell != endc; ++cell)
2773 if (!cell->has_children())
2774 cell->set_level_subdomain_id(cell->subdomain_id());
2777 Assert(cell->child(0)->level_subdomain_id() !=
2780 cell->set_level_subdomain_id(
2781 cell->child(0)->level_subdomain_id());
2788 template <
int dim,
int spacedim>
2791 std::vector<types::subdomain_id> & subdomain)
2798 cell != triangulation.
end();
2800 subdomain[cell->active_cell_index()] = cell->subdomain_id();
2805 template <
int dim,
int spacedim>
2811 unsigned int count = 0;
2814 cell != triangulation.
end();
2816 if (cell->subdomain_id() == subdomain)
2824 template <
int dim,
int spacedim>
2829 std::vector<bool> locally_owned_vertices =
2839 for (typename ::internal::ActiveCellIterator<
2844 cell != triangulation.
end();
2846 if (cell->is_artificial() ||
2847 (cell->is_ghost() &&
2848 (cell->subdomain_id() < tr->locally_owned_subdomain())))
2849 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
2851 locally_owned_vertices[cell->vertex_index(v)] =
false;
2853 return locally_owned_vertices;
2860 template <
int dim,
int spacedim>
2869 return (vertices[1] - vertices[0]).norm();
2871 return std::max((vertices[3] - vertices[0]).norm(),
2872 (vertices[2] - vertices[1]).norm());
2874 return std::max(std::max((vertices[7] - vertices[0]).norm(),
2875 (vertices[6] - vertices[1]).norm()),
2876 std::max((vertices[2] - vertices[5]).norm(),
2877 (vertices[3] - vertices[4]).norm()));
2886 template <
int dim,
int spacedim>
2891 double min_diameter = std::numeric_limits<double>::max();
2893 if (!cell->is_artificial())
2895 std::min(min_diameter, diameter<dim, spacedim>(cell, mapping));
2897 double global_min_diameter = 0;
2899 #ifdef DEAL_II_WITH_MPI 2903 global_min_diameter =
2907 global_min_diameter = min_diameter;
2909 return global_min_diameter;
2914 template <
int dim,
int spacedim>
2919 double max_diameter = 0.;
2921 if (!cell->is_artificial())
2922 max_diameter = std::max(max_diameter,
diameter(cell, mapping));
2924 double global_max_diameter = 0;
2926 #ifdef DEAL_II_WITH_MPI 2930 global_max_diameter =
2934 global_max_diameter = max_diameter;
2936 return global_max_diameter;
2943 namespace FixUpDistortedChildCells
2966 template <
typename Iterator,
int spacedim>
2968 objective_function(
const Iterator &
object,
2971 const unsigned int structdim =
2972 Iterator::AccessorType::structure_dimension;
2973 Assert(spacedim == Iterator::AccessorType::dimension,
2979 Assert(object->refinement_case() ==
2987 Tensor<spacedim - structdim, spacedim>
2990 for (
unsigned int i = 0; i < GeometryInfo<structdim>::vertices_per_cell;
2992 parent_vertices[i] = object->vertex(i);
2995 parent_vertices, parent_alternating_forms);
2997 const Tensor<spacedim - structdim, spacedim>
2998 average_parent_alternating_form =
2999 std::accumulate(parent_alternating_forms,
3000 parent_alternating_forms +
3014 Tensor<spacedim - structdim, spacedim> child_alternating_forms
3018 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3019 for (
unsigned int i = 0;
3020 i < GeometryInfo<structdim>::vertices_per_cell;
3022 child_vertices[c][i] = object->child(c)->vertex(i);
3030 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3032 1] = object_mid_point;
3034 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3036 child_vertices[c], child_alternating_forms[c]);
3048 double objective = 0;
3049 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3050 for (
unsigned int i = 0;
3051 i < GeometryInfo<structdim>::vertices_per_cell;
3054 (child_alternating_forms[c][i] -
3055 average_parent_alternating_form / std::pow(2., 1. * structdim))
3067 template <
typename Iterator>
3069 get_face_midpoint(
const Iterator &
object,
3070 const unsigned int f,
3071 std::integral_constant<int, 1>)
3073 return object->vertex(f);
3083 template <
typename Iterator>
3085 get_face_midpoint(
const Iterator &
object,
3086 const unsigned int f,
3087 std::integral_constant<int, 2>)
3089 return object->line(f)->center();
3099 template <
typename Iterator>
3101 get_face_midpoint(
const Iterator &
object,
3102 const unsigned int f,
3103 std::integral_constant<int, 3>)
3105 return object->face(f)->center();
3132 template <
typename Iterator>
3134 minimal_diameter(
const Iterator &
object)
3136 const unsigned int structdim =
3137 Iterator::AccessorType::structure_dimension;
3139 double diameter =
object->diameter();
3140 for (
unsigned int f = 0; f < GeometryInfo<structdim>::faces_per_cell;
3142 for (
unsigned int e = f + 1;
3143 e < GeometryInfo<structdim>::faces_per_cell;
3145 diameter = std::min(
3147 get_face_midpoint(
object,
3149 std::integral_constant<int, structdim>())
3150 .distance(get_face_midpoint(
3151 object, e, std::integral_constant<int, structdim>())));
3162 template <
typename Iterator>
3164 fix_up_object(
const Iterator &
object)
3166 const unsigned int structdim =
3167 Iterator::AccessorType::structure_dimension;
3168 const unsigned int spacedim = Iterator::AccessorType::space_dimension;
3175 Assert(object->refinement_case() ==
3185 unsigned int iteration = 0;
3186 const double diameter = minimal_diameter(
object);
3189 double current_value = objective_function(
object, object_mid_point);
3190 double initial_delta = 0;
3197 const double step_length = diameter / 4 / (iteration + 1);
3202 for (
unsigned int d = 0; d < spacedim; ++d)
3204 const double eps = step_length / 10;
3210 (objective_function(
3218 if (gradient.
norm() == 0)
3227 std::min(2 * current_value / (gradient * gradient),
3228 step_length / gradient.norm()) *
3233 const double previous_value = current_value;
3234 current_value = objective_function(
object, object_mid_point);
3237 initial_delta = (previous_value - current_value);
3240 if ((iteration >= 1) &&
3241 ((previous_value - current_value < 0) ||
3242 (std::fabs(previous_value - current_value) <
3243 0.001 * initial_delta)))
3248 while (iteration < 20);
3264 double old_min_product, new_min_product;
3268 for (
unsigned int i = 0; i < GeometryInfo<structdim>::vertices_per_cell;
3270 parent_vertices[i] = object->vertex(i);
3272 Tensor<spacedim - structdim, spacedim>
3275 parent_vertices, parent_alternating_forms);
3281 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3282 for (
unsigned int i = 0;
3283 i < GeometryInfo<structdim>::vertices_per_cell;
3285 child_vertices[c][i] = object->child(c)->vertex(i);
3287 Tensor<spacedim - structdim, spacedim> child_alternating_forms
3291 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3293 child_vertices[c], child_alternating_forms[c]);
3296 child_alternating_forms[0][0] * parent_alternating_forms[0];
3297 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3298 for (
unsigned int i = 0;
3299 i < GeometryInfo<structdim>::vertices_per_cell;
3301 for (
unsigned int j = 0;
3302 j < GeometryInfo<structdim>::vertices_per_cell;
3304 old_min_product = std::min<double>(old_min_product,
3305 child_alternating_forms[c][i] *
3306 parent_alternating_forms[j]);
3314 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3316 1] = object_mid_point;
3318 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3320 child_vertices[c], child_alternating_forms[c]);
3323 child_alternating_forms[0][0] * parent_alternating_forms[0];
3324 for (
unsigned int c = 0; c <
object->n_children(); ++c)
3325 for (
unsigned int i = 0;
3326 i < GeometryInfo<structdim>::vertices_per_cell;
3328 for (
unsigned int j = 0;
3329 j < GeometryInfo<structdim>::vertices_per_cell;
3331 new_min_product = std::min<double>(new_min_product,
3332 child_alternating_forms[c][i] *
3333 parent_alternating_forms[j]);
3341 if (new_min_product >= old_min_product)
3342 object->child(0)->vertex(
3349 return (std::max(new_min_product, old_min_product) > 0);
3355 fix_up_faces(const ::Triangulation<1, 1>::cell_iterator &,
3356 std::integral_constant<int, 1>,
3357 std::integral_constant<int, 1>)
3365 template <
int dim,
int spacedim>
3368 const typename ::Triangulation<dim, spacedim>::cell_iterator
3370 std::integral_constant<int, dim>,
3371 std::integral_constant<int, spacedim>)
3381 for (
unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3384 Assert(cell->face(f)->refinement_case() ==
3388 bool subface_is_more_refined =
false;
3389 for (
unsigned int g = 0;
3390 g < GeometryInfo<dim>::max_children_per_face;
3392 if (cell->face(f)->child(g)->has_children())
3394 subface_is_more_refined =
true;
3398 if (subface_is_more_refined ==
true)
3402 fix_up_object(cell->face(f));
3409 template <
int dim,
int spacedim>
3419 for (
typename std::list<
3430 "This function is only valid for a list of cells that " 3431 "have children (i.e., no cell in the list may be active)."));
3433 internal::FixUpDistortedChildCells ::fix_up_faces(
3435 std::integral_constant<int, dim>(),
3436 std::integral_constant<int, spacedim>());
3439 if (!internal::FixUpDistortedChildCells::fix_up_object(cell))
3443 return unfixable_subset;
3448 template <
int dim,
int spacedim>
3451 const bool reset_boundary_ids)
3454 std::vector<types::manifold_id> dst_manifold_ids(src_boundary_ids.size());
3455 auto m_it = dst_manifold_ids.begin();
3456 for (
auto b : src_boundary_ids)
3461 const std::vector<types::boundary_id> reset_boundary_id =
3462 reset_boundary_ids ?
3463 std::vector<types::boundary_id>(src_boundary_ids.size(), 0) :
3473 template <
int dim,
int spacedim>
3476 const std::vector<types::boundary_id> &src_boundary_ids,
3477 const std::vector<types::manifold_id> &dst_manifold_ids,
3479 const std::vector<types::boundary_id> &reset_boundary_ids_)
3482 const auto reset_boundary_ids =
3483 reset_boundary_ids_.size() ? reset_boundary_ids_ : src_boundary_ids;
3495 for (
unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3496 if (cell->face(f)->at_boundary())
3497 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_face; ++e)
3499 const auto bid = cell->face(f)->line(e)->boundary_id();
3500 const auto ind = std::find(src_boundary_ids.begin(),
3501 src_boundary_ids.end(),
3503 src_boundary_ids.begin();
3504 if ((
unsigned int)ind < src_boundary_ids.size())
3505 cell->face(f)->line(e)->set_manifold_id(
3506 dst_manifold_ids[ind]);
3514 for (
unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3515 if (cell->face(f)->at_boundary())
3517 const auto bid = cell->face(f)->boundary_id();
3519 std::find(src_boundary_ids.begin(), src_boundary_ids.end(), bid) -
3520 src_boundary_ids.begin();
3522 if ((
unsigned int)ind < src_boundary_ids.size())
3525 cell->face(f)->set_manifold_id(dst_manifold_ids[ind]);
3527 cell->face(f)->set_boundary_id(reset_boundary_ids[ind]);
3531 for (
unsigned int e = 0; e < GeometryInfo<dim>::lines_per_face;
3534 const auto bid = cell->face(f)->line(e)->boundary_id();
3535 const auto ind = std::find(src_boundary_ids.begin(),
3536 src_boundary_ids.end(),
3538 src_boundary_ids.begin();
3539 if ((
unsigned int)ind < src_boundary_ids.size())
3540 cell->face(f)->line(e)->set_boundary_id(
3541 reset_boundary_ids[ind]);
3547 template <
int dim,
int spacedim>
3550 const bool compute_face_ids)
3556 for (; cell != endc; ++cell)
3558 cell->set_manifold_id(cell->material_id());
3559 if (compute_face_ids ==
true)
3561 for (
unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3563 if (cell->at_boundary(f) ==
false)
3564 cell->face(f)->set_manifold_id(
3565 std::min(cell->material_id(),
3566 cell->neighbor(f)->material_id()));
3568 cell->face(f)->set_manifold_id(cell->material_id());
3574 template <
int dim,
int spacedim>
3575 std::pair<unsigned int, double>
3579 double max_ratio = 1;
3580 unsigned int index = 0;
3582 for (
unsigned int i = 0; i < dim; ++i)
3583 for (
unsigned int j = i + 1; j < dim; ++j)
3585 unsigned int ax = i % dim;
3586 unsigned int next_ax = j % dim;
3589 cell->extent_in_direction(ax) / cell->extent_in_direction(next_ax);
3591 if (ratio > max_ratio)
3596 else if (1.0 / ratio > max_ratio)
3598 max_ratio = 1.0 / ratio;
3602 return std::make_pair(index, max_ratio);
3606 template <
int dim,
int spacedim>
3609 const bool isotropic,
3610 const unsigned int max_iterations)
3612 unsigned int iter = 0;
3613 bool continue_refinement =
true;
3619 while (continue_refinement && (iter < max_iterations))
3623 continue_refinement =
false;
3626 for (
unsigned int j = 0; j < GeometryInfo<dim>::faces_per_cell; j++)
3627 if (cell->at_boundary(j) ==
false &&
3628 cell->neighbor(j)->has_children())
3632 cell->set_refine_flag();
3633 continue_refinement =
true;
3636 continue_refinement |= cell->flag_for_face_refinement(j);
3643 template <
int dim,
int spacedim>
3646 const double max_ratio,
3647 const unsigned int max_iterations)
3649 unsigned int iter = 0;
3650 bool continue_refinement =
true;
3656 while (continue_refinement && (iter < max_iterations))
3659 continue_refinement =
false;
3662 std::pair<unsigned int, double> info =
3663 GridTools::get_longest_direction<dim, spacedim>(cell);
3664 if (info.second > max_ratio)
3666 cell->set_refine_flag(
3668 continue_refinement =
true;
3676 template <
int dim,
int spacedim>
3679 const double limit_angle_fraction)
3687 "have hanging nodes."));
3690 bool has_cells_with_more_than_dim_faces_on_boundary =
true;
3691 bool has_cells_with_dim_faces_on_boundary =
false;
3693 unsigned int refinement_cycles = 0;
3695 while (has_cells_with_more_than_dim_faces_on_boundary)
3697 has_cells_with_more_than_dim_faces_on_boundary =
false;
3701 unsigned int boundary_face_counter = 0;
3702 for (
unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3703 if (cell->face(f)->at_boundary())
3704 boundary_face_counter++;
3705 if (boundary_face_counter > dim)
3707 has_cells_with_more_than_dim_faces_on_boundary =
true;
3710 else if (boundary_face_counter == dim)
3711 has_cells_with_dim_faces_on_boundary =
true;
3713 if (has_cells_with_more_than_dim_faces_on_boundary)
3716 refinement_cycles++;
3720 if (has_cells_with_dim_faces_on_boundary)
3723 refinement_cycles++;
3727 while (refinement_cycles > 0)
3730 cell->set_coarsen_flag();
3732 refinement_cycles--;
3738 std::vector<Point<spacedim>> vertices = tria.
get_vertices();
3740 std::vector<bool> faces_to_remove(tria.
n_raw_faces(),
false);
3742 std::vector<CellData<dim>> cells_to_add;
3746 const unsigned int v0 = 0, v1 = 1, v2 = (dim > 1 ? 2 : 0),
3747 v3 = (dim > 1 ? 3 : 0);
3751 double angle_fraction = 0;
3757 p0[spacedim > 1 ? 1 : 0] = 1;
3761 if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
3763 p0 = cell->vertex(v0) - cell->vertex(v2);
3764 p1 = cell->vertex(v3) - cell->vertex(v2);
3765 vertex_at_corner = v2;
3767 else if (cell->face(v3)->at_boundary() &&
3768 cell->face(v1)->at_boundary())
3770 p0 = cell->vertex(v2) - cell->vertex(v3);
3771 p1 = cell->vertex(v1) - cell->vertex(v3);
3772 vertex_at_corner = v3;
3774 else if (cell->face(1)->at_boundary() &&
3775 cell->face(2)->at_boundary())
3777 p0 = cell->vertex(v0) - cell->vertex(v1);
3778 p1 = cell->vertex(v3) - cell->vertex(v1);
3779 vertex_at_corner = v1;
3781 else if (cell->face(2)->at_boundary() &&
3782 cell->face(0)->at_boundary())
3784 p0 = cell->vertex(v2) - cell->vertex(v0);
3785 p1 = cell->vertex(v1) - cell->vertex(v0);
3786 vertex_at_corner = v0;
3790 angle_fraction = std::acos(p0 * p1) /
numbers::PI;
3797 if (angle_fraction > limit_angle_fraction)
3799 auto flags_removal = [&](
unsigned int f1,
3802 unsigned int n2) ->
void {
3803 cells_to_remove[cell->active_cell_index()] =
true;
3804 cells_to_remove[cell->neighbor(n1)->active_cell_index()] =
true;
3805 cells_to_remove[cell->neighbor(n2)->active_cell_index()] =
true;
3807 faces_to_remove[cell->face(f1)->index()] =
true;
3808 faces_to_remove[cell->face(f2)->index()] =
true;
3810 faces_to_remove[cell->neighbor(n1)->face(f1)->index()] =
true;
3811 faces_to_remove[cell->neighbor(n2)->face(f2)->index()] =
true;
3814 auto cell_creation = [&](
const unsigned int vv0,
3815 const unsigned int vv1,
3816 const unsigned int f0,
3817 const unsigned int f1,
3819 const unsigned int n0,
3820 const unsigned int v0n0,
3821 const unsigned int v1n0,
3823 const unsigned int n1,
3824 const unsigned int v0n1,
3825 const unsigned int v1n1) {
3829 c1.
vertices[v0] = cell->vertex_index(vv0);
3830 c1.
vertices[v1] = cell->vertex_index(vv1);
3831 c1.
vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
3832 c1.
vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
3837 c2.
vertices[v0] = cell->vertex_index(vv0);
3838 c2.
vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
3839 c2.
vertices[v2] = cell->vertex_index(vv1);
3840 c2.
vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
3845 l1.
vertices[0] = cell->vertex_index(vv0);
3846 l1.
vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
3852 l2.
vertices[0] = cell->vertex_index(vv0);
3853 l2.
vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
3859 cells_to_add.push_back(c1);
3860 cells_to_add.push_back(c2);
3865 switch (vertex_at_corner)
3868 flags_removal(0, 2, 3, 1);
3869 cell_creation(0, 3, 0, 2, 3, 2, 3, 1, 1, 3);
3872 flags_removal(1, 2, 3, 0);
3873 cell_creation(1, 2, 2, 1, 0, 0, 2, 3, 3, 2);
3876 flags_removal(3, 0, 1, 2);
3877 cell_creation(2, 1, 3, 0, 1, 3, 1, 2, 0, 1);
3880 flags_removal(3, 1, 0, 2);
3881 cell_creation(3, 0, 1, 3, 2, 1, 0, 0, 2, 0);
3894 if (cells_to_add.size() == 0)
3896 while (refinement_cycles > 0)
3899 cell->set_coarsen_flag();
3901 refinement_cycles--;
3909 if (cells_to_remove[cell->active_cell_index()] ==
false)
3912 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
3914 c.
vertices[v] = cell->vertex_index(v);
3917 cells_to_add.push_back(c);
3925 for (; face != endf; ++face)
3926 if ((face->at_boundary() ||
3928 faces_to_remove[face->index()] ==
false)
3930 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_face; ++l)
3935 for (
unsigned int v = 0;
3936 v < GeometryInfo<1>::vertices_per_cell;
3938 line.
vertices[v] = face->vertex_index(v);
3944 for (
unsigned int v = 0;
3945 v < GeometryInfo<1>::vertices_per_cell;
3947 line.
vertices[v] = face->line(l)->vertex_index(v);
3956 for (
unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell;
3958 quad.
vertices[v] = face->vertex_index(v);
3966 subcelldata_to_add);
3971 std::map<types::manifold_id, std::unique_ptr<Manifold<dim, spacedim>>>
3974 for (
auto manifold_id : manifold_ids)
3976 manifolds[manifold_id] = tria.
get_manifold(manifold_id).clone();
3983 for (
auto manifold_id : manifold_ids)
3985 tria.
set_manifold(manifold_id, *manifolds[manifold_id]);
3990 template <
int dim,
int spacedim>
3993 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
3994 std::vector<std::vector<Point<dim>>>,
3995 std::vector<std::vector<unsigned int>>>
4006 const unsigned int np = points.size();
4009 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
4010 std::vector<std::vector<Point<dim>>>,
4011 std::vector<std::vector<unsigned int>>>
4016 return cell_qpoint_map;
4019 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4028 std::get<0>(cell_qpoint_map).emplace_back(my_pair.first);
4029 std::get<1>(cell_qpoint_map).emplace_back(1, my_pair.second);
4030 std::get<2>(cell_qpoint_map).emplace_back(1, 0);
4034 return cell_qpoint_map;
4036 Point<spacedim> cell_center = std::get<0>(cell_qpoint_map)[0]->center();
4037 double cell_diameter = std::get<0>(cell_qpoint_map)[0]->
diameter() *
4038 (0.5 + std::numeric_limits<double>::epsilon());
4041 for (
unsigned int p = 1; p < np; ++p)
4045 if (cell_center.
distance(points[p]) < cell_diameter)
4047 cache, points[p], std::get<0>(cell_qpoint_map).back());
4052 if (my_pair.first == std::get<0>(cell_qpoint_map).back())
4055 std::get<1>(cell_qpoint_map).back().emplace_back(my_pair.second);
4056 std::get<2>(cell_qpoint_map).back().emplace_back(p);
4061 typename std::vector<typename Triangulation<dim, spacedim>::
4062 active_cell_iterator>::iterator cells_it =
4063 std::find(std::get<0>(cell_qpoint_map).begin(),
4064 std::get<0>(cell_qpoint_map).end() - 1,
4067 if (cells_it == std::get<0>(cell_qpoint_map).end() - 1)
4070 std::get<0>(cell_qpoint_map).emplace_back(my_pair.first);
4071 std::get<1>(cell_qpoint_map).emplace_back(1, my_pair.second);
4072 std::get<2>(cell_qpoint_map).emplace_back(1, p);
4074 cell_center = std::get<0>(cell_qpoint_map).back()->center();
4076 std::get<0>(cell_qpoint_map).back()->diameter() *
4077 (0.5 + std::numeric_limits<double>::epsilon());
4081 unsigned int current_cell =
4082 cells_it - std::get<0>(cell_qpoint_map).begin();
4085 std::get<1>(cell_qpoint_map)[current_cell].emplace_back(
4087 std::get<2>(cell_qpoint_map)[current_cell].emplace_back(p);
4093 Assert(std::get<0>(cell_qpoint_map).size() ==
4094 std::get<2>(cell_qpoint_map).size(),
4096 std::get<2>(cell_qpoint_map).size()));
4098 Assert(std::get<0>(cell_qpoint_map).size() ==
4099 std::get<1>(cell_qpoint_map).size(),
4101 std::get<1>(cell_qpoint_map).size()));
4104 unsigned int c = std::get<0>(cell_qpoint_map).size();
4105 unsigned int qps = 0;
4110 for (
unsigned int n = 0; n < c; ++n)
4112 Assert(std::get<1>(cell_qpoint_map)[n].size() ==
4113 std::get<2>(cell_qpoint_map)[n].size(),
4115 std::get<2>(cell_qpoint_map)[n].size()));
4116 qps += std::get<1>(cell_qpoint_map)[n].size();
4121 return cell_qpoint_map;
4129 namespace distributed_cptloc
4132 template <
int dim,
int spacedim>
4141 return k->active_cell_index();
4149 template <
int dim,
int spacedim>
4152 std::pair<std::vector<Point<dim>>, std::vector<unsigned int>>,
4153 cell_hash<dim, spacedim>>
4154 compute_point_locations_unmap(
4159 const unsigned int np = points.size();
4162 typename Triangulation<dim, spacedim>::active_cell_iterator,
4163 std::pair<std::vector<Point<dim>>, std::vector<unsigned int>>,
4164 cell_hash<dim, spacedim>>
4169 return cell_qpoint_map;
4174 auto last_cell = cell_qpoint_map.emplace(
4175 std::make_pair(my_pair.first,
4176 std::make_pair(std::vector<
Point<dim>>{my_pair.second},
4177 std::vector<unsigned int>{0})));
4180 return cell_qpoint_map;
4183 double cell_diameter = my_pair.first->diameter() *
4184 (0.5 + std::numeric_limits<double>::epsilon());
4187 for (
unsigned int p = 1; p < np; ++p)
4191 if (cell_center.
distance(points[p]) < cell_diameter)
4193 cache, points[p], last_cell.first->first);
4198 if (last_cell.first->first == my_pair.first)
4200 last_cell.first->second.first.emplace_back(my_pair.second);
4201 last_cell.first->second.second.emplace_back(p);
4206 last_cell = cell_qpoint_map.emplace(std::make_pair(
4208 std::make_pair(std::vector<
Point<dim>>{my_pair.second},
4209 std::vector<unsigned int>{p})));
4211 if (last_cell.second ==
false)
4214 last_cell.first->second.first.emplace_back(my_pair.second);
4215 last_cell.first->second.second.emplace_back(p);
4220 cell_center = my_pair.first->center();
4222 my_pair.first->diameter() *
4223 (0.5 + std::numeric_limits<double>::epsilon());
4229 unsigned int qps = 0;
4234 for (
const auto &m : cell_qpoint_map)
4236 Assert(m.second.second.size() == m.second.first.size(),
4238 m.second.first.size()));
4239 qps += m.second.second.size();
4243 return cell_qpoint_map;
4256 template <
int dim,
int spacedim>
4258 merge_cptloc_outputs(
4260 typename Triangulation<dim, spacedim>::active_cell_iterator,
4262 std::vector<unsigned int>,
4264 std::vector<unsigned int>>,
4265 cell_hash<dim, spacedim>> &output_unmap,
4267 typename Triangulation<dim, spacedim>::active_cell_iterator>
4269 const std::vector<std::vector<
Point<dim>>> & in_qpoints,
4270 const std::vector<std::vector<unsigned int>> & in_maps,
4272 const unsigned int in_rank)
4275 for (
unsigned int c = 0; c < in_cells.size(); ++c)
4278 auto current_c = output_unmap.emplace(
4279 std::make_pair(in_cells[c],
4280 std::make_tuple(in_qpoints[c],
4283 std::vector<unsigned int>(
4284 in_points[c].size(), in_rank))));
4286 if (current_c.second ==
false)
4290 auto &cell_qpts = std::get<0>(current_c.first->second);
4291 auto &cell_maps = std::get<1>(current_c.first->second);
4292 auto &cell_pts = std::get<2>(current_c.first->second);
4293 auto &cell_ranks = std::get<3>(current_c.first->second);
4294 cell_qpts.insert(cell_qpts.end(),
4295 in_qpoints[c].begin(),
4296 in_qpoints[c].end());
4297 cell_maps.insert(cell_maps.end(),
4300 cell_pts.insert(cell_pts.end(),
4301 in_points[c].begin(),
4302 in_points[c].end());
4303 std::vector<unsigned int> ranks_tmp(in_points[c].size(),
4305 cell_ranks.insert(cell_ranks.end(),
4323 template <
int dim,
int spacedim>
4325 compute_and_classify_points(
4328 const std::vector<unsigned int> & local_points_idx,
4330 typename Triangulation<dim, spacedim>::active_cell_iterator,
4332 std::vector<unsigned int>,
4334 std::vector<unsigned int>>,
4335 cell_hash<dim, spacedim>> &output_unmap,
4336 std::map<
unsigned int,
4337 std::tuple<std::vector<CellId>,
4339 std::vector<std::vector<unsigned int>>,
4342 std::vector<unsigned int> &classified_pts)
4344 auto cpt_loc_pts = compute_point_locations_unmap(cache, local_points);
4348 for (
auto const &cell_tuples : cpt_loc_pts)
4350 auto &cell_loc = cell_tuples.first;
4351 auto &q_loc = std::get<0>(cell_tuples.second);
4352 auto &indices_loc = std::get<1>(cell_tuples.second);
4353 if (cell_loc->is_locally_owned())
4356 std::vector<Point<spacedim>> cell_points(indices_loc.size());
4357 std::vector<unsigned int> cell_points_idx(indices_loc.size());
4358 for (
unsigned int i = 0; i < indices_loc.size(); ++i)
4361 cell_points[i] = local_points[indices_loc[i]];
4366 cell_points_idx[i] = local_points_idx[indices_loc[i]];
4367 classified_pts.emplace_back(
4368 local_points_idx[indices_loc[i]]);
4370 output_unmap.emplace(
4371 std::make_pair(cell_loc,
4372 std::make_tuple(q_loc,
4375 std::vector<unsigned int>(
4377 cell_loc->subdomain_id()))));
4379 else if (cell_loc->is_ghost())
4383 std::vector<Point<spacedim>> cell_points(indices_loc.size());
4384 std::vector<unsigned int> cell_points_idx(indices_loc.size());
4385 for (
unsigned int i = 0; i < indices_loc.size(); ++i)
4387 cell_points[i] = local_points[indices_loc[i]];
4388 cell_points_idx[i] = local_points_idx[indices_loc[i]];
4389 classified_pts.emplace_back(
4390 local_points_idx[indices_loc[i]]);
4396 auto &map_tuple_owner = ghost_loc_pts[cell_loc->subdomain_id()];
4398 std::get<0>(map_tuple_owner).emplace_back(cell_loc->id());
4399 std::get<1>(map_tuple_owner).emplace_back(q_loc);
4400 std::get<2>(map_tuple_owner).emplace_back(cell_points_idx);
4401 std::get<3>(map_tuple_owner).emplace_back(cell_points);
4414 template <
int dim,
int spacedim>
4416 compute_and_merge_from_map(
4418 const std::map<
unsigned int,
4420 std::vector<unsigned int>>> &map_pts,
4422 typename Triangulation<dim, spacedim>::active_cell_iterator,
4424 std::vector<unsigned int>,
4426 std::vector<unsigned int>>,
4427 cell_hash<dim, spacedim>> &output_unmap,
4428 const bool & check_owned)
4430 bool no_check = !check_owned;
4434 for (
auto const &rank_and_points : map_pts)
4437 const auto &received_process = rank_and_points.first;
4438 const auto &received_points = rank_and_points.second.first;
4439 const auto &received_map = rank_and_points.second.second;
4444 typename Triangulation<dim, spacedim>::active_cell_iterator>
4446 std::vector<std::vector<Point<dim>>> in_qpoints;
4447 std::vector<std::vector<unsigned int>> in_maps;
4448 std::vector<std::vector<Point<spacedim>>> in_points;
4451 compute_point_locations_unmap(cache,
4452 rank_and_points.second.first);
4453 for (
const auto &map_c_pt_idx : cpt_loc_pts)
4456 const auto &proc_cell = map_c_pt_idx.first;
4457 const auto &proc_qpoints = map_c_pt_idx.second.first;
4458 const auto &proc_maps = map_c_pt_idx.second.second;
4462 if (no_check || proc_cell->is_locally_owned())
4464 in_cell.emplace_back(proc_cell);
4465 in_qpoints.emplace_back(proc_qpoints);
4467 unsigned int loc_size = proc_qpoints.size();
4468 std::vector<unsigned int> cell_maps(loc_size);
4469 std::vector<Point<spacedim>> cell_points(loc_size);
4470 for (
unsigned int pt = 0; pt < loc_size; ++pt)
4472 cell_maps[pt] = received_map[proc_maps[pt]];
4473 cell_points[pt] = received_points[proc_maps[pt]];
4475 in_maps.emplace_back(cell_maps);
4476 in_points.emplace_back(cell_points);
4481 internal::distributed_cptloc::merge_cptloc_outputs(
4495 template <
int dim,
int spacedim>
4498 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
4499 std::vector<std::vector<Point<dim>>>,
4500 std::vector<std::vector<unsigned int>>,
4501 std::vector<std::vector<Point<spacedim>>>,
4502 std::vector<std::vector<unsigned int>>>
4511 #ifndef DEAL_II_WITH_MPI 4514 (void)global_bboxes;
4517 "GridTools::distributed_compute_point_locations() requires MPI."));
4519 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
4520 std::vector<std::vector<Point<dim>>>,
4521 std::vector<std::vector<unsigned int>>,
4522 std::vector<std::vector<Point<spacedim>>>,
4523 std::vector<std::vector<unsigned int>>>
4528 const auto &tria_mpi =
4536 "GridTools::distributed_compute_point_locations() requires a parallel triangulation."));
4537 auto mpi_communicator = tria_mpi->get_communicator();
4540 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
4541 std::vector<std::vector<Point<dim>>>,
4542 std::vector<std::vector<unsigned int>>,
4543 std::vector<std::vector<Point<spacedim>>>,
4544 std::vector<std::vector<unsigned int>>>
4550 std::tuple<std::vector<Point<dim>>,
4551 std::vector<unsigned int>,
4552 std::vector<Point<spacedim>>,
4553 std::vector<unsigned int>>,
4554 internal::distributed_cptloc::cell_hash<dim, spacedim>>
4562 std::tuple<std::vector<std::vector<unsigned int>>,
4563 std::map<unsigned int, unsigned int>,
4564 std::map<unsigned int, std::vector<unsigned int>>>
4570 const auto &guess_loc_idx = std::get<0>(guessed_points)[my_rank];
4571 const unsigned int n_local_guess = guess_loc_idx.size();
4573 std::vector<Point<spacedim>> guess_local_pts(n_local_guess);
4574 for (
unsigned int i = 0; i < n_local_guess; ++i)
4575 guess_local_pts[i] = local_points[guess_loc_idx[i]];
4578 std::map<
unsigned int,
4579 std::tuple<std::vector<CellId>,
4580 std::vector<std::vector<Point<dim>>>,
4581 std::vector<std::vector<unsigned int>>,
4582 std::vector<std::vector<Point<spacedim>>>>>
4586 std::vector<unsigned int> classified_pts;
4590 &internal::distributed_cptloc::compute_and_classify_points<dim, spacedim>,
4600 const auto &other_owned_idx = std::get<1>(guessed_points);
4601 std::map<
unsigned int,
4602 std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>>
4605 for (
const auto &indices : other_owned_idx)
4606 if (indices.second != my_rank)
4609 auto ¤t_pts = other_owned_pts[indices.second];
4611 current_pts.first.emplace_back(local_points[indices.first]);
4612 current_pts.second.emplace_back(indices.first);
4616 auto owned_rank_pts =
4623 &internal::distributed_cptloc::compute_and_merge_from_map<dim, spacedim>,
4637 std::map<
unsigned int,
4638 std::pair<std::vector<Point<spacedim>>, std::vector<unsigned int>>>
4643 const std::map<unsigned int, std::vector<unsigned int>> &other_check_idx =
4644 std::get<2>(guessed_points);
4651 std::sort(classified_pts.begin(), classified_pts.end());
4653 for (
const auto &pt_to_guesses : other_check_idx)
4655 const auto &point_idx = pt_to_guesses.first;
4656 const auto &probable_owners_rks = pt_to_guesses.second;
4657 if (!std::binary_search(classified_pts.begin(),
4658 classified_pts.end(),
4662 for (
unsigned int i = 0; i < probable_owners_rks.size(); ++i)
4663 if (probable_owners_rks[i] != my_rank)
4666 auto ¤t_pts = other_check_pts[probable_owners_rks[i]];
4668 current_pts.first.emplace_back(local_points[point_idx]);
4670 current_pts.second.emplace_back(point_idx);
4678 owned_pts_tsk.
join();
4681 for (
const auto &rank_vals : cpt_ghost)
4684 const auto &cell_ids = std::get<0>(rank_vals.second);
4685 unsigned int n_cells = cell_ids.size();
4686 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
4688 for (
unsigned int c = 0; c < n_cells; ++c)
4691 internal::distributed_cptloc::merge_cptloc_outputs(
4694 std::get<1>(rank_vals.second),
4695 std::get<2>(rank_vals.second),
4696 std::get<3>(rank_vals.second),
4702 internal::distributed_cptloc::compute_and_merge_from_map(cache,
4709 unsigned int size_output = temporary_unmap.size();
4710 auto &out_cells = std::get<0>(output_tuple);
4711 auto &out_qpoints = std::get<1>(output_tuple);
4712 auto &out_maps = std::get<2>(output_tuple);
4713 auto &out_points = std::get<3>(output_tuple);
4714 auto &out_ranks = std::get<4>(output_tuple);
4716 out_cells.resize(size_output);
4717 out_qpoints.resize(size_output);
4718 out_maps.resize(size_output);
4719 out_points.resize(size_output);
4720 out_ranks.resize(size_output);
4723 for (
const auto &rank_and_tuple : temporary_unmap)
4725 out_cells[c] = rank_and_tuple.first;
4726 out_qpoints[c] = std::get<0>(rank_and_tuple.second);
4727 out_maps[c] = std::get<1>(rank_and_tuple.second);
4728 out_points[c] = std::get<2>(rank_and_tuple.second);
4729 out_ranks[c] = std::get<3>(rank_and_tuple.second);
4733 return output_tuple;
4738 template <
int dim,
int spacedim>
4739 std::map<unsigned int, Point<spacedim>>
4743 std::map<unsigned int, Point<spacedim>> result;
4746 if (!cell->is_artificial())
4749 for (
unsigned int i = 0; i < vs.size(); ++i)
4750 result[cell->vertex_index(i)] = vs[i];
4757 template <
int spacedim>
4762 auto id_and_v = std::min_element(
4765 [&](
const std::pair<const unsigned int, Point<spacedim>> &p1,
4766 const std::pair<const unsigned int, Point<spacedim>> &p2) ->
bool {
4767 return p1.second.
distance(p) < p2.second.distance(p);
4769 return id_and_v->first;
4773 template <
int dim,
int spacedim>
4774 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
4781 const std::vector<bool> &marked_vertices)
4786 const auto &vertex_to_cell_centers =
4793 vertex_to_cell_centers,
4798 template <
int spacedim>
4799 std::vector<std::vector<BoundingBox<spacedim>>>
4800 exchange_local_bounding_boxes(
4802 MPI_Comm mpi_communicator)
4804 #ifndef DEAL_II_WITH_MPI 4806 (void)mpi_communicator;
4809 "GridTools::exchange_local_bounding_boxes() requires MPI."));
4813 unsigned int n_bboxes = local_bboxes.size();
4815 int n_local_data = 2 * spacedim * n_bboxes;
4817 std::vector<double> loc_data_array(n_local_data);
4818 for (
unsigned int i = 0; i < n_bboxes; ++i)
4819 for (
unsigned int d = 0; d < spacedim; ++d)
4822 loc_data_array[2 * i * spacedim + d] =
4823 local_bboxes[i].get_boundary_points().first[d];
4824 loc_data_array[2 * i * spacedim + spacedim + d] =
4825 local_bboxes[i].get_boundary_points().second[d];
4832 std::vector<int> size_all_data(n_procs);
4835 int ierr = MPI_Allgather(&n_local_data,
4838 &(size_all_data[0]),
4846 std::vector<int> rdispls(n_procs);
4848 for (
unsigned int i = 1; i < n_procs; ++i)
4849 rdispls[i] = rdispls[i - 1] + size_all_data[i - 1];
4853 std::vector<double> data_array(rdispls.back() + size_all_data.back());
4855 ierr = MPI_Allgatherv(&(loc_data_array[0]),
4859 &(size_all_data[0]),
4866 std::vector<std::vector<BoundingBox<spacedim>>> global_bboxes(n_procs);
4867 unsigned int begin_idx = 0;
4868 for (
unsigned int i = 0; i < n_procs; ++i)
4871 unsigned int n_bbox_i = size_all_data[i] / (spacedim * 2);
4872 global_bboxes[i].resize(n_bbox_i);
4873 for (
unsigned int bbox = 0; bbox < n_bbox_i; ++bbox)
4876 for (
unsigned int d = 0; d < spacedim; ++d)
4878 p1[d] = data_array[begin_idx + 2 * bbox * spacedim + d];
4880 data_array[begin_idx + 2 * bbox * spacedim + spacedim + d];
4883 global_bboxes[i][bbox] = loc_bbox;
4886 begin_idx += size_all_data[i];
4888 return global_bboxes;
4889 #endif // DEAL_II_WITH_MPI 4897 #include "grid_tools.inst" 4899 DEAL_II_NAMESPACE_CLOSE
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_raw_faces() const
static void reorder_cells(std::vector< CellData< dim >> &original_cells, const bool use_new_style_ordering=false)
static const unsigned int invalid_unsigned_int
virtual std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > get_vertices(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
active_cell_iterator begin_active(const unsigned int level=0) const
#define AssertDimension(dim1, dim2)
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
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 >())
const std::vector< Point< spacedim > > & get_vertices() const
types::subdomain_id locally_owned_subdomain() const override
cell_iterator end() const
active_face_iterator begin_active_face() const
Task< RT > new_task(const std::function< RT()> &function)
bool check_consistency(const unsigned int dim) const
IteratorRange< active_cell_iterator > active_cell_iterators() const
void add(const size_type i, const size_type j)
const std::vector< bool > & get_used_vertices() const
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
numbers::NumberTraits< Number >::real_type norm() 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
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
cell_iterator begin(const unsigned int level=0) const
unsigned int n_active_cells() const
unsigned long long int global_dof_index
unsigned int n_levels() const
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
static double distance_to_unit_cell(const Point< dim > &p)
virtual void execute_coarsening_and_refinement()
virtual MPI_Comm get_communicator() const
unsigned int size() const
static::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
unsigned int n_cells() const
T sum(const T &t, const MPI_Comm &mpi_communicator)
virtual void create_triangulation(const std::vector< Point< spacedim >> &vertices, const std::vector< CellData< dim >> &cells, const SubCellData &subcelldata)
#define Assert(cond, exc)
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
types::global_dof_index n_dofs() const
std::vector< types::manifold_id > get_manifold_ids() 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)
Abstract base class for mapping classes.
std::list< typename Triangulation< dim, spacedim >::cell_iterator > distorted_cells
virtual bool has_hanging_nodes() const
types::material_id material_id
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
std::vector< types::boundary_id > get_boundary_ids() const
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
IteratorState::IteratorStates state() const
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
std::vector< unsigned int > invert_permutation(const std::vector< unsigned int > &permutation)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
cell_iterator end() const
bool vertex_used(const unsigned int index) const
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access >> &cell)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
types::manifold_id manifold_id
void add_constraints(const ConstraintList &new_constraints)
numbers::NumberTraits< Number >::real_type square() const
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
void swap(Vector< Number > &u, Vector< Number > &v)
static void alternating_form_at_vertices(const Point< spacedim >(&vertices)[vertices_per_cell], Tensor< spacedim-dim, spacedim >(&forms)[vertices_per_cell])
const types::subdomain_id artificial_subdomain_id
const types::manifold_id invalid_manifold_id
#define AssertThrowMPI(error_code)
face_iterator end_face() const
T min(const T &t, const MPI_Comm &mpi_communicator)
std::vector< CellData< 2 > > boundary_quads
virtual void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
void refine_global(const unsigned int times=1)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
virtual bool preserves_vertex_locations() const =0
static::ExceptionBase & ExcNotImplemented()
Iterator points to a valid object.
static::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
active_cell_iterator begin_active(const unsigned int level=0) const
unsigned long long int global_vertex_index
void make_sparsity_pattern(const DoFHandlerType &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)
unsigned int n_vertices() const
std::map< unsigned int, T > some_to_some(const MPI_Comm &comm, const std::map< unsigned int, T > &objects_to_send)
T max(const T &t, const MPI_Comm &mpi_communicator)
double JxW(const unsigned int quadrature_point) const
void apply_constraints(VectorType &v, const bool matrix_is_symmetric) const
unsigned int vertices[GeometryInfo< structdim >::vertices_per_cell]
static::ExceptionBase & ExcInternalError()
Triangulation< dim, spacedim > & get_triangulation()