16 #ifndef dealii_grid_tools_h
17 #define dealii_grid_tools_h
51 #include <boost/archive/binary_iarchive.hpp>
52 #include <boost/archive/binary_oarchive.hpp>
53 #include <boost/random/mersenne_twister.hpp>
54 #include <boost/serialization/array.hpp>
55 #include <boost/serialization/vector.hpp>
57 #ifdef DEAL_II_WITH_ZLIB
58 # include <boost/iostreams/device/back_inserter.hpp>
59 # include <boost/iostreams/filter/gzip.hpp>
60 # include <boost/iostreams/filtering_stream.hpp>
61 # include <boost/iostreams/stream.hpp>
85 class MappingCollection;
92 template <
int dim,
int spacedim>
99 template <
int dim,
int spacedim,
class MeshType>
104 using type =
typename MeshType::active_cell_iterator;
111 template <
int dim,
int spacedim>
142 template <
int dim,
int spacedim>
172 template <
int dim,
int spacedim>
176 (ReferenceCells::get_hypercube<dim>()
178 .
template get_default_linear_mapping<dim, spacedim>()
180 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
194 template <
int dim,
int spacedim>
199 (ReferenceCells::get_hypercube<dim>()
201 .
template get_default_linear_mapping<dim, spacedim>()
203 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
217 template <
int dim,
int spacedim>
222 (ReferenceCells::get_hypercube<dim>()
224 .
template get_default_linear_mapping<dim, spacedim>()
226 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
293 template <
int dim,
int spacedim>
355 template <
int dim,
int spacedim>
376 template <
typename Iterator>
379 const Iterator &
object,
394 template <
int dim,
int spacedim>
396 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>,
SubCellData>
421 template <
int dim,
int spacedim>
445 template <
int dim,
int spacedim>
450 std::vector<unsigned int> & considered_vertices,
451 const double tol = 1
e-12);
471 template <
int dim,
int spacedim>
486 template <
int dim,
int spacedim>
574 template <
int dim,
typename Transformation,
int spacedim>
585 template <
int dim,
int spacedim>
641 const unsigned int axis,
706 const bool solve_for_absolute_positions =
false);
713 template <
int dim,
int spacedim>
714 std::map<unsigned int, Point<spacedim>>
724 template <
int dim,
int spacedim>
726 scale(
const double scaling_factor,
743 template <
int dim,
int spacedim>
748 const bool keep_boundary =
true,
749 const unsigned int seed = boost::random::mt19937::default_seed);
784 template <
int dim,
int spacedim>
787 const bool isotropic =
false,
788 const unsigned int max_iterations = 100);
814 template <
int dim,
int spacedim>
817 const double max_ratio = 1.6180339887,
818 const unsigned int max_iterations = 5);
909 template <
int dim,
int spacedim>
912 const double limit_angle_fraction = .75);
974 template <
int dim,
int spacedim>
977 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
978 std::vector<std::vector<Point<dim>>>,
979 std::vector<std::vector<unsigned int>>>
1023 template <
int dim,
int spacedim>
1026 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1027 std::vector<std::vector<Point<dim>>>,
1028 std::vector<std::vector<unsigned int>>,
1029 std::vector<unsigned int>>
1109 template <
int dim,
int spacedim>
1112 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1113 std::vector<std::vector<Point<dim>>>,
1114 std::vector<std::vector<unsigned int>>,
1115 std::vector<std::vector<Point<spacedim>>>,
1116 std::vector<std::vector<unsigned int>>>
1124 const double tolerance = 1
e-10);
1142 template <
int dim,
int spacedim>
1151 std::vector<std::tuple<std::pair<int, int>,
1181 std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1206 template <
int dim,
int spacedim>
1212 const std::vector<bool> & marked_vertices,
1213 const double tolerance,
1214 const bool perform_handshake,
1215 const bool enforce_unique_mapping =
false);
1255 template <
int dim,
int spacedim>
1256 std::map<unsigned int, Point<spacedim>>
1260 (ReferenceCells::get_hypercube<dim>()
1262 .
template get_default_linear_mapping<dim, spacedim>()
1264 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1277 template <
int spacedim>
1305 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1309 const std::vector<bool> & marked_vertices = {});
1334 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1337 const MeshType<dim, spacedim> &mesh,
1339 const std::vector<bool> & marked_vertices = {});
1361 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1363 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1366 typename ::internal::
1367 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1370 const unsigned int vertex_index);
1434 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1436 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1438 std::pair<typename ::internal::
1439 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1443 const MeshType<dim, spacedim> &mesh,
1445 const std::vector<bool> &marked_vertices = {},
1446 const double tolerance = 1.e-10);
1455 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1457 typename MeshType<dim, spacedim>::active_cell_iterator
1459 typename ::internal::
1460 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1464 const std::vector<bool> &marked_vertices = {},
1465 const double tolerance = 1.e-10);
1473 template <
int dim,
int spacedim>
1474 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1480 const double tolerance = 1.e-10);
1533 template <
int dim,
int spacedim>
1534 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1537 const Cache<dim, spacedim> &cache,
1541 const std::vector<bool> &marked_vertices = {},
1542 const double tolerance = 1.e-10);
1557 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1559 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1561 std::pair<typename ::internal::
1562 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1567 const MeshType<dim, spacedim> &mesh,
1570 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
1573 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1574 typename MeshType<dim, spacedim>::active_cell_iterator(),
1575 const std::vector<bool> & marked_vertices = {},
1578 const double tolerance = 1.e-10,
1580 std::pair<BoundingBox<spacedim>,
1582 *relevant_cell_bounding_boxes_rtree =
nullptr);
1605 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1607 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1610 std::vector<std::pair<
1611 typename ::internal::
1612 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1617 const MeshType<dim, spacedim> &mesh,
1619 const double tolerance,
1620 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
1629 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1631 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1634 std::vector<std::pair<
1635 typename ::internal::
1636 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1641 const MeshType<dim, spacedim> &mesh,
1643 const double tolerance = 1
e-10,
1644 const std::vector<bool> & marked_vertices = {});
1667 template <
class MeshType>
1668 std::vector<typename MeshType::active_cell_iterator>
1695 template <
class MeshType>
1698 const typename MeshType::active_cell_iterator & cell,
1699 std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1750 template <
class MeshType>
1751 std::vector<typename MeshType::active_cell_iterator>
1753 const MeshType &mesh,
1754 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1765 template <
class MeshType>
1766 std::vector<typename MeshType::cell_iterator>
1768 const MeshType &mesh,
1769 const std::function<
bool(
const typename MeshType::cell_iterator &)>
1771 const unsigned int level);
1786 template <
class MeshType>
1787 std::vector<typename MeshType::active_cell_iterator>
1838 template <
class MeshType>
1839 std::vector<typename MeshType::active_cell_iterator>
1841 const MeshType &mesh,
1842 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1844 const double layer_thickness);
1868 template <
class MeshType>
1869 std::vector<typename MeshType::active_cell_iterator>
1871 const double layer_thickness);
1888 template <
class MeshType>
1891 const MeshType &mesh,
1892 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1947 template <
class MeshType>
1948 std::vector<BoundingBox<MeshType::space_dimension>>
1950 const MeshType &mesh,
1951 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1953 const unsigned int refinement_level = 0,
1954 const bool allow_merge =
false,
1984 template <
int spacedim>
1986 std::tuple<std::vector<std::vector<unsigned int>>,
1987 std::map<unsigned int, unsigned int>,
1988 std::map<unsigned int, std::vector<unsigned int>>>
2031 template <
int spacedim>
2033 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2034 std::map<unsigned int, unsigned int>,
2035 std::map<unsigned int, std::vector<unsigned int>>>
2052 template <
int dim,
int spacedim>
2054 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2069 template <
int dim,
int spacedim>
2070 std::vector<std::vector<Tensor<1, spacedim>>>
2085 template <
int dim,
int spacedim>
2091 (ReferenceCells::get_hypercube<dim>()
2093 .
template get_default_linear_mapping<dim, spacedim>()
2095 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2110 template <
int dim,
int spacedim>
2111 std::map<unsigned int, types::global_vertex_index>
2126 template <
int dim,
int spacedim>
2127 std::pair<unsigned int, double>
2145 template <
int dim,
int spacedim>
2159 template <
int dim,
int spacedim>
2173 template <
int dim,
int spacedim>
2177 const unsigned int level,
2200 template <
int dim,
int spacedim>
2217 template <
int dim,
int spacedim>
2220 const std::vector<unsigned int> &cell_weights,
2270 template <
int dim,
int spacedim>
2288 template <
int dim,
int spacedim>
2291 const std::vector<unsigned int> &cell_weights,
2311 template <
int dim,
int spacedim>
2315 const bool group_siblings =
true);
2328 template <
int dim,
int spacedim>
2339 template <
int dim,
int spacedim>
2340 std::vector<types::subdomain_id>
2342 const std::vector<CellId> & cell_ids);
2354 template <
int dim,
int spacedim>
2357 std::vector<types::subdomain_id> & subdomain);
2373 template <
int dim,
int spacedim>
2408 template <
int dim,
int spacedim>
2451 template <
typename MeshType>
2452 std::list<std::pair<
typename MeshType::cell_iterator,
2453 typename MeshType::cell_iterator>>
2465 template <
int dim,
int spacedim>
2479 template <
typename MeshType>
2504 template <
int dim,
int spacedim>
2560 template <
class MeshType>
2561 std::vector<typename MeshType::active_cell_iterator>
2586 template <
class Container>
2587 std::vector<typename Container::cell_iterator>
2589 const std::vector<typename Container::active_cell_iterator> &patch_cells);
2657 template <
class Container>
2660 const std::vector<typename Container::active_cell_iterator> &patch,
2662 &local_triangulation,
2665 Container::space_dimension>::active_cell_iterator,
2666 typename Container::active_cell_iterator> &patch_to_global_tria_map);
2699 template <
int dim,
int spacedim>
2702 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2718 template <
typename CellIterator>
2818 template <
typename FaceIterator>
2821 std::bitset<3> & orientation,
2822 const FaceIterator & face1,
2823 const FaceIterator & face2,
2824 const unsigned int direction,
2833 template <
typename FaceIterator>
2836 const FaceIterator & face1,
2837 const FaceIterator & face2,
2838 const unsigned int direction,
2900 template <
typename MeshType>
2903 const MeshType & mesh,
2906 const unsigned int direction,
2936 template <
typename MeshType>
2939 const MeshType & mesh,
2941 const unsigned int direction,
2944 const ::Tensor<1, MeshType::space_dimension> &offset =
2974 template <
int dim,
int spacedim>
2977 const bool reset_boundary_ids =
false);
3000 template <
int dim,
int spacedim>
3003 const std::vector<types::boundary_id> &src_boundary_ids,
3004 const std::vector<types::manifold_id> &dst_manifold_ids,
3006 const std::vector<types::boundary_id> &reset_boundary_ids = {});
3037 template <
int dim,
int spacedim>
3040 const bool compute_face_ids =
false);
3066 template <
int dim,
int spacedim>
3071 const std::set<types::manifold_id> &)> &disambiguation_function =
3072 [](
const std::set<types::manifold_id> &manifold_ids) {
3073 if (manifold_ids.size() == 1)
3074 return *manifold_ids.
begin();
3078 bool overwrite_only_flat_manifold_ids =
true);
3165 template <
typename DataType,
typename MeshType>
3168 const MeshType & mesh,
3169 const std::function<std_cxx17::optional<DataType>(
3170 const typename MeshType::active_cell_iterator &)> &
pack,
3171 const std::function<
void(
const typename MeshType::active_cell_iterator &,
3172 const DataType &)> &
unpack,
3173 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
3187 template <
typename DataType,
typename MeshType>
3190 const MeshType & mesh,
3191 const std::function<std_cxx17::optional<DataType>(
3192 const typename MeshType::level_cell_iterator &)> &
pack,
3193 const std::function<
void(
const typename MeshType::level_cell_iterator &,
3194 const DataType &)> &
unpack,
3195 const std::function<
bool(
const typename MeshType::level_cell_iterator &)> &
3210 template <
int spacedim>
3211 std::vector<std::vector<BoundingBox<spacedim>>>
3214 const MPI_Comm & mpi_communicator);
3248 template <
int spacedim>
3252 const MPI_Comm & mpi_communicator);
3271 template <
int dim,
int spacedim>
3284 template <
int dim,
int spacedim>
3285 std::map<unsigned int, std::set<::types::subdomain_id>>
3304 template <
int dim,
typename T>
3325 template <
class Archive>
3327 save(Archive &ar,
const unsigned int)
const
3331 Assert(cell_ids.size() == data.size(),
3334 const std::size_t
n_cells = cell_ids.size();
3336 for (
const auto &
id : cell_ids)
3339 ar & binary_cell_id;
3351 template <
class Archive>
3353 load(Archive &ar,
const unsigned int)
3361 for (
unsigned int c = 0; c <
n_cells; ++c)
3365 cell_ids.emplace_back(value);
3377 template <
class Archive>
3383 BOOST_SERIALIZATION_SPLIT_MEMBER()
3406 template <
int dim,
typename VectorType>
3429 const VectorType & ls_vector,
3430 const double iso_level,
3442 const VectorType & ls_vector,
3443 const double iso_level,
3461 const double iso_level,
3473 const std::vector<
Point<2>> & points,
3474 const std::vector<unsigned int> mask,
3475 const double iso_level,
3484 const std::vector<
Point<3>> & points,
3485 const std::vector<unsigned int> mask,
3486 const double iso_level,
3521 <<
"The number of partitions you gave is " << arg1
3522 <<
", but must be greater than zero.");
3528 <<
"The subdomain id " << arg1
3529 <<
" has no cells associated with it.");
3540 <<
"The scaling factor must be positive, but it is " << arg1
3548 <<
"The given vertex with index " << arg1
3549 <<
" is not used in the given triangulation.");
3566 "The edges of the mesh are not consistently orientable.");
3603 template <
int dim,
typename Predicate,
int spacedim>
3608 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
3619 for (; cell != endc; ++cell)
3620 for (
const unsigned int v : cell->vertex_indices())
3621 if (treated_vertices[cell->vertex_index(v)] ==
false)
3624 cell->vertex(v) = predicate(cell->vertex(v));
3626 treated_vertices[cell->vertex_index(v)] =
true;
3636 for (; cell != endc; ++cell)
3637 for (
const unsigned int face : cell->face_indices())
3638 if (cell->face(face)->has_children() &&
3639 !cell->face(face)->at_boundary())
3641 Assert(cell->reference_cell() ==
3642 ReferenceCells::get_hypercube<dim>(),
3646 cell->face(face)->child(0)->vertex(1) =
3647 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3656 for (; cell != endc; ++cell)
3657 for (
const unsigned int face : cell->face_indices())
3658 if (cell->face(face)->has_children() &&
3659 !cell->face(face)->at_boundary())
3661 Assert(cell->reference_cell() ==
3662 ReferenceCells::get_hypercube<dim>(),
3666 cell->face(face)->child(0)->vertex(1) =
3667 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3669 cell->face(face)->child(0)->vertex(2) =
3670 (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3672 cell->face(face)->child(1)->vertex(3) =
3673 (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3675 cell->face(face)->child(2)->vertex(3) =
3676 (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3680 cell->face(face)->child(0)->vertex(3) =
3681 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3682 cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3693 template <
class MeshType>
3694 std::vector<typename MeshType::active_cell_iterator>
3697 std::vector<typename MeshType::active_cell_iterator> child_cells;
3699 if (cell->has_children())
3701 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3702 if (cell->child(child)->has_children())
3704 const std::vector<typename MeshType::active_cell_iterator>
3705 children = get_active_child_cells<MeshType>(cell->child(child));
3706 child_cells.insert(child_cells.end(),
3711 child_cells.push_back(cell->child(child));
3719 template <
class MeshType>
3722 const typename MeshType::active_cell_iterator & cell,
3723 std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3725 active_neighbors.clear();
3726 for (
const unsigned int n : cell->face_indices())
3727 if (!cell->at_boundary(n))
3729 if (MeshType::dimension == 1)
3736 typename MeshType::cell_iterator neighbor_child =
3738 if (!neighbor_child->is_active())
3740 while (neighbor_child->has_children())
3741 neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3743 Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3746 active_neighbors.push_back(neighbor_child);
3750 if (cell->face(n)->has_children())
3754 for (
unsigned int c = 0;
3755 c < cell->face(n)->n_active_descendants();
3757 active_neighbors.push_back(
3758 cell->neighbor_child_on_subface(n, c));
3764 active_neighbors.push_back(cell->neighbor(n));
3774 namespace ProjectToObject
3788 struct CrossDerivative
3790 const unsigned int direction_0;
3791 const unsigned int direction_1;
3793 CrossDerivative(
const unsigned int d0,
const unsigned int d1);
3796 inline CrossDerivative::CrossDerivative(
const unsigned int d0,
3797 const unsigned int d1)
3808 template <
typename F>
3810 centered_first_difference(
const double center,
3814 return (f(
center + step) - f(
center - step)) / (2.0 * step);
3823 template <
typename F>
3825 centered_second_difference(
const double center,
3844 template <
int structdim,
typename F>
3847 const CrossDerivative cross_derivative,
3853 simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3854 simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3855 return (-4.0 * f(
center) - 1.0 * f(
center + simplex_vector) -
3856 1.0 / 3.0 * f(
center - simplex_vector) +
3857 16.0 / 3.0 * f(
center + 0.5 * simplex_vector)) /
3869 template <
int spacedim,
int structdim,
typename F>
3872 const unsigned int row_n,
3873 const unsigned int dependent_direction,
3880 dependent_direction <
3882 ExcMessage(
"This function assumes that the last weight is a "
3883 "dependent variable (and hence we cannot take its "
3884 "derivative directly)."));
3885 Assert(row_n != dependent_direction,
3887 "We cannot differentiate with respect to the variable "
3888 "that is assumed to be dependent."));
3892 {row_n, dependent_direction},
center, step, f);
3894 for (
unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
3896 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
3905 template <
typename Iterator,
int spacedim,
int structdim>
3907 project_to_d_linear_object(
const Iterator &
object,
3942 for (
unsigned int d = 0;
d < structdim; ++
d)
3947 x_k +=
object->vertex(i) *
3953 for (
const unsigned int i :
3956 (x_k - trial_point) * object->vertex(i) *
3961 for (
const unsigned int i :
3963 for (
const unsigned int j :
3971 H_k += (
object->vertex(i) *
object->vertex(j)) * tmp;
3978 for (
const unsigned int i :
3980 x_k +=
object->vertex(i) *
3983 if (delta_xi.
norm() < 1
e-7)
3999 template <
int structdim>
4006 static const std::size_t n_vertices_per_cell =
4008 n_independent_components;
4009 std::array<double, n_vertices_per_cell> copied_weights;
4010 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
4012 copied_weights[i] = v[i];
4013 if (v[i] < 0.0 || v[i] > 1.0)
4018 std::sort(copied_weights.begin(), copied_weights.end());
4020 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4025 template <
typename Iterator>
4028 const Iterator &
object,
4031 const int spacedim = Iterator::AccessorType::space_dimension;
4032 const int structdim = Iterator::AccessorType::structure_dimension;
4036 if (structdim >= spacedim)
4037 return projected_point;
4038 else if (structdim == 1 || structdim == 2)
4040 using namespace internal::ProjectToObject;
4045 const int dim = Iterator::AccessorType::dimension;
4048 &manifold) !=
nullptr)
4051 project_to_d_linear_object<Iterator, spacedim, structdim>(
4052 object, trial_point);
4097 const double step_size =
object->diameter() / 64.0;
4099 constexpr
unsigned int n_vertices_per_cell =
4102 std::array<Point<spacedim>, n_vertices_per_cell>
vertices;
4103 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4105 vertices[vertex_n] = object->vertex(vertex_n);
4107 auto get_point_from_weights =
4110 return object->get_manifold().get_new_point(
4118 double guess_weights_sum = 0.0;
4119 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4122 const double distance =
4124 if (distance == 0.0)
4126 guess_weights = 0.0;
4127 guess_weights[vertex_n] = 1.0;
4128 guess_weights_sum = 1.0;
4133 guess_weights[vertex_n] = 1.0 / distance;
4134 guess_weights_sum += guess_weights[vertex_n];
4137 guess_weights /= guess_weights_sum;
4138 Assert(internal::weights_are_ok<structdim>(guess_weights),
4147 for (
unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4149 const unsigned int dependent_direction =
4150 n_vertices_per_cell - 1;
4152 for (
unsigned int row_n = 0; row_n < n_vertices_per_cell;
4155 if (row_n != dependent_direction)
4157 current_gradient[row_n] =
4158 gradient_entry<spacedim, structdim>(
4160 dependent_direction,
4164 get_point_from_weights);
4166 current_gradient[dependent_direction] -=
4167 current_gradient[row_n];
4185 double gradient_weight = -0.5;
4186 auto gradient_weight_objective_function =
4187 [&](
const double gradient_weight_guess) ->
double {
4188 return (trial_point -
4189 get_point_from_weights(guess_weights +
4190 gradient_weight_guess *
4195 for (
unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4197 const double update_numerator = centered_first_difference(
4200 gradient_weight_objective_function);
4201 const double update_denominator =
4202 centered_second_difference(
4205 gradient_weight_objective_function);
4209 if (std::abs(update_denominator) == 0.0)
4212 gradient_weight - update_numerator / update_denominator;
4217 if (std::abs(gradient_weight) > 10)
4219 gradient_weight = -10.0;
4229 guess_weights + gradient_weight * current_gradient;
4231 double new_gradient_weight = gradient_weight;
4232 for (
unsigned int iteration_count = 0; iteration_count < 40;
4235 if (internal::weights_are_ok<structdim>(tentative_weights))
4238 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
4240 if (tentative_weights[i] < 0.0)
4242 tentative_weights -=
4243 (tentative_weights[i] / current_gradient[i]) *
4246 if (tentative_weights[i] < 0.0 ||
4247 1.0 < tentative_weights[i])
4249 new_gradient_weight /= 2.0;
4252 new_gradient_weight * current_gradient;
4259 if (!internal::weights_are_ok<structdim>(tentative_weights))
4264 if (get_point_from_weights(tentative_weights)
4265 .distance(trial_point) <
4266 get_point_from_weights(guess_weights).distance(trial_point))
4267 guess_weights = tentative_weights;
4270 Assert(internal::weights_are_ok<structdim>(guess_weights),
4273 Assert(internal::weights_are_ok<structdim>(guess_weights),
4275 projected_point = get_point_from_weights(guess_weights);
4285 for (
unsigned int line_n = 0;
4286 line_n < GeometryInfo<structdim>::lines_per_cell;
4289 line_projections[line_n] =
4292 std::sort(line_projections.begin(),
4293 line_projections.end(),
4295 return a.distance(trial_point) <
4296 b.distance(trial_point);
4298 if (line_projections[0].distance(trial_point) <
4299 projected_point.
distance(trial_point))
4300 projected_point = line_projections[0];
4306 return projected_point;
4309 return projected_point;
4316 template <
typename DataType,
4318 typename MeshCellIteratorType>
4321 const MeshType &mesh,
4322 const std::function<
4323 std_cxx17::optional<DataType>(
const MeshCellIteratorType &)> &
pack,
4324 const std::function<
void(
const MeshCellIteratorType &,
const DataType &)>
4326 const std::function<
bool(
const MeshCellIteratorType &)> & cell_filter,
4327 const std::function<
void(
4328 const std::function<
void(
const MeshCellIteratorType &,
4330 const std::function<std::set<types::subdomain_id>(
4332 MeshType::space_dimension> &)>
4333 &compute_ghost_owners)
4335 # ifndef DEAL_II_WITH_MPI
4340 (void)process_cells;
4341 (void)compute_ghost_owners;
4344 constexpr
int dim = MeshType::dimension;
4345 constexpr
int spacedim = MeshType::space_dimension;
4348 &mesh.get_triangulation());
4352 "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4354 if (
const auto tria =
dynamic_cast<
4356 &mesh.get_triangulation()))
4359 tria->with_artificial_cells(),
4361 "The functions GridTools::exchange_cell_data_to_ghosts() and "
4362 "GridTools::exchange_cell_data_to_level_ghosts() can only "
4363 "operate on a single layer of ghost cells. However, you have "
4364 "given a Triangulation object of type "
4365 "parallel::shared::Triangulation without artificial cells "
4366 "resulting in arbitrary numbers of ghost layers."));
4370 std::set<::types::subdomain_id> ghost_owners =
4371 compute_ghost_owners(*
tria);
4373 std::vector<typename CellId::binary_type>>
4376 for (
const auto ghost_owner : ghost_owners)
4377 neighbor_cell_list[ghost_owner] = {};
4379 process_cells([&](
const auto &cell,
const auto key) ->
void {
4380 if (cell_filter(cell))
4382 constexpr
int spacedim = MeshType::space_dimension;
4383 neighbor_cell_list[key].emplace_back(
4384 cell->id().template to_binary<spacedim>());
4388 Assert(ghost_owners.size() == neighbor_cell_list.size(),
4400 const int mpi_tag_reply =
4404 std::vector<MPI_Request> requests(ghost_owners.size());
4406 unsigned int idx = 0;
4407 for (
const auto &it : neighbor_cell_list)
4410 const int ierr = MPI_Isend(it.second.data(),
4411 it.second.size() *
sizeof(it.second[0]),
4423 std::vector<MPI_Request> reply_requests(ghost_owners.size());
4424 std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4426 for (
unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4429 int ierr = MPI_Probe(MPI_ANY_SOURCE,
4436 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4443 std::vector<typename CellId::binary_type> cells_with_requests(
4445 std::vector<DataType> data_to_send;
4446 data_to_send.reserve(
n_cells);
4447 std::vector<bool> cell_carries_data(
n_cells,
false);
4449 ierr = MPI_Recv(cells_with_requests.data(),
4459 for (
unsigned int c = 0; c < static_cast<unsigned int>(
n_cells); ++c)
4464 MeshCellIteratorType mesh_it(
tria,
4469 const std_cxx17::optional<DataType> data =
pack(mesh_it);
4472 data_to_send.emplace_back(*data);
4473 cell_carries_data[c] =
true;
4481 sendbuffers[idx].resize(
sizeof(std::size_t));
4486 std::size_t size_of_send =
4490 std::memcpy(sendbuffers[idx].data(),
4492 sizeof(std::size_t));
4496 if (data_to_send.size() <
n_cells)
4502 ierr = MPI_Isend(sendbuffers[idx].data(),
4503 sendbuffers[idx].size(),
4508 &reply_requests[idx]);
4513 std::vector<char> receive;
4514 for (
unsigned int id = 0;
id < neighbor_cell_list.size(); ++id)
4517 int ierr = MPI_Probe(MPI_ANY_SOURCE,
4524 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4527 receive.resize(len);
4529 ierr = MPI_Recv(receive.data(),
4540 auto data_iterator = receive.begin();
4541 std::size_t size_of_received_data =
4542 Utilities::unpack<std::size_t>(data_iterator,
4543 data_iterator +
sizeof(std::size_t));
4544 data_iterator +=
sizeof(std::size_t);
4547 auto received_data = Utilities::unpack<std::vector<DataType>>(
4549 data_iterator + size_of_received_data,
4551 data_iterator += size_of_received_data;
4556 const std::vector<typename CellId::binary_type> &this_cell_list =
4557 neighbor_cell_list[status.MPI_SOURCE];
4559 std::vector<bool> cells_with_data;
4560 if (received_data.size() < this_cell_list.size())
4562 cells_with_data = Utilities::unpack<std::vector<bool>>(
4563 data_iterator, receive.end(),
false);
4569 auto received_data_iterator = received_data.begin();
4570 for (
unsigned int c = 0; c < this_cell_list.size(); ++c)
4571 if (cells_with_data.empty() || cells_with_data[c])
4576 MeshCellIteratorType cell(
tria,
4581 unpack(cell, *received_data_iterator);
4582 ++received_data_iterator;
4588 if (requests.size() > 0)
4591 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4594 if (reply_requests.size() > 0)
4596 const int ierr = MPI_Waitall(reply_requests.size(),
4597 reply_requests.data(),
4598 MPI_STATUSES_IGNORE);
4608 template <
typename DataType,
typename MeshType>
4611 const MeshType & mesh,
4612 const std::function<std_cxx17::optional<DataType>(
4613 const typename MeshType::active_cell_iterator &)> &
pack,
4614 const std::function<
void(
const typename MeshType::active_cell_iterator &,
4615 const DataType &)> &
unpack,
4616 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
4619 # ifndef DEAL_II_WITH_MPI
4626 internal::exchange_cell_data<DataType,
4628 typename MeshType::active_cell_iterator>(
4633 [&](
const auto &process) {
4634 for (
const auto &cell : mesh.active_cell_iterators())
4635 if (cell->is_ghost())
4636 process(cell, cell->subdomain_id());
4638 [](
const auto &
tria) {
return tria.ghost_owners(); });
4644 template <
typename DataType,
typename MeshType>
4647 const MeshType & mesh,
4648 const std::function<std_cxx17::optional<DataType>(
4649 const typename MeshType::level_cell_iterator &)> &
pack,
4650 const std::function<
void(
const typename MeshType::level_cell_iterator &,
4651 const DataType &)> &
unpack,
4652 const std::function<
bool(
const typename MeshType::level_cell_iterator &)>
4655 # ifndef DEAL_II_WITH_MPI
4662 internal::exchange_cell_data<DataType,
4664 typename MeshType::level_cell_iterator>(
4669 [&](
const auto &process) {
4670 for (
const auto &cell : mesh.cell_iterators())
4671 if (cell->is_ghost_on_level())
4672 process(cell, cell->level_subdomain_id());
4674 [](
const auto &
tria) {
return tria.level_ghost_owners(); });
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
std::array< unsigned int, 4 > binary_type
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
numbers::NumberTraits< Number >::real_type norm() const
virtual MPI_Comm get_communicator() const
cell_iterator create_cell_iterator(const CellId &cell_id) const
cell_iterator begin(const unsigned int level=0) const
typename MeshType::active_cell_iterator type
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
static ::ExceptionBase & ExcMeshNotOrientable()
static ::ExceptionBase & ExcScalingFactorNotPositive(double arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsMPI()
static ::ExceptionBase & ExcVertexNotUsed(unsigned int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcTriangulationHasBeenRefined()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNonExistentSubdomain(int arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void collect_coinciding_vertices(const Triangulation< dim, spacedim > &tria, std::map< unsigned int, std::vector< unsigned int >> &coinciding_vertex_groups, std::map< unsigned int, unsigned int > &vertex_to_coinciding_vertex_group)
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
void copy_material_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool compute_face_ids=false)
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::map< unsigned int, std::set<::types::subdomain_id > > compute_vertices_with_ghost_neighbors(const Triangulation< dim, spacedim > &tria)
void exchange_cell_data_to_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::active_cell_iterator &)> &pack, const std::function< void(const typename MeshType::active_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::active_cell_iterator &)> &cell_filter=always_return< typename MeshType::active_cell_iterator, bool >{true})
std::vector< std::vector< BoundingBox< spacedim > > > exchange_local_bounding_boxes(const std::vector< BoundingBox< spacedim >> &local_bboxes, const MPI_Comm &mpi_communicator)
void exchange_cell_data_to_level_ghosts(const MeshType &mesh, const std::function< std_cxx17::optional< DataType >(const typename MeshType::level_cell_iterator &)> &pack, const std::function< void(const typename MeshType::level_cell_iterator &, const DataType &)> &unpack, const std::function< bool(const typename MeshType::level_cell_iterator &)> &cell_filter=always_return< typename MeshType::level_cell_iterator, bool >{ true})
RTree< std::pair< BoundingBox< spacedim >, unsigned int > > build_global_description_tree(const std::vector< BoundingBox< spacedim >> &local_description, const MPI_Comm &mpi_communicator)
void assign_co_dimensional_manifold_indicators(Triangulation< dim, spacedim > &tria, const std::function< types::manifold_id(const std::set< types::manifold_id > &)> &disambiguation_function=[](const std::set< types::manifold_id > &manifold_ids) { if(manifold_ids.size()==1) return *manifold_ids.begin();else return numbers::flat_manifold_id;}, bool overwrite_only_flat_manifold_ids=true)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
T sum(const T &t, const MPI_Comm &mpi_communicator)
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
T unpack(const std::vector< char > &buffer, const bool allow_compression=true)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
unsigned int global_dof_index
unsigned int subdomain_id
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
boost::geometry::index::rtree< LeafType, IndexType, IndexableGetter > RTree
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
std::map< unsigned int, unsigned int > vertex_to_coinciding_vertex_group
std::map< unsigned int, std::vector< unsigned int > > coinciding_vertex_groups
const ::Triangulation< dim, spacedim > & tria