16 #ifndef dealii_grid_tools_h
17 #define dealii_grid_tools_h
49 #include <boost/archive/binary_iarchive.hpp>
50 #include <boost/archive/binary_oarchive.hpp>
51 #include <boost/random/mersenne_twister.hpp>
52 #include <boost/serialization/array.hpp>
53 #include <boost/serialization/vector.hpp>
55 #ifdef DEAL_II_WITH_ZLIB
56 # include <boost/iostreams/device/back_inserter.hpp>
57 # include <boost/iostreams/filter/gzip.hpp>
58 # include <boost/iostreams/filtering_stream.hpp>
59 # include <boost/iostreams/stream.hpp>
66 #ifdef DEAL_II_HAVE_CXX20
79 template <
int dim,
int spacedim>
87 class MappingCollection;
94 template <
int dim,
int spacedim>
101 template <
int dim,
int spacedim,
class MeshType>
106 using type =
typename MeshType::active_cell_iterator;
113 template <
int dim,
int spacedim>
144 template <
int dim,
int spacedim>
171 template <
int dim,
int spacedim>
202 template <
int dim,
int 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>()
240 template <
int dim,
int spacedim>
245 (ReferenceCells::get_hypercube<dim>()
247 .
template get_default_linear_mapping<dim, spacedim>()
249 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
316 template <
int dim,
int spacedim>
378 template <
int dim,
int spacedim>
399 template <
typename Iterator>
402 const Iterator &
object,
417 template <
int dim,
int spacedim>
419 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>,
SubCellData>
444 template <
int dim,
int spacedim>
468 template <
int dim,
int spacedim>
473 std::vector<unsigned int> & considered_vertices,
474 const double tol = 1
e-12);
486 const double tol = 1
e-12);
506 template <
int dim,
int spacedim>
521 template <
int dim,
int spacedim>
609 template <
int dim,
typename Transformation,
int spacedim>
612 std::assignable_from<
624 template <
int dim,
int spacedim>
640 template <
int dim,
int spacedim>
680 const
unsigned int axis,
744 const
Function<dim,
double> *coefficient =
nullptr,
745 const
bool solve_for_absolute_positions = false);
752 template <
int dim,
int spacedim>
753 std::map<
unsigned int,
Point<spacedim>>
763 template <
int dim,
int spacedim>
765 scale(const
double scaling_factor,
782 template <
int dim,
int spacedim>
787 const
bool keep_boundary = true,
788 const
unsigned int seed =
boost::
random::mt19937::default_seed);
823 template <
int dim,
int spacedim>
826 const
bool isotropic = false,
827 const
unsigned int max_iterations = 100);
853 template <
int dim,
int spacedim>
856 const
double max_ratio = 1.6180339887,
857 const
unsigned int max_iterations = 5);
948 template <
int dim,
int spacedim>
951 const
double limit_angle_fraction = .75);
1013 template <
int dim,
int spacedim>
1016 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1017 std::vector<std::vector<Point<dim>>>,
1018 std::vector<std::vector<unsigned int>>>
1062 template <
int dim,
int spacedim>
1065 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1066 std::vector<std::vector<Point<dim>>>,
1067 std::vector<std::vector<unsigned int>>,
1068 std::vector<unsigned int>>
1148 template <
int dim,
int spacedim>
1151 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1152 std::vector<std::vector<Point<dim>>>,
1153 std::vector<std::vector<unsigned int>>,
1154 std::vector<std::vector<Point<spacedim>>>,
1155 std::vector<std::vector<unsigned int>>>
1163 const double tolerance = 1
e-10);
1181 template <
int dim,
int spacedim>
1190 std::vector<std::tuple<std::pair<int, int>,
1220 std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1245 template <
int dim,
int spacedim>
1251 const std::vector<bool> & marked_vertices,
1252 const double tolerance,
1253 const bool perform_handshake,
1254 const bool enforce_unique_mapping =
false);
1294 template <
int dim,
int spacedim>
1295 std::map<unsigned int, Point<spacedim>>
1299 (ReferenceCells::get_hypercube<dim>()
1301 .
template get_default_linear_mapping<dim, spacedim>()
1303 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1316 template <
int spacedim>
1344 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1348 const std::vector<bool> & marked_vertices = {});
1373 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1376 const MeshType<dim, spacedim> &mesh,
1378 const std::vector<bool> & marked_vertices = {});
1400 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1402 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1405 typename ::internal::
1406 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1409 const unsigned int vertex_index);
1473 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1475 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1477 std::pair<typename ::internal::
1478 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1482 const MeshType<dim, spacedim> &mesh,
1484 const std::vector<bool> &marked_vertices = {},
1485 const double tolerance = 1.e-10);
1494 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1496 typename MeshType<dim, spacedim>::active_cell_iterator
1498 typename ::internal::
1499 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1503 const std::vector<bool> &marked_vertices = {},
1504 const double tolerance = 1.e-10);
1512 template <
int dim,
int spacedim>
1513 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1519 const double tolerance = 1.e-10);
1572 template <
int dim,
int spacedim>
1573 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1576 const Cache<dim, spacedim> &cache,
1580 const std::vector<bool> &marked_vertices = {},
1581 const double tolerance = 1.e-10);
1596 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1598 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1600 std::pair<typename ::internal::
1601 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1606 const MeshType<dim, spacedim> &mesh,
1609 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
1612 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1613 typename MeshType<dim, spacedim>::active_cell_iterator(),
1614 const std::vector<bool> & marked_vertices = {},
1617 const double tolerance = 1.e-10,
1619 std::pair<BoundingBox<spacedim>,
1621 *relevant_cell_bounding_boxes_rtree =
nullptr);
1644 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1646 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1649 std::vector<std::pair<
1650 typename ::internal::
1651 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1656 const MeshType<dim, spacedim> &mesh,
1658 const double tolerance,
1659 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
1668 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1670 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1673 std::vector<std::pair<
1674 typename ::internal::
1675 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1680 const MeshType<dim, spacedim> &mesh,
1682 const double tolerance = 1
e-10,
1683 const std::vector<bool> & marked_vertices = {});
1706 template <
class MeshType>
1707 std::vector<typename MeshType::active_cell_iterator>
1734 template <
class MeshType>
1737 const typename MeshType::active_cell_iterator & cell,
1738 std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1789 template <
class MeshType>
1790 std::vector<typename MeshType::active_cell_iterator>
1792 const MeshType &mesh,
1793 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1804 template <
class MeshType>
1805 std::vector<typename MeshType::cell_iterator>
1807 const MeshType &mesh,
1808 const std::function<
bool(
const typename MeshType::cell_iterator &)>
1810 const unsigned int level);
1825 template <
class MeshType>
1826 std::vector<typename MeshType::active_cell_iterator>
1877 template <
class MeshType>
1878 std::vector<typename MeshType::active_cell_iterator>
1880 const MeshType &mesh,
1881 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
1883 const double layer_thickness);
1907 template <
class MeshType>
1908 std::vector<typename MeshType::active_cell_iterator>
1910 const double layer_thickness);
1927 template <
class MeshType>
1930 const MeshType &mesh,
1931 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
2004 template <
class MeshType>
2005 std::vector<BoundingBox<MeshType::space_dimension>>
2007 const MeshType &mesh,
2008 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
2010 const unsigned int refinement_level = 0,
2011 const bool allow_merge =
false,
2041 template <
int spacedim>
2043 std::tuple<std::vector<std::vector<unsigned int>>,
2044 std::map<unsigned int, unsigned int>,
2045 std::map<unsigned int, std::vector<unsigned int>>>
2088 template <
int spacedim>
2090 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2091 std::map<unsigned int, unsigned int>,
2092 std::map<unsigned int, std::vector<unsigned int>>>
2109 template <
int dim,
int spacedim>
2111 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2126 template <
int dim,
int spacedim>
2127 std::vector<std::vector<Tensor<1, spacedim>>>
2142 template <
int dim,
int spacedim>
2148 (ReferenceCells::get_hypercube<dim>()
2150 .
template get_default_linear_mapping<dim, spacedim>()
2152 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2167 template <
int dim,
int spacedim>
2168 std::map<unsigned int, types::global_vertex_index>
2183 template <
int dim,
int spacedim>
2184 std::pair<unsigned int, double>
2202 template <
int dim,
int spacedim>
2216 template <
int dim,
int spacedim>
2230 template <
int dim,
int spacedim>
2234 const unsigned int level,
2257 template <
int dim,
int spacedim>
2274 template <
int dim,
int spacedim>
2277 const std::vector<unsigned int> &cell_weights,
2327 template <
int dim,
int spacedim>
2345 template <
int dim,
int spacedim>
2348 const std::vector<unsigned int> &cell_weights,
2368 template <
int dim,
int spacedim>
2372 const bool group_siblings =
true);
2385 template <
int dim,
int spacedim>
2396 template <
int dim,
int spacedim>
2397 std::vector<types::subdomain_id>
2399 const std::vector<CellId> & cell_ids);
2411 template <
int dim,
int spacedim>
2414 std::vector<types::subdomain_id> & subdomain);
2430 template <
int dim,
int spacedim>
2465 template <
int dim,
int spacedim>
2508 template <
typename MeshType>
2509 std::list<std::pair<
typename MeshType::cell_iterator,
2510 typename MeshType::cell_iterator>>
2522 template <
int dim,
int spacedim>
2536 template <
typename MeshType>
2561 template <
int dim,
int spacedim>
2617 template <
class MeshType>
2618 std::vector<typename MeshType::active_cell_iterator>
2643 template <
class Container>
2644 std::vector<typename Container::cell_iterator>
2646 const std::vector<typename Container::active_cell_iterator> &patch_cells);
2714 template <
class Container>
2717 const std::vector<typename Container::active_cell_iterator> &patch,
2719 &local_triangulation,
2722 Container::space_dimension>::active_cell_iterator,
2723 typename Container::active_cell_iterator> &patch_to_global_tria_map);
2756 template <
int dim,
int spacedim>
2759 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2775 template <
typename CellIterator>
2881 template <
typename FaceIterator>
2884 std::bitset<3> & orientation,
2885 const FaceIterator & face1,
2886 const FaceIterator & face2,
2887 const unsigned int direction,
2896 template <
typename FaceIterator>
2899 const FaceIterator & face1,
2900 const FaceIterator & face2,
2901 const unsigned int direction,
2963 template <
typename MeshType>
2966 const MeshType & mesh,
2969 const unsigned int direction,
2999 template <
typename MeshType>
3002 const MeshType & mesh,
3004 const unsigned int direction,
3007 const ::Tensor<1, MeshType::space_dimension> &offset =
3037 template <
int dim,
int spacedim>
3040 const bool reset_boundary_ids =
false);
3063 template <
int dim,
int spacedim>
3066 const std::vector<types::boundary_id> &src_boundary_ids,
3067 const std::vector<types::manifold_id> &dst_manifold_ids,
3069 const std::vector<types::boundary_id> &reset_boundary_ids = {});
3100 template <
int dim,
int spacedim>
3103 const bool compute_face_ids =
false);
3129 template <
int dim,
int spacedim>
3134 const std::set<types::manifold_id> &)> &disambiguation_function =
3135 [](
const std::set<types::manifold_id> &manifold_ids) {
3136 if (manifold_ids.size() == 1)
3137 return *manifold_ids.
begin();
3141 bool overwrite_only_flat_manifold_ids =
true);
3228 template <
typename DataType,
typename MeshType>
3231 const MeshType & mesh,
3232 const std::function<std_cxx17::optional<DataType>(
3233 const typename MeshType::active_cell_iterator &)> &
pack,
3234 const std::function<
void(
const typename MeshType::active_cell_iterator &,
3235 const DataType &)> &
unpack,
3236 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
3250 template <
typename DataType,
typename MeshType>
3253 const MeshType & mesh,
3254 const std::function<std_cxx17::optional<DataType>(
3255 const typename MeshType::level_cell_iterator &)> &
pack,
3256 const std::function<
void(
const typename MeshType::level_cell_iterator &,
3257 const DataType &)> &
unpack,
3258 const std::function<
bool(
const typename MeshType::level_cell_iterator &)> &
3273 template <
int spacedim>
3274 std::vector<std::vector<BoundingBox<spacedim>>>
3277 const MPI_Comm & mpi_communicator);
3311 template <
int spacedim>
3315 const MPI_Comm & mpi_communicator);
3334 template <
int dim,
int spacedim>
3347 template <
int dim,
int spacedim>
3348 std::map<unsigned int, std::set<::types::subdomain_id>>
3367 template <
int dim,
typename T>
3388 template <
class Archive>
3390 save(Archive &ar,
const unsigned int)
const
3394 Assert(cell_ids.size() == data.size(),
3397 const std::size_t
n_cells = cell_ids.size();
3399 for (
const auto &
id : cell_ids)
3402 ar & binary_cell_id;
3414 template <
class Archive>
3416 load(Archive &ar,
const unsigned int)
3424 for (
unsigned int c = 0; c <
n_cells; ++c)
3428 cell_ids.emplace_back(value);
3440 template <
class Archive>
3446 BOOST_SERIALIZATION_SPLIT_MEMBER()
3472 template <
int dim,
typename VectorType>
3501 const VectorType & ls_vector,
3502 const double iso_level,
3512 const VectorType & ls_vector,
3513 const double iso_level,
3527 const VectorType & ls_vector,
3528 const double iso_level,
3538 const VectorType & ls_vector,
3539 const double iso_level,
3556 const double iso_level,
3559 const bool write_back_cell_data =
true)
const;
3567 const std::vector<unsigned int> &,
3584 const std::vector<
Point<2>> & points,
3585 const std::vector<unsigned int> &mask,
3586 const double iso_level,
3589 const bool write_back_cell_data)
const;
3596 const std::vector<
Point<3>> & points,
3597 const std::vector<unsigned int> &mask,
3598 const double iso_level,
3601 const bool write_back_cell_data)
const;
3634 <<
"The number of partitions you gave is " << arg1
3635 <<
", but must be greater than zero.");
3641 <<
"The subdomain id " << arg1
3642 <<
" has no cells associated with it.");
3653 <<
"The scaling factor must be positive, but it is " << arg1
3661 <<
"The given vertex with index " << arg1
3662 <<
" is not used in the given triangulation.");
3679 "The edges of the mesh are not consistently orientable.");
3716 template <
int dim,
typename Transformation,
int spacedim>
3719 std::assignable_from<
3722 void
transform(const Transformation & transformation,
3725 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
3736 for (; cell != endc; ++cell)
3737 for (
const unsigned int v : cell->vertex_indices())
3738 if (treated_vertices[cell->vertex_index(v)] ==
false)
3741 cell->vertex(v) = transformation(cell->vertex(v));
3743 treated_vertices[cell->vertex_index(v)] =
true;
3753 for (; cell != endc; ++cell)
3754 for (
const unsigned int face : cell->face_indices())
3755 if (cell->face(face)->has_children() &&
3756 !cell->face(face)->at_boundary())
3758 Assert(cell->reference_cell() ==
3759 ReferenceCells::get_hypercube<dim>(),
3763 cell->face(face)->child(0)->vertex(1) =
3764 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3773 for (; cell != endc; ++cell)
3774 for (
const unsigned int face : cell->face_indices())
3775 if (cell->face(face)->has_children() &&
3776 !cell->face(face)->at_boundary())
3778 Assert(cell->reference_cell() ==
3779 ReferenceCells::get_hypercube<dim>(),
3783 cell->face(face)->child(0)->vertex(1) =
3784 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3786 cell->face(face)->child(0)->vertex(2) =
3787 (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3789 cell->face(face)->child(1)->vertex(3) =
3790 (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3792 cell->face(face)->child(2)->vertex(3) =
3793 (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3797 cell->face(face)->child(0)->vertex(3) =
3798 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3799 cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3810 template <
class MeshType>
3811 std::vector<typename MeshType::active_cell_iterator>
3814 std::vector<typename MeshType::active_cell_iterator> child_cells;
3816 if (cell->has_children())
3818 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3819 if (cell->child(child)->has_children())
3821 const std::vector<typename MeshType::active_cell_iterator>
3822 children = get_active_child_cells<MeshType>(cell->child(child));
3823 child_cells.insert(child_cells.end(),
3828 child_cells.push_back(cell->child(child));
3836 template <
class MeshType>
3839 const typename MeshType::active_cell_iterator & cell,
3840 std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3842 active_neighbors.clear();
3843 for (
const unsigned int n : cell->face_indices())
3844 if (!cell->at_boundary(n))
3846 if (MeshType::dimension == 1)
3853 typename MeshType::cell_iterator neighbor_child =
3855 if (!neighbor_child->is_active())
3857 while (neighbor_child->has_children())
3858 neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3860 Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3863 active_neighbors.push_back(neighbor_child);
3867 if (cell->face(n)->has_children())
3871 for (
unsigned int c = 0;
3872 c < cell->face(n)->n_active_descendants();
3874 active_neighbors.push_back(
3875 cell->neighbor_child_on_subface(n, c));
3881 active_neighbors.push_back(cell->neighbor(n));
3889 template <
typename CellIterator>
3893 return sizeof(*this) +
matrix.memory_consumption();
3900 namespace ProjectToObject
3914 struct CrossDerivative
3916 const unsigned int direction_0;
3917 const unsigned int direction_1;
3919 CrossDerivative(
const unsigned int d0,
const unsigned int d1);
3922 inline CrossDerivative::CrossDerivative(
const unsigned int d0,
3923 const unsigned int d1)
3934 template <
typename F>
3936 centered_first_difference(
const double center,
3940 return (f(
center + step) - f(
center - step)) / (2.0 * step);
3949 template <
typename F>
3951 centered_second_difference(
const double center,
3970 template <
int structdim,
typename F>
3973 const CrossDerivative cross_derivative,
3979 simplex_vector[cross_derivative.direction_0] = 0.5 * step;
3980 simplex_vector[cross_derivative.direction_1] = -0.5 * step;
3981 return (-4.0 * f(
center) - 1.0 * f(
center + simplex_vector) -
3982 1.0 / 3.0 * f(
center - simplex_vector) +
3983 16.0 / 3.0 * f(
center + 0.5 * simplex_vector)) /
3995 template <
int spacedim,
int structdim,
typename F>
3998 const unsigned int row_n,
3999 const unsigned int dependent_direction,
4006 dependent_direction <
4008 ExcMessage(
"This function assumes that the last weight is a "
4009 "dependent variable (and hence we cannot take its "
4010 "derivative directly)."));
4011 Assert(row_n != dependent_direction,
4013 "We cannot differentiate with respect to the variable "
4014 "that is assumed to be dependent."));
4018 {row_n, dependent_direction},
center, step, f);
4020 for (
unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
4022 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
4031 template <
typename Iterator,
int spacedim,
int structdim>
4033 project_to_d_linear_object(
const Iterator &
object,
4068 for (
unsigned int d = 0;
d < structdim; ++
d)
4073 x_k +=
object->vertex(i) *
4079 for (
const unsigned int i :
4082 (x_k - trial_point) * object->vertex(i) *
4087 for (
const unsigned int i :
4089 for (
const unsigned int j :
4097 H_k += (
object->vertex(i) *
object->vertex(j)) * tmp;
4104 for (
const unsigned int i :
4106 x_k +=
object->vertex(i) *
4109 if (delta_xi.
norm() < 1
e-7)
4125 template <
int structdim>
4132 static const std::size_t n_vertices_per_cell =
4134 n_independent_components;
4135 std::array<double, n_vertices_per_cell> copied_weights;
4136 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
4138 copied_weights[i] = v[i];
4139 if (v[i] < 0.0 || v[i] > 1.0)
4144 std::sort(copied_weights.begin(), copied_weights.end());
4146 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4151 template <
typename Iterator>
4154 const Iterator &
object,
4157 const int spacedim = Iterator::AccessorType::space_dimension;
4158 const int structdim = Iterator::AccessorType::structure_dimension;
4162 if (structdim >= spacedim)
4163 return projected_point;
4164 else if (structdim == 1 || structdim == 2)
4166 using namespace internal::ProjectToObject;
4171 const int dim = Iterator::AccessorType::dimension;
4174 &manifold) !=
nullptr)
4177 project_to_d_linear_object<Iterator, spacedim, structdim>(
4178 object, trial_point);
4223 const double step_size =
object->diameter() / 64.0;
4225 constexpr
unsigned int n_vertices_per_cell =
4228 std::array<Point<spacedim>, n_vertices_per_cell>
vertices;
4229 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4231 vertices[vertex_n] = object->vertex(vertex_n);
4233 auto get_point_from_weights =
4236 return object->get_manifold().get_new_point(
4244 double guess_weights_sum = 0.0;
4245 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4248 const double distance =
4250 if (distance == 0.0)
4252 guess_weights = 0.0;
4253 guess_weights[vertex_n] = 1.0;
4254 guess_weights_sum = 1.0;
4259 guess_weights[vertex_n] = 1.0 / distance;
4260 guess_weights_sum += guess_weights[vertex_n];
4263 guess_weights /= guess_weights_sum;
4264 Assert(internal::weights_are_ok<structdim>(guess_weights),
4273 for (
unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4275 const unsigned int dependent_direction =
4276 n_vertices_per_cell - 1;
4278 for (
unsigned int row_n = 0; row_n < n_vertices_per_cell;
4281 if (row_n != dependent_direction)
4283 current_gradient[row_n] =
4284 gradient_entry<spacedim, structdim>(
4286 dependent_direction,
4290 get_point_from_weights);
4292 current_gradient[dependent_direction] -=
4293 current_gradient[row_n];
4311 double gradient_weight = -0.5;
4312 auto gradient_weight_objective_function =
4313 [&](
const double gradient_weight_guess) ->
double {
4314 return (trial_point -
4315 get_point_from_weights(guess_weights +
4316 gradient_weight_guess *
4321 for (
unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4323 const double update_numerator = centered_first_difference(
4326 gradient_weight_objective_function);
4327 const double update_denominator =
4328 centered_second_difference(
4331 gradient_weight_objective_function);
4335 if (std::abs(update_denominator) == 0.0)
4338 gradient_weight - update_numerator / update_denominator;
4343 if (std::abs(gradient_weight) > 10)
4345 gradient_weight = -10.0;
4355 guess_weights + gradient_weight * current_gradient;
4357 double new_gradient_weight = gradient_weight;
4358 for (
unsigned int iteration_count = 0; iteration_count < 40;
4361 if (internal::weights_are_ok<structdim>(tentative_weights))
4364 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
4366 if (tentative_weights[i] < 0.0)
4368 tentative_weights -=
4369 (tentative_weights[i] / current_gradient[i]) *
4372 if (tentative_weights[i] < 0.0 ||
4373 1.0 < tentative_weights[i])
4375 new_gradient_weight /= 2.0;
4378 new_gradient_weight * current_gradient;
4385 if (!internal::weights_are_ok<structdim>(tentative_weights))
4390 if (get_point_from_weights(tentative_weights)
4391 .distance(trial_point) <
4392 get_point_from_weights(guess_weights).distance(trial_point))
4393 guess_weights = tentative_weights;
4396 Assert(internal::weights_are_ok<structdim>(guess_weights),
4399 Assert(internal::weights_are_ok<structdim>(guess_weights),
4401 projected_point = get_point_from_weights(guess_weights);
4411 for (
unsigned int line_n = 0;
4412 line_n < GeometryInfo<structdim>::lines_per_cell;
4415 line_projections[line_n] =
4418 std::sort(line_projections.begin(),
4419 line_projections.end(),
4421 return a.distance(trial_point) <
4422 b.distance(trial_point);
4424 if (line_projections[0].distance(trial_point) <
4425 projected_point.
distance(trial_point))
4426 projected_point = line_projections[0];
4432 return projected_point;
4435 return projected_point;
4442 template <
typename DataType,
4444 typename MeshCellIteratorType>
4447 const MeshType &mesh,
4448 const std::function<
4449 std_cxx17::optional<DataType>(
const MeshCellIteratorType &)> &
pack,
4450 const std::function<
void(
const MeshCellIteratorType &,
const DataType &)>
4452 const std::function<
bool(
const MeshCellIteratorType &)> & cell_filter,
4453 const std::function<
void(
4454 const std::function<
void(
const MeshCellIteratorType &,
4456 const std::function<std::set<types::subdomain_id>(
4458 MeshType::space_dimension> &)>
4459 &compute_ghost_owners)
4461 # ifndef DEAL_II_WITH_MPI
4466 (void)process_cells;
4467 (void)compute_ghost_owners;
4470 constexpr
int dim = MeshType::dimension;
4471 constexpr
int spacedim = MeshType::space_dimension;
4474 &mesh.get_triangulation());
4478 "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4480 if (
const auto tria =
dynamic_cast<
4482 &mesh.get_triangulation()))
4485 tria->with_artificial_cells(),
4487 "The functions GridTools::exchange_cell_data_to_ghosts() and "
4488 "GridTools::exchange_cell_data_to_level_ghosts() can only "
4489 "operate on a single layer of ghost cells. However, you have "
4490 "given a Triangulation object of type "
4491 "parallel::shared::Triangulation without artificial cells "
4492 "resulting in arbitrary numbers of ghost layers."));
4496 std::set<::types::subdomain_id> ghost_owners =
4497 compute_ghost_owners(*
tria);
4499 std::vector<typename CellId::binary_type>>
4502 for (
const auto ghost_owner : ghost_owners)
4503 neighbor_cell_list[ghost_owner] = {};
4505 process_cells([&](
const auto &cell,
const auto key) ->
void {
4506 if (cell_filter(cell))
4508 constexpr
int spacedim = MeshType::space_dimension;
4509 neighbor_cell_list[key].emplace_back(
4510 cell->id().template to_binary<spacedim>());
4514 Assert(ghost_owners.size() == neighbor_cell_list.size(),
4526 const int mpi_tag_reply =
4530 std::vector<MPI_Request> requests(ghost_owners.size());
4532 unsigned int idx = 0;
4533 for (
const auto &it : neighbor_cell_list)
4536 const int ierr = MPI_Isend(it.second.data(),
4537 it.second.size() *
sizeof(it.second[0]),
4549 std::vector<MPI_Request> reply_requests(ghost_owners.size());
4550 std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4552 for (
unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4555 int ierr = MPI_Probe(MPI_ANY_SOURCE,
4562 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4569 std::vector<typename CellId::binary_type> cells_with_requests(
4571 std::vector<DataType> data_to_send;
4572 data_to_send.reserve(
n_cells);
4573 std::vector<bool> cell_carries_data(
n_cells,
false);
4575 ierr = MPI_Recv(cells_with_requests.data(),
4585 for (
unsigned int c = 0; c < static_cast<unsigned int>(
n_cells); ++c)
4590 MeshCellIteratorType mesh_it(
tria,
4595 const std_cxx17::optional<DataType> data =
pack(mesh_it);
4598 data_to_send.emplace_back(*data);
4599 cell_carries_data[c] =
true;
4607 sendbuffers[idx].resize(
sizeof(std::size_t));
4612 std::size_t size_of_send =
4616 std::memcpy(sendbuffers[idx].data(),
4618 sizeof(std::size_t));
4622 if (data_to_send.size() <
n_cells)
4628 ierr = MPI_Isend(sendbuffers[idx].data(),
4629 sendbuffers[idx].size(),
4634 &reply_requests[idx]);
4639 std::vector<char> receive;
4640 for (
unsigned int id = 0;
id < neighbor_cell_list.size(); ++id)
4643 int ierr = MPI_Probe(MPI_ANY_SOURCE,
4650 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4653 receive.resize(len);
4655 ierr = MPI_Recv(receive.data(),
4666 auto data_iterator = receive.begin();
4667 std::size_t size_of_received_data =
4668 Utilities::unpack<std::size_t>(data_iterator,
4669 data_iterator +
sizeof(std::size_t));
4670 data_iterator +=
sizeof(std::size_t);
4673 auto received_data = Utilities::unpack<std::vector<DataType>>(
4675 data_iterator + size_of_received_data,
4677 data_iterator += size_of_received_data;
4682 const std::vector<typename CellId::binary_type> &this_cell_list =
4683 neighbor_cell_list[status.MPI_SOURCE];
4685 std::vector<bool> cells_with_data;
4686 if (received_data.size() < this_cell_list.size())
4688 cells_with_data = Utilities::unpack<std::vector<bool>>(
4689 data_iterator, receive.end(),
false);
4695 auto received_data_iterator = received_data.begin();
4696 for (
unsigned int c = 0; c < this_cell_list.size(); ++c)
4697 if (cells_with_data.empty() || cells_with_data[c])
4702 MeshCellIteratorType cell(
tria,
4707 unpack(cell, *received_data_iterator);
4708 ++received_data_iterator;
4714 if (requests.size() > 0)
4717 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4720 if (reply_requests.size() > 0)
4722 const int ierr = MPI_Waitall(reply_requests.size(),
4723 reply_requests.data(),
4724 MPI_STATUSES_IGNORE);
4734 template <
typename DataType,
typename MeshType>
4737 const MeshType & mesh,
4738 const std::function<std_cxx17::optional<DataType>(
4739 const typename MeshType::active_cell_iterator &)> &
pack,
4740 const std::function<
void(
const typename MeshType::active_cell_iterator &,
4741 const DataType &)> &
unpack,
4742 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
4745 # ifndef DEAL_II_WITH_MPI
4752 internal::exchange_cell_data<DataType,
4754 typename MeshType::active_cell_iterator>(
4759 [&](
const auto &process) {
4760 for (
const auto &cell : mesh.active_cell_iterators())
4761 if (cell->is_ghost())
4762 process(cell, cell->subdomain_id());
4764 [](
const auto &
tria) {
return tria.ghost_owners(); });
4770 template <
typename DataType,
typename MeshType>
4773 const MeshType & mesh,
4774 const std::function<std_cxx17::optional<DataType>(
4775 const typename MeshType::level_cell_iterator &)> &
pack,
4776 const std::function<
void(
const typename MeshType::level_cell_iterator &,
4777 const DataType &)> &
unpack,
4778 const std::function<
bool(
const typename MeshType::level_cell_iterator &)>
4781 # ifndef DEAL_II_WITH_MPI
4788 internal::exchange_cell_data<DataType,
4790 typename MeshType::level_cell_iterator>(
4795 [&](
const auto &process) {
4796 for (
const auto &cell : mesh.cell_iterators())
4797 if (cell->is_ghost_on_level())
4798 process(cell, cell->level_subdomain_id());
4800 [](
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 DEAL_II_HOST SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
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_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
__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)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
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={})
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)
void random(DoFHandler< dim, spacedim > &dof_handler)
@ 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