16 #ifndef dealii_grid_tools_h 17 # define dealii_grid_tools_h 20 # include <deal.II/base/config.h> 22 # include <deal.II/base/bounding_box.h> 23 # include <deal.II/base/geometry_info.h> 25 # include <deal.II/dofs/dof_handler.h> 27 # include <deal.II/fe/mapping.h> 28 # include <deal.II/fe/mapping_q1.h> 30 # include <deal.II/grid/manifold.h> 31 # include <deal.II/grid/tria.h> 32 # include <deal.II/grid/tria_accessor.h> 33 # include <deal.II/grid/tria_iterator.h> 35 # include <deal.II/hp/dof_handler.h> 37 # include <deal.II/lac/sparsity_tools.h> 39 # include <boost/archive/binary_iarchive.hpp> 40 # include <boost/archive/binary_oarchive.hpp> 41 # include <boost/optional.hpp> 42 # include <boost/serialization/array.hpp> 43 # include <boost/serialization/vector.hpp> 45 # ifdef DEAL_II_WITH_ZLIB 46 # include <boost/iostreams/device/back_inserter.hpp> 47 # include <boost/iostreams/filter/gzip.hpp> 48 # include <boost/iostreams/filtering_stream.hpp> 49 # include <boost/iostreams/stream.hpp> 57 DEAL_II_NAMESPACE_OPEN
71 class MappingCollection;
78 template <
int dim,
int spacedim,
class MeshType>
79 class ActiveCellIterator
83 using type =
typename MeshType::active_cell_iterator;
90 template <
int dim,
int spacedim>
91 class ActiveCellIterator<dim, spacedim, ::
DoFHandler<dim, spacedim>>
98 template <
int dim,
int spacedim>
99 class ActiveCellIterator<dim, spacedim, ::
hp::DoFHandler<dim, spacedim>>
118 template <
int dim,
int spacedim>
132 template <
int dim,
int spacedim>
162 template <
int dim,
int spacedim>
178 template <
int dim,
int spacedim>
194 template <
int dim,
int spacedim>
220 template <
int dim,
typename T>
237 template <
int dim,
int spacedim>
260 template <
typename Iterator>
263 const Iterator &
object,
288 template <
int dim,
int spacedim>
310 template <
int dim,
int spacedim>
315 std::vector<unsigned int> & considered_vertices,
316 const double tol = 1e-12);
347 template <
int dim,
typename Transformation,
int spacedim>
349 transform(
const Transformation & transformation,
357 template <
int dim,
int spacedim>
387 rotate(
const double angle,
388 const unsigned int axis,
453 const bool solve_for_absolute_positions =
false);
460 template <
int dim,
int spacedim>
461 std::map<unsigned int, Point<spacedim>>
471 template <
int dim,
int spacedim>
473 scale(
const double scaling_factor,
486 template <
int dim,
int spacedim>
490 const bool keep_boundary =
true);
527 template <
int dim,
int spacedim>
530 const bool isotropic =
false,
531 const unsigned int max_iterations = 100);
559 template <
int dim,
int spacedim>
562 const double max_ratio = 1.6180339887,
563 const unsigned int max_iterations = 5);
656 template <
int dim,
int spacedim>
659 const double limit_angle_fraction = .75);
708 template <
int dim,
int spacedim>
711 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
712 std::vector<std::vector<Point<dim>>>,
713 std::vector<std::vector<unsigned int>>>
788 template <
int dim,
int spacedim>
791 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
792 std::vector<std::vector<Point<dim>>>,
793 std::vector<std::vector<unsigned int>>,
794 std::vector<std::vector<Point<spacedim>>>,
795 std::vector<std::vector<unsigned int>>>
842 template <
int dim,
int spacedim>
843 std::map<unsigned int, Point<spacedim>>
859 template <
int spacedim>
890 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
893 const MeshType<dim, spacedim> &mesh,
895 const std::vector<bool> & marked_vertices = std::vector<bool>());
922 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
926 const MeshType<dim, spacedim> &mesh,
928 const std::vector<bool> & marked_vertices = std::vector<bool>());
955 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
957 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
960 typename ::internal::
961 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
964 const unsigned int vertex_index);
991 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
993 typename MeshType<dim, spacedim>::active_cell_iterator
995 typename ::internal::
996 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
999 const MeshType<dim, spacedim> &mesh,
1001 const std::vector<bool> & marked_vertices = std::vector<bool>());
1089 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1091 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1093 std::pair<typename ::internal::
1094 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1099 const MeshType<dim, spacedim> &mesh,
1101 const std::vector<bool> & marked_vertices = std::vector<bool>());
1112 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1114 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1116 std::pair<typename ::internal::
1117 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1122 const MeshType<dim, spacedim> &mesh,
1125 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
1126 & vertex_to_cell_map,
1128 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1129 typename MeshType<dim, spacedim>::active_cell_iterator(),
1130 const std::vector<bool> &marked_vertices = std::vector<bool>());
1153 template <
int dim,
int spacedim>
1154 std::pair<typename hp::DoFHandler<dim, spacedim>::active_cell_iterator,
1167 template <
int dim,
int spacedim>
1168 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1175 const std::vector<bool> &marked_vertices = std::vector<bool>());
1194 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1196 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1199 std::vector<std::pair<
1200 typename ::internal::
1201 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1206 const MeshType<dim, spacedim> &mesh,
1208 const double tolerance = 1e-12,
1209 const std::vector<bool> & marked_vertices = std::vector<bool>());
1232 template <
class MeshType>
1233 std::vector<typename MeshType::active_cell_iterator>
1246 template <
class MeshType>
1249 const typename MeshType::active_cell_iterator & cell,
1250 std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1304 template <
class MeshType>
1305 std::vector<typename MeshType::active_cell_iterator>
1307 const MeshType &mesh,
1308 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1319 template <
class MeshType>
1320 std::vector<typename MeshType::cell_iterator>
1322 const MeshType &mesh,
1323 const std::function<
bool(
const typename MeshType::cell_iterator &)>
1325 const unsigned int level);
1343 template <
class MeshType>
1344 std::vector<typename MeshType::active_cell_iterator>
1398 template <
class MeshType>
1399 std::vector<typename MeshType::active_cell_iterator>
1401 const MeshType &mesh,
1402 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1404 const double layer_thickness);
1431 template <
class MeshType>
1432 std::vector<typename MeshType::active_cell_iterator>
1434 const double layer_thickness);
1451 template <
class MeshType>
1454 const MeshType &mesh,
1455 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1510 template <
class MeshType>
1511 std::vector<BoundingBox<MeshType::space_dimension>>
1513 const MeshType &mesh,
1514 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1516 const unsigned int refinement_level = 0,
1517 const bool allow_merge =
false,
1549 template <
int spacedim>
1551 std::tuple<std::vector<std::vector<unsigned int>>,
1552 std::map<unsigned int, unsigned int>,
1553 std::map<unsigned int, std::vector<unsigned int>>>
1570 template <
int dim,
int spacedim>
1572 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1589 template <
int dim,
int spacedim>
1590 std::vector<std::vector<Tensor<1, spacedim>>>
1604 template <
int dim,
int spacedim>
1621 template <
int dim,
int spacedim>
1622 std::map<unsigned int, types::global_vertex_index>
1639 template <
int dim,
int spacedim>
1640 std::pair<unsigned int, double>
1658 template <
int dim,
int spacedim>
1672 template <
int dim,
int spacedim>
1686 template <
int dim,
int spacedim>
1690 const unsigned int level,
1713 template <
int dim,
int spacedim>
1730 template <
int dim,
int spacedim>
1733 const std::vector<unsigned int> &cell_weights,
1783 template <
int dim,
int spacedim>
1801 template <
int dim,
int spacedim>
1804 const std::vector<unsigned int> &cell_weights,
1817 template <
int dim,
int spacedim>
1833 template <
int dim,
int spacedim>
1847 template <
int dim,
int spacedim>
1850 std::vector<types::subdomain_id> & subdomain);
1866 template <
int dim,
int spacedim>
1902 template <
int dim,
int spacedim>
1940 template <
typename MeshType>
1941 std::list<std::pair<
typename MeshType::cell_iterator,
1942 typename MeshType::cell_iterator>>
1954 template <
int dim,
int spacedim>
1968 template <
typename MeshType>
1993 template <
int dim,
int spacedim>
2049 template <
class MeshType>
2050 std::vector<typename MeshType::active_cell_iterator>
2077 template <
class Container>
2078 std::vector<typename Container::cell_iterator>
2080 const std::vector<typename Container::active_cell_iterator> &patch_cells);
2150 template <
class Container>
2153 const std::vector<typename Container::active_cell_iterator> &patch,
2155 &local_triangulation,
2158 Container::space_dimension>::active_cell_iterator,
2159 typename Container::active_cell_iterator> &patch_to_global_tria_map);
2197 template <
class DoFHandlerType>
2199 std::vector<typename DoFHandlerType::active_cell_iterator>>
2215 template <
typename CellIterator>
2221 CellIterator cell[2];
2227 unsigned int face_idx[2];
2317 template <
typename FaceIterator>
2319 std::bitset<3> & orientation,
2320 const FaceIterator & face1,
2321 const FaceIterator & face2,
2322 const int direction,
2331 template <
typename FaceIterator>
2334 const FaceIterator & face1,
2335 const FaceIterator & face2,
2336 const int direction,
2400 template <
typename MeshType>
2403 const MeshType & mesh,
2406 const int direction,
2436 template <
typename MeshType>
2439 const MeshType & mesh,
2441 const int direction,
2444 const ::Tensor<1, MeshType::space_dimension> &offset =
2476 template <
int dim,
int spacedim>
2479 const bool reset_boundary_ids =
false);
2504 template <
int dim,
int spacedim>
2507 const std::vector<types::boundary_id> &src_boundary_ids,
2508 const std::vector<types::manifold_id> &dst_manifold_ids,
2510 const std::vector<types::boundary_id> &reset_boundary_ids = {});
2543 template <
int dim,
int spacedim>
2546 const bool compute_face_ids =
false);
2633 template <
typename DataType,
typename MeshType>
2636 const MeshType & mesh,
2637 const std::function<boost::optional<DataType>(
2638 const typename MeshType::active_cell_iterator &)> &pack,
2639 const std::function<
void(
const typename MeshType::active_cell_iterator &,
2640 const DataType &)> & unpack);
2653 template <
int spacedim>
2654 std::vector<std::vector<BoundingBox<spacedim>>>
2655 exchange_local_bounding_boxes(
2657 MPI_Comm mpi_communicator);
2671 template <
int dim,
typename T>
2691 template <
class Archive>
2693 save(Archive &ar,
const unsigned int version)
const;
2699 template <
class Archive>
2701 load(Archive &ar,
const unsigned int version);
2703 BOOST_SERIALIZATION_SPLIT_MEMBER()
2716 <<
"The number of partitions you gave is " << arg1
2717 <<
", but must be greater than zero.");
2723 <<
"The subdomain id " << arg1
2724 <<
" has no cells associated with it.");
2735 <<
"The scaling factor must be positive, but it is " << arg1
2743 <<
"The point <" << arg1
2744 <<
"> could not be found inside any of the " 2745 <<
"coarse grid cells.");
2752 <<
"The point <" << arg1
2753 <<
"> could not be found inside any of the " 2754 <<
"subcells of a coarse grid cell.");
2761 <<
"The given vertex with index " << arg1
2762 <<
" is not used in the given triangulation.");
2777 template <
int dim,
typename T>
2782 return std::numeric_limits<double>::quiet_NaN();
2785 template <
int dim,
typename Predicate,
int spacedim>
2790 std::vector<bool> treated_vertices(triangulation.
n_vertices(),
false);
2800 endc = triangulation.
end();
2801 for (; cell != endc; ++cell)
2802 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2803 if (treated_vertices[cell->vertex_index(v)] ==
false)
2806 cell->vertex(v) = predicate(cell->vertex(v));
2808 treated_vertices[cell->vertex_index(v)] =
true;
2817 endc = triangulation.
end();
2818 for (; cell != endc; ++cell)
2819 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2821 if (cell->face(face)->has_children() &&
2822 !cell->face(face)->at_boundary())
2825 cell->face(face)->child(0)->vertex(1) =
2826 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
2834 endc = triangulation.
end();
2835 for (; cell != endc; ++cell)
2836 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
2838 if (cell->face(face)->has_children() &&
2839 !cell->face(face)->at_boundary())
2842 cell->face(face)->child(0)->vertex(1) =
2843 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
2845 cell->face(face)->child(0)->vertex(2) =
2846 (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
2848 cell->face(face)->child(1)->vertex(3) =
2849 (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
2851 cell->face(face)->child(2)->vertex(3) =
2852 (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
2856 cell->face(face)->child(0)->vertex(3) =
2857 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
2858 cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
2864 triangulation.
signals.mesh_movement();
2869 template <
class MeshType>
2870 std::vector<typename MeshType::active_cell_iterator>
2873 std::vector<typename MeshType::active_cell_iterator> child_cells;
2875 if (cell->has_children())
2877 for (
unsigned int child = 0; child < cell->n_children(); ++child)
2878 if (cell->child(child)->has_children())
2880 const std::vector<typename MeshType::active_cell_iterator>
2881 children = get_active_child_cells<MeshType>(cell->child(child));
2882 child_cells.insert(child_cells.end(),
2887 child_cells.push_back(cell->child(child));
2895 template <
class MeshType>
2898 const typename MeshType::active_cell_iterator & cell,
2899 std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
2901 active_neighbors.clear();
2902 for (
unsigned int n = 0;
2903 n < GeometryInfo<MeshType::dimension>::faces_per_cell;
2905 if (!cell->at_boundary(n))
2907 if (MeshType::dimension == 1)
2914 typename MeshType::cell_iterator neighbor_child =
2916 if (!neighbor_child->active())
2918 while (neighbor_child->has_children())
2919 neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
2921 Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
2924 active_neighbors.push_back(neighbor_child);
2928 if (cell->face(n)->has_children())
2932 for (
unsigned int c = 0;
2933 c < cell->face(n)->number_of_children();
2935 active_neighbors.push_back(
2936 cell->neighbor_child_on_subface(n, c));
2942 active_neighbors.push_back(cell->neighbor(n));
2952 namespace ProjectToObject
2966 struct CrossDerivative
2968 const unsigned int direction_0;
2969 const unsigned int direction_1;
2971 CrossDerivative(
const unsigned int d0,
const unsigned int d1);
2974 inline CrossDerivative::CrossDerivative(
const unsigned int d0,
2975 const unsigned int d1)
2986 template <
typename F>
2988 centered_first_difference(
const double center,
2990 const F &f) -> decltype(f(center) - f(center))
2992 return (f(center + step) - f(center - step)) / (2.0 * step);
3001 template <
typename F>
3003 centered_second_difference(
const double center,
3005 const F &f) -> decltype(f(center) - f(center))
3007 return (f(center + step) - 2.0 * f(center) + f(center - step)) /
3022 template <
int structdim,
typename F>
3025 const CrossDerivative cross_derivative,
3028 const F &f) -> decltype(f(center) - f(center))
3031 simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3032 simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3033 return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
3034 1.0 / 3.0 * f(center - simplex_vector) +
3035 16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
3047 template <
int spacedim,
int structdim,
typename F>
3050 const unsigned int row_n,
3051 const unsigned int dependent_direction,
3058 dependent_direction <
3060 ExcMessage(
"This function assumes that the last weight is a " 3061 "dependent variable (and hence we cannot take its " 3062 "derivative directly)."));
3063 Assert(row_n != dependent_direction,
3065 "We cannot differentiate with respect to the variable " 3066 "that is assumed to be dependent."));
3070 {row_n, dependent_direction}, center, step, f);
3072 for (
unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3074 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3083 template <
typename Iterator,
int spacedim,
int structdim>
3085 project_to_d_linear_object(
const Iterator &
object,
3120 for (
unsigned int d = 0; d < structdim; ++d)
3124 for (
unsigned int i = 0; i < GeometryInfo<structdim>::vertices_per_cell;
3126 x_k += object->vertex(i) *
3132 for (
unsigned int i = 0;
3133 i < GeometryInfo<structdim>::vertices_per_cell;
3136 (x_k - trial_point) *
object->vertex(i) *
3141 for (
unsigned int i = 0;
3142 i < GeometryInfo<structdim>::vertices_per_cell;
3144 for (
unsigned int j = 0;
3145 j < GeometryInfo<structdim>::vertices_per_cell;
3153 H_k += (
object->vertex(i) *
object->vertex(j)) * tmp;
3160 for (
unsigned int i = 0;
3161 i < GeometryInfo<structdim>::vertices_per_cell;
3163 x_k += object->vertex(i) *
3166 if (delta_xi.
norm() < 1e-7)
3182 template <
int structdim>
3189 static const std::size_t n_vertices_per_cell =
3191 n_independent_components;
3192 std::array<double, n_vertices_per_cell> copied_weights;
3193 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
3195 copied_weights[i] = v[i];
3196 if (v[i] < 0.0 || v[i] > 1.0)
3201 std::sort(copied_weights.begin(), copied_weights.end());
3203 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
3204 return std::abs(sum - 1.0) < 1e-10;
3208 template <
typename Iterator>
3211 const Iterator &
object,
3214 const int spacedim = Iterator::AccessorType::space_dimension;
3215 const int structdim = Iterator::AccessorType::structure_dimension;
3219 if (structdim >= spacedim)
3220 return projected_point;
3221 else if (structdim == 1 || structdim == 2)
3223 using namespace internal::ProjectToObject;
3228 const int dim = Iterator::AccessorType::dimension;
3231 &manifold) !=
nullptr)
3234 project_to_d_linear_object<Iterator, spacedim, structdim>(
3235 object, trial_point);
3280 const double step_size =
object->diameter() / 64.0;
3282 constexpr
unsigned int n_vertices_per_cell =
3285 std::array<Point<spacedim>, n_vertices_per_cell> vertices;
3286 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3288 vertices[vertex_n] = object->vertex(vertex_n);
3290 auto get_point_from_weights =
3293 return object->get_manifold().get_new_point(
3294 make_array_view(vertices.begin(), vertices.end()),
3295 make_array_view(&weights[0],
3296 &weights[n_vertices_per_cell - 1] + 1));
3302 double guess_weights_sum = 0.0;
3303 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
3306 const double distance =
3307 vertices[vertex_n].distance(trial_point);
3308 if (distance == 0.0)
3310 guess_weights = 0.0;
3311 guess_weights[vertex_n] = 1.0;
3312 guess_weights_sum = 1.0;
3317 guess_weights[vertex_n] = 1.0 / distance;
3318 guess_weights_sum += guess_weights[vertex_n];
3321 guess_weights /= guess_weights_sum;
3322 Assert(internal::weights_are_ok<structdim>(guess_weights),
3331 for (
unsigned int outer_n = 0; outer_n < 40; ++outer_n)
3333 const unsigned int dependent_direction =
3334 n_vertices_per_cell - 1;
3336 for (
unsigned int row_n = 0; row_n < n_vertices_per_cell;
3339 if (row_n != dependent_direction)
3341 current_gradient[row_n] =
3342 gradient_entry<spacedim, structdim>(
3344 dependent_direction,
3348 get_point_from_weights);
3350 current_gradient[dependent_direction] -=
3351 current_gradient[row_n];
3369 double gradient_weight = -0.5;
3370 auto gradient_weight_objective_function =
3371 [&](
const double gradient_weight_guess) ->
double {
3372 return (trial_point -
3373 get_point_from_weights(guess_weights +
3374 gradient_weight_guess *
3379 for (
unsigned int inner_n = 0; inner_n < 10; ++inner_n)
3381 const double update_numerator = centered_first_difference(
3384 gradient_weight_objective_function);
3385 const double update_denominator =
3386 centered_second_difference(
3389 gradient_weight_objective_function);
3393 if (std::abs(update_denominator) == 0.0)
3396 gradient_weight - update_numerator / update_denominator;
3401 if (std::abs(gradient_weight) > 10)
3403 gradient_weight = -10.0;
3413 guess_weights + gradient_weight * current_gradient;
3415 double new_gradient_weight = gradient_weight;
3416 for (
unsigned int iteration_count = 0; iteration_count < 40;
3419 if (internal::weights_are_ok<structdim>(tentative_weights))
3422 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
3424 if (tentative_weights[i] < 0.0)
3426 tentative_weights -=
3427 (tentative_weights[i] / current_gradient[i]) *
3430 if (tentative_weights[i] < 0.0 ||
3431 1.0 < tentative_weights[i])
3433 new_gradient_weight /= 2.0;
3436 new_gradient_weight * current_gradient;
3443 if (!internal::weights_are_ok<structdim>(tentative_weights))
3448 if (get_point_from_weights(tentative_weights)
3449 .distance(trial_point) <
3450 get_point_from_weights(guess_weights).distance(trial_point))
3451 guess_weights = tentative_weights;
3454 Assert(internal::weights_are_ok<structdim>(guess_weights),
3457 Assert(internal::weights_are_ok<structdim>(guess_weights),
3459 projected_point = get_point_from_weights(guess_weights);
3469 for (
unsigned int line_n = 0;
3470 line_n < GeometryInfo<structdim>::lines_per_cell;
3473 line_projections[line_n] =
3476 std::sort(line_projections.begin(),
3477 line_projections.end(),
3480 b.distance(trial_point);
3482 if (line_projections[0].distance(trial_point) <
3483 projected_point.
distance(trial_point))
3484 projected_point = line_projections[0];
3490 return projected_point;
3493 return projected_point;
3498 template <
int dim,
typename T>
3499 template <
class Archive>
3502 const unsigned int )
const 3507 const size_t n_cells =
cell_ids.size();
3512 ar & binary_cell_id;
3520 template <
int dim,
typename T>
3521 template <
class Archive>
3524 const unsigned int )
3530 for (
unsigned int c = 0; c < n_cells; ++c)
3541 template <
typename DataType,
typename MeshType>
3544 const MeshType & mesh,
3545 const std::function<boost::optional<DataType>(
3546 const typename MeshType::active_cell_iterator &)> &pack,
3547 const std::function<
void(
const typename MeshType::active_cell_iterator &,
3548 const DataType &)> & unpack)
3550 # ifndef DEAL_II_WITH_MPI 3556 "GridTools::exchange_cell_data_to_ghosts() requires MPI."));
3558 constexpr
int dim = MeshType::dimension;
3559 constexpr
int spacedim = MeshType::space_dimension;
3561 &mesh.get_triangulation());
3565 "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
3568 using DestinationToBufferMap =
3571 DestinationToBufferMap destination_to_data_buffer_map;
3573 std::map<unsigned int, std::set<::types::subdomain_id>>
3574 vertices_with_ghost_neighbors =
3575 tria->compute_vertices_with_ghost_neighbors();
3578 if (cell->is_locally_owned())
3580 std::set<::types::subdomain_id> send_to;
3581 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
3584 const std::map<
unsigned int,
3585 std::set<::types::subdomain_id>>::
3586 const_iterator neighbor_subdomains_of_vertex =
3587 vertices_with_ghost_neighbors.find(cell->vertex_index(v));
3589 if (neighbor_subdomains_of_vertex ==
3590 vertices_with_ghost_neighbors.end())
3593 Assert(neighbor_subdomains_of_vertex->second.size() != 0,
3596 send_to.insert(neighbor_subdomains_of_vertex->second.begin(),
3597 neighbor_subdomains_of_vertex->second.end());
3600 if (send_to.size() > 0)
3603 typename MeshType::active_cell_iterator mesh_it(tria,
3608 const boost::optional<DataType>
data = pack(mesh_it);
3612 const CellId cellid = cell->id();
3614 for (
auto it : send_to)
3616 const ::types::subdomain_id subdomain = it;
3620 typename DestinationToBufferMap::iterator p =
3621 destination_to_data_buffer_map
3622 .insert(std::make_pair(
3623 subdomain, CellDataTransferBuffer<dim, DataType>()))
3626 p->second.cell_ids.emplace_back(cellid);
3627 p->second.data.emplace_back(data.get());
3635 std::set<::types::subdomain_id> ghost_owners = tria->ghost_owners();
3636 const unsigned int n_ghost_owners = ghost_owners.size();
3637 std::vector<std::vector<char>> sendbuffers(n_ghost_owners);
3638 std::vector<MPI_Request> requests(n_ghost_owners);
3640 unsigned int idx = 0;
3641 for (
auto it = ghost_owners.begin(); it != ghost_owners.end(); ++it, ++idx)
3643 CellDataTransferBuffer<dim, DataType> &data =
3644 destination_to_data_buffer_map[*it];
3650 const int ierr = MPI_Isend(sendbuffers[idx].
data(),
3651 sendbuffers[idx].size(),
3655 tria->get_communicator(),
3661 std::vector<char> receive;
3662 for (
unsigned int idx = 0; idx < n_ghost_owners; ++idx)
3667 MPI_Probe(MPI_ANY_SOURCE, 786, tria->get_communicator(), &status);
3669 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
3672 receive.resize(len);
3674 char *ptr = receive.data();
3675 ierr = MPI_Recv(ptr,
3680 tria->get_communicator(),
3685 Utilities::unpack<CellDataTransferBuffer<dim, DataType>>(receive);
3687 DataType *data = cellinfo.
data.data();
3688 for (
unsigned int c = 0; c < cellinfo.cell_ids.size(); ++c, ++
data)
3691 tria_cell = cellinfo.cell_ids[c].to_cell(*tria);
3693 const typename MeshType::active_cell_iterator cell(
3694 tria, tria_cell->level(), tria_cell->index(), &mesh);
3696 unpack(cell, *data);
3702 if (requests.size())
3705 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
3708 # endif // DEAL_II_WITH_MPI 3714 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={})
static::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
static const unsigned int invalid_unsigned_int
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
cell_iterator end() const
static::ExceptionBase & ExcPointNotFoundInCoarseGrid(Point< N > arg1)
IteratorRange< active_cell_iterator > active_cell_iterators() const
SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
static::ExceptionBase & ExcNonExistentSubdomain(int arg1)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
numbers::NumberTraits< Number >::real_type norm() const
std::array< unsigned int, 4 > binary_type
unsigned long long int global_dof_index
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
static::ExceptionBase & ExcTriangulationHasBeenRefined()
static::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
#define DeclException1(Exception1, type1, outsequence)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
#define DeclException0(Exception0)
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, Predicate &predicate, const unsigned int grainsize)
#define AssertThrowMPI(error_code)
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static::ExceptionBase & ExcNotImplemented()
static::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
active_cell_iterator begin_active(const unsigned int level=0) const
unsigned int n_vertices() const
static::ExceptionBase & ExcPointNotFound(Point< N > arg1)
static::ExceptionBase & ExcInternalError()