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>
107 using type =
typename MeshType::active_cell_iterator;
114 template <
int dim,
int spacedim>
145 template <
int dim,
int spacedim>
172 template <
int dim,
int spacedim>
203 template <
int dim,
int spacedim>
218 template <
int dim,
int spacedim>
223 (ReferenceCells::get_hypercube<dim>()
225 .
template get_default_linear_mapping<dim, spacedim>()
227 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
241 template <
int dim,
int spacedim>
246 (ReferenceCells::get_hypercube<dim>()
248 .
template get_default_linear_mapping<dim, spacedim>()
250 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
317 template <
int dim,
int spacedim>
379 template <
int dim,
int spacedim>
400 template <
typename Iterator>
403 const Iterator &
object,
418 template <
int dim,
int spacedim>
420 tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>,
SubCellData>
445 template <
int dim,
int spacedim>
469 template <
int dim,
int spacedim>
474 std::vector<unsigned int> & considered_vertices,
475 const double tol = 1
e-12);
487 const double tol = 1
e-12);
507 template <
int dim,
int spacedim>
522 template <
int dim,
int spacedim>
610 template <
int dim,
typename Transformation,
int spacedim>
613 std::assignable_from<
625 template <
int dim,
int spacedim>
641 template <
int dim,
int spacedim>
681 const
unsigned int axis,
745 const
Function<dim,
double> *coefficient =
nullptr,
746 const
bool solve_for_absolute_positions = false);
753 template <
int dim,
int spacedim>
754 std::map<
unsigned int,
Point<spacedim>>
764 template <
int dim,
int spacedim>
766 scale(const
double scaling_factor,
783 template <
int dim,
int spacedim>
788 const
bool keep_boundary = true,
789 const
unsigned int seed =
boost::
random::mt19937::default_seed);
824 template <
int dim,
int spacedim>
827 const
bool isotropic = false,
828 const
unsigned int max_iterations = 100);
854 template <
int dim,
int spacedim>
857 const
double max_ratio = 1.6180339887,
858 const
unsigned int max_iterations = 5);
949 template <
int dim,
int spacedim>
952 const
double limit_angle_fraction = .75);
1014 template <
int dim,
int spacedim>
1017 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1018 std::vector<std::vector<Point<dim>>>,
1019 std::vector<std::vector<unsigned int>>>
1063 template <
int dim,
int spacedim>
1066 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1067 std::vector<std::vector<Point<dim>>>,
1068 std::vector<std::vector<unsigned int>>,
1069 std::vector<unsigned int>>
1149 template <
int dim,
int spacedim>
1152 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
1153 std::vector<std::vector<Point<dim>>>,
1154 std::vector<std::vector<unsigned int>>,
1155 std::vector<std::vector<Point<spacedim>>>,
1156 std::vector<std::vector<unsigned int>>>
1164 const double tolerance = 1
e-10);
1182 template <
int dim,
int spacedim>
1191 std::vector<std::tuple<std::pair<int, int>,
1221 std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
1246 template <
int dim,
int spacedim>
1252 const std::vector<bool> & marked_vertices,
1253 const double tolerance,
1254 const bool perform_handshake,
1255 const bool enforce_unique_mapping =
false);
1295 template <
int dim,
int spacedim>
1296 std::map<unsigned int, Point<spacedim>>
1300 (ReferenceCells::get_hypercube<dim>()
1302 .
template get_default_linear_mapping<dim, spacedim>()
1304 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
1317 template <
int spacedim>
1345 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1349 const MeshType<dim, spacedim> &mesh,
1350 const
Point<spacedim> & p,
1351 const std::vector<
bool> & marked_vertices = {});
1376 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1380 const
Mapping<dim, spacedim> & mapping,
1381 const MeshType<dim, spacedim> &mesh,
1382 const
Point<spacedim> & p,
1383 const std::vector<
bool> & marked_vertices = {});
1405 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1409 std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
1412 typename ::internal::
1413 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
1416 const unsigned int vertex_index);
1480 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1484 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1486 std::pair<typename ::internal::
1487 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1491 const MeshType<dim, spacedim> &mesh,
1493 const std::vector<bool> &marked_vertices = {},
1494 const double tolerance = 1.e-10);
1503 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1507 typename MeshType<dim, spacedim>::active_cell_iterator
1509 typename ::internal::
1510 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type
1514 const std::vector<bool> &marked_vertices = {},
1515 const double tolerance = 1.e-10);
1523 template <
int dim,
int spacedim>
1524 std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
1530 const double tolerance = 1.e-10);
1583 template <
int dim,
int spacedim>
1584 std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
1587 const Cache<dim, spacedim> &cache,
1591 const std::vector<bool> &marked_vertices = {},
1592 const double tolerance = 1.e-10);
1607 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1611 std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
Point<dim>>
1613 std::pair<typename ::internal::
1614 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1619 const MeshType<dim, spacedim> &mesh,
1622 std::set<
typename MeshType<dim, spacedim>::active_cell_iterator>>
1625 &vertex_to_cell_centers,
1626 const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint =
1627 typename MeshType<dim, spacedim>::active_cell_iterator(),
1628 const std::vector<bool> &marked_vertices = {},
1631 const double tolerance = 1.e-10,
1633 std::pair<BoundingBox<spacedim>,
1635 *relevant_cell_bounding_boxes_rtree =
nullptr);
1658 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1662 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1665 std::vector<std::pair<
1666 typename ::internal::
1667 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1672 const MeshType<dim, spacedim> &mesh,
1674 const double tolerance,
1675 const std::pair<
typename MeshType<dim, spacedim>::active_cell_iterator,
1684 template <
int dim,
template <
int,
int>
class MeshType,
int spacedim>
1688 std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
1691 std::vector<std::pair<
1692 typename ::internal::
1693 ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
1698 const MeshType<dim, spacedim> &mesh,
1700 const double tolerance = 1
e-10,
1701 const std::vector<bool> & marked_vertices = {});
1724 template <
class MeshType>
1727 const typename MeshType::cell_iterator &cell);
1753 template <class MeshType>
1756 const typename MeshType::active_cell_iterator & cell,
1757 std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
1808 template <class MeshType>
1812 const MeshType &mesh,
1813 const std::function<
bool(const typename MeshType::active_cell_iterator &)>
1824 template <class MeshType>
1828 const MeshType &mesh,
1829 const std::function<
bool(const typename MeshType::cell_iterator &)>
1831 const
unsigned int level);
1846 template <class MeshType>
1900 template <class MeshType>
1904 const MeshType &mesh,
1905 const std::function<
bool(const typename MeshType::active_cell_iterator &)>
1907 const
double layer_thickness);
1931 template <class MeshType>
1935 const MeshType &mesh,
1936 const
double layer_thickness);
1953 template <class MeshType>
1956 Point<MeshType::space_dimension>,
1959 const std::function<
bool(
1960 const typename MeshType::
1961 active_cell_iterator &)>
2034 template <class MeshType>
2038 const MeshType &mesh,
2039 const std::function<
bool(const typename MeshType::active_cell_iterator &)>
2041 const
unsigned int refinement_level = 0,
2042 const
bool allow_merge = false,
2072 template <
int spacedim>
2074 std::tuple<std::vector<std::vector<unsigned int>>,
2075 std::map<unsigned int, unsigned int>,
2076 std::map<unsigned int, std::vector<unsigned int>>>
2119 template <
int spacedim>
2121 std::tuple<std::map<unsigned int, std::vector<unsigned int>>,
2122 std::map<unsigned int, unsigned int>,
2123 std::map<unsigned int, std::vector<unsigned int>>>
2140 template <
int dim,
int spacedim>
2142 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
2157 template <
int dim,
int spacedim>
2158 std::vector<std::vector<Tensor<1, spacedim>>>
2173 template <
int dim,
int spacedim>
2179 (ReferenceCells::get_hypercube<dim>()
2181 .
template get_default_linear_mapping<dim, spacedim>()
2183 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
2198 template <
int dim,
int spacedim>
2199 std::map<unsigned int, types::global_vertex_index>
2214 template <
int dim,
int spacedim>
2215 std::pair<unsigned int, double>
2233 template <
int dim,
int spacedim>
2247 template <
int dim,
int spacedim>
2261 template <
int dim,
int spacedim>
2265 const unsigned int level,
2288 template <
int dim,
int spacedim>
2305 template <
int dim,
int spacedim>
2308 const std::vector<unsigned int> &cell_weights,
2358 template <
int dim,
int spacedim>
2376 template <
int dim,
int spacedim>
2379 const std::vector<unsigned int> &cell_weights,
2399 template <
int dim,
int spacedim>
2403 const bool group_siblings =
true);
2416 template <
int dim,
int spacedim>
2427 template <
int dim,
int spacedim>
2428 std::vector<types::subdomain_id>
2430 const std::vector<CellId> & cell_ids);
2442 template <
int dim,
int spacedim>
2445 std::vector<types::subdomain_id> & subdomain);
2461 template <
int dim,
int spacedim>
2496 template <
int dim,
int spacedim>
2539 template <
typename MeshType>
2541 std::list<std::pair<
2542 typename MeshType::cell_iterator,
2557 template <
int dim,
int spacedim>
2571 template <
typename MeshType>
2596 template <
int dim,
int spacedim>
2652 template <
class MeshType>
2655 const typename MeshType::active_cell_iterator &cell);
2679 template <
class Container>
2680 std::vector<typename Container::cell_iterator>
2682 const std::vector<typename Container::active_cell_iterator> &patch_cells);
2750 template <
class Container>
2753 const std::vector<typename Container::active_cell_iterator> &patch,
2755 &local_triangulation,
2758 Container::space_dimension>::active_cell_iterator,
2759 typename Container::active_cell_iterator> &patch_to_global_tria_map);
2792 template <
int dim,
int spacedim>
2795 std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
2811 template <
typename CellIterator>
2917 template <
typename FaceIterator>
2920 std::bitset<3> & orientation,
2921 const FaceIterator & face1,
2922 const FaceIterator & face2,
2923 const unsigned int direction,
2932 template <
typename FaceIterator>
2935 const FaceIterator & face1,
2936 const FaceIterator & face2,
2937 const unsigned int direction,
2999 template <
typename MeshType>
3002 const MeshType & mesh,
3005 const unsigned int direction,
3035 template <
typename MeshType>
3038 const MeshType & mesh,
3040 const
unsigned int direction,
3043 const ::
Tensor<1, MeshType::space_dimension> &offset =
3044 ::
Tensor<1, MeshType::space_dimension>(),
3073 template <
int dim,
int spacedim>
3076 const
bool reset_boundary_ids = false);
3099 template <
int dim,
int spacedim>
3136 template <
int dim,
int spacedim>
3139 const bool compute_face_ids =
false);
3165 template <
int dim,
int spacedim>
3170 const std::set<types::manifold_id> &)> &disambiguation_function =
3171 [](
const std::set<types::manifold_id> &manifold_ids) {
3172 if (manifold_ids.size() == 1)
3173 return *manifold_ids.
begin();
3177 bool overwrite_only_flat_manifold_ids =
true);
3264 template <
typename DataType,
typename MeshType>
3267 const MeshType & mesh,
3268 const std::function<
std_cxx17::optional<DataType>(
3269 const typename MeshType::active_cell_iterator &)> &
pack,
3270 const std::function<
void(const typename MeshType::active_cell_iterator &,
3271 const DataType &)> &
unpack,
3272 const std::function<
bool(const typename MeshType::active_cell_iterator &)>
3274 always_return<typename MeshType::active_cell_iterator,
bool>{
true});
3286 template <
typename DataType,
typename MeshType>
3289 const MeshType & mesh,
3290 const std::function<
std_cxx17::optional<DataType>(
3291 const typename MeshType::level_cell_iterator &)> &
pack,
3292 const std::function<
void(const typename MeshType::level_cell_iterator &,
3293 const DataType &)> &
unpack,
3294 const std::function<
bool(const typename MeshType::level_cell_iterator &)> &
3295 cell_filter =
always_return<typename MeshType::level_cell_iterator,
bool>{
3309 template <
int spacedim>
3310 std::vector<std::vector<BoundingBox<spacedim>>>
3313 const MPI_Comm & mpi_communicator);
3347 template <
int spacedim>
3351 const MPI_Comm & mpi_communicator);
3370 template <
int dim,
int spacedim>
3396 template <
int dim,
int spacedim>
3397 std::map<unsigned int, std::set<::types::subdomain_id>>
3416 template <
int dim,
typename T>
3437 template <
class Archive>
3439 save(Archive &ar,
const unsigned int)
const
3443 Assert(cell_ids.size() == data.size(),
3446 const std::size_t
n_cells = cell_ids.size();
3448 for (
const auto &
id : cell_ids)
3451 ar & binary_cell_id;
3463 template <
class Archive>
3465 load(Archive &ar,
const unsigned int)
3473 for (
unsigned int c = 0; c <
n_cells; ++c)
3477 cell_ids.emplace_back(value);
3489 template <
class Archive>
3495 BOOST_SERIALIZATION_SPLIT_MEMBER()
3521 template <
int dim,
typename VectorType>
3550 const VectorType & ls_vector,
3551 const double iso_level,
3561 const VectorType & ls_vector,
3562 const double iso_level,
3576 const VectorType & ls_vector,
3577 const double iso_level,
3587 const VectorType & ls_vector,
3588 const double iso_level,
3605 const double iso_level,
3608 const bool write_back_cell_data =
true)
const;
3616 const std::vector<unsigned int> &,
3633 const std::vector<
Point<2>> & points,
3634 const std::vector<unsigned int> &mask,
3635 const double iso_level,
3638 const bool write_back_cell_data)
const;
3645 const std::vector<
Point<3>> & points,
3646 const std::vector<unsigned int> &mask,
3647 const double iso_level,
3650 const bool write_back_cell_data)
const;
3683 <<
"The number of partitions you gave is " << arg1
3684 <<
", but must be greater than zero.");
3690 <<
"The subdomain id " << arg1
3691 <<
" has no cells associated with it.");
3702 <<
"The scaling factor must be positive, but it is " << arg1
3710 <<
"The given vertex with index " << arg1
3711 <<
" is not used in the given triangulation.");
3728 "The edges of the mesh are not consistently orientable.");
3765 template <
int dim,
typename Transformation,
int spacedim>
3768 std::assignable_from<
3771 void
transform(const Transformation & transformation,
3774 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
3785 for (; cell != endc; ++cell)
3786 for (
const unsigned int v : cell->vertex_indices())
3787 if (treated_vertices[cell->vertex_index(v)] ==
false)
3790 cell->vertex(v) = transformation(cell->vertex(v));
3792 treated_vertices[cell->vertex_index(v)] =
true;
3802 for (; cell != endc; ++cell)
3803 for (
const unsigned int face : cell->face_indices())
3804 if (cell->face(face)->has_children() &&
3805 !cell->face(face)->at_boundary())
3807 Assert(cell->reference_cell() ==
3808 ReferenceCells::get_hypercube<dim>(),
3812 cell->face(face)->child(0)->vertex(1) =
3813 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3822 for (; cell != endc; ++cell)
3823 for (
const unsigned int face : cell->face_indices())
3824 if (cell->face(face)->has_children() &&
3825 !cell->face(face)->at_boundary())
3827 Assert(cell->reference_cell() ==
3828 ReferenceCells::get_hypercube<dim>(),
3832 cell->face(face)->child(0)->vertex(1) =
3833 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
3835 cell->face(face)->child(0)->vertex(2) =
3836 (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) /
3838 cell->face(face)->child(1)->vertex(3) =
3839 (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) /
3841 cell->face(face)->child(2)->vertex(3) =
3842 (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3846 cell->face(face)->child(0)->vertex(3) =
3847 (cell->face(face)->vertex(0) + cell->face(face)->vertex(1) +
3848 cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) /
3859 template <
class MeshType>
3862 const typename MeshType::cell_iterator &cell)
3864 std::vector<typename MeshType::active_cell_iterator> child_cells;
3866 if (cell->has_children())
3868 for (
unsigned int child = 0; child < cell->n_children(); ++child)
3869 if (cell->child(child)->has_children())
3871 const std::vector<typename MeshType::active_cell_iterator>
3872 children = get_active_child_cells<MeshType>(cell->child(child));
3873 child_cells.insert(child_cells.end(),
3878 child_cells.push_back(cell->child(child));
3886 template <
class MeshType>
3889 const typename MeshType::active_cell_iterator & cell,
3890 std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
3892 active_neighbors.clear();
3893 for (
const unsigned int n : cell->face_indices())
3894 if (!cell->at_boundary(n))
3896 if (MeshType::dimension == 1)
3903 typename MeshType::cell_iterator neighbor_child =
3905 if (!neighbor_child->is_active())
3907 while (neighbor_child->has_children())
3908 neighbor_child = neighbor_child->child(n == 0 ? 1 : 0);
3910 Assert(neighbor_child->neighbor(n == 0 ? 1 : 0) == cell,
3913 active_neighbors.push_back(neighbor_child);
3917 if (cell->face(n)->has_children())
3921 for (
unsigned int c = 0;
3922 c < cell->face(n)->n_active_descendants();
3924 active_neighbors.push_back(
3925 cell->neighbor_child_on_subface(n, c));
3931 active_neighbors.push_back(cell->neighbor(n));
3939 template <
typename CellIterator>
3943 return sizeof(*this) +
matrix.memory_consumption();
3950 namespace ProjectToObject
3964 struct CrossDerivative
3966 const unsigned int direction_0;
3967 const unsigned int direction_1;
3969 CrossDerivative(
const unsigned int d0,
const unsigned int d1);
3972 inline CrossDerivative::CrossDerivative(
const unsigned int d0,
3973 const unsigned int d1)
3984 template <
typename F>
3986 centered_first_difference(
const double center,
3990 return (f(
center + step) - f(
center - step)) / (2.0 * step);
3999 template <
typename F>
4001 centered_second_difference(
const double center,
4020 template <
int structdim,
typename F>
4023 const CrossDerivative cross_derivative,
4029 simplex_vector[cross_derivative.direction_0] = 0.5 * step;
4030 simplex_vector[cross_derivative.direction_1] = -0.5 * step;
4031 return (-4.0 * f(
center) - 1.0 * f(
center + simplex_vector) -
4032 1.0 / 3.0 * f(
center - simplex_vector) +
4033 16.0 / 3.0 * f(
center + 0.5 * simplex_vector)) /
4045 template <
int spacedim,
int structdim,
typename F>
4048 const unsigned int row_n,
4049 const unsigned int dependent_direction,
4056 dependent_direction <
4058 ExcMessage(
"This function assumes that the last weight is a "
4059 "dependent variable (and hence we cannot take its "
4060 "derivative directly)."));
4061 Assert(row_n != dependent_direction,
4063 "We cannot differentiate with respect to the variable "
4064 "that is assumed to be dependent."));
4068 {row_n, dependent_direction},
center, step, f);
4070 for (
unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
4072 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
4081 template <
typename Iterator,
int spacedim,
int structdim>
4083 project_to_d_linear_object(
const Iterator &
object,
4118 for (
unsigned int d = 0;
d < structdim; ++
d)
4123 x_k +=
object->vertex(i) *
4129 for (
const unsigned int i :
4132 (x_k - trial_point) * object->vertex(i) *
4137 for (
const unsigned int i :
4139 for (
const unsigned int j :
4147 H_k += (
object->vertex(i) *
object->vertex(j)) * tmp;
4154 for (
const unsigned int i :
4156 x_k +=
object->vertex(i) *
4159 if (delta_xi.
norm() < 1
e-7)
4175 template <
int structdim>
4182 static const std::size_t n_vertices_per_cell =
4184 n_independent_components;
4185 std::array<double, n_vertices_per_cell> copied_weights;
4186 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
4188 copied_weights[i] = v[i];
4189 if (v[i] < 0.0 || v[i] > 1.0)
4194 std::sort(copied_weights.begin(), copied_weights.end());
4196 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
4201 template <
typename Iterator>
4204 const Iterator &
object,
4207 const int spacedim = Iterator::AccessorType::space_dimension;
4208 const int structdim = Iterator::AccessorType::structure_dimension;
4212 if (structdim >= spacedim)
4213 return projected_point;
4214 else if (structdim == 1 || structdim == 2)
4216 using namespace internal::ProjectToObject;
4221 const int dim = Iterator::AccessorType::dimension;
4224 &manifold) !=
nullptr)
4227 project_to_d_linear_object<Iterator, spacedim, structdim>(
4228 object, trial_point);
4273 const double step_size =
object->diameter() / 64.0;
4275 constexpr
unsigned int n_vertices_per_cell =
4278 std::array<Point<spacedim>, n_vertices_per_cell>
vertices;
4279 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4281 vertices[vertex_n] = object->vertex(vertex_n);
4283 auto get_point_from_weights =
4286 return object->get_manifold().get_new_point(
4294 double guess_weights_sum = 0.0;
4295 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
4298 const double distance =
4300 if (distance == 0.0)
4302 guess_weights = 0.0;
4303 guess_weights[vertex_n] = 1.0;
4304 guess_weights_sum = 1.0;
4309 guess_weights[vertex_n] = 1.0 / distance;
4310 guess_weights_sum += guess_weights[vertex_n];
4313 guess_weights /= guess_weights_sum;
4314 Assert(internal::weights_are_ok<structdim>(guess_weights),
4323 for (
unsigned int outer_n = 0; outer_n < 40; ++outer_n)
4325 const unsigned int dependent_direction =
4326 n_vertices_per_cell - 1;
4328 for (
unsigned int row_n = 0; row_n < n_vertices_per_cell;
4331 if (row_n != dependent_direction)
4333 current_gradient[row_n] =
4334 gradient_entry<spacedim, structdim>(
4336 dependent_direction,
4340 get_point_from_weights);
4342 current_gradient[dependent_direction] -=
4343 current_gradient[row_n];
4361 double gradient_weight = -0.5;
4362 auto gradient_weight_objective_function =
4363 [&](
const double gradient_weight_guess) ->
double {
4364 return (trial_point -
4365 get_point_from_weights(guess_weights +
4366 gradient_weight_guess *
4371 for (
unsigned int inner_n = 0; inner_n < 10; ++inner_n)
4373 const double update_numerator = centered_first_difference(
4376 gradient_weight_objective_function);
4377 const double update_denominator =
4378 centered_second_difference(
4381 gradient_weight_objective_function);
4385 if (std::abs(update_denominator) == 0.0)
4388 gradient_weight - update_numerator / update_denominator;
4393 if (std::abs(gradient_weight) > 10)
4395 gradient_weight = -10.0;
4405 guess_weights + gradient_weight * current_gradient;
4407 double new_gradient_weight = gradient_weight;
4408 for (
unsigned int iteration_count = 0; iteration_count < 40;
4411 if (internal::weights_are_ok<structdim>(tentative_weights))
4414 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
4416 if (tentative_weights[i] < 0.0)
4418 tentative_weights -=
4419 (tentative_weights[i] / current_gradient[i]) *
4422 if (tentative_weights[i] < 0.0 ||
4423 1.0 < tentative_weights[i])
4425 new_gradient_weight /= 2.0;
4428 new_gradient_weight * current_gradient;
4435 if (!internal::weights_are_ok<structdim>(tentative_weights))
4440 if (get_point_from_weights(tentative_weights)
4441 .distance(trial_point) <
4442 get_point_from_weights(guess_weights).distance(trial_point))
4443 guess_weights = tentative_weights;
4446 Assert(internal::weights_are_ok<structdim>(guess_weights),
4449 Assert(internal::weights_are_ok<structdim>(guess_weights),
4451 projected_point = get_point_from_weights(guess_weights);
4461 for (
unsigned int line_n = 0;
4462 line_n < GeometryInfo<structdim>::lines_per_cell;
4465 line_projections[line_n] =
4468 std::sort(line_projections.begin(),
4469 line_projections.end(),
4471 return a.distance(trial_point) <
4472 b.distance(trial_point);
4474 if (line_projections[0].distance(trial_point) <
4475 projected_point.
distance(trial_point))
4476 projected_point = line_projections[0];
4482 return projected_point;
4485 return projected_point;
4492 template <
typename DataType,
4494 typename MeshCellIteratorType>
4496 inline void exchange_cell_data(
4497 const MeshType &mesh,
4498 const std::function<
4499 std_cxx17::optional<DataType>(
const MeshCellIteratorType &)> &
pack,
4500 const std::function<
void(
const MeshCellIteratorType &,
const DataType &)>
4502 const std::function<
bool(
const MeshCellIteratorType &)> & cell_filter,
4503 const std::function<
void(
4504 const std::function<
void(
const MeshCellIteratorType &,
4506 const std::function<std::set<types::subdomain_id>(
4508 MeshType::space_dimension> &)>
4509 &compute_ghost_owners)
4511 # ifndef DEAL_II_WITH_MPI
4516 (void)process_cells;
4517 (void)compute_ghost_owners;
4520 constexpr
int dim = MeshType::dimension;
4521 constexpr
int spacedim = MeshType::space_dimension;
4524 &mesh.get_triangulation());
4528 "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
4530 if (
const auto tria =
dynamic_cast<
4532 &mesh.get_triangulation()))
4535 tria->with_artificial_cells(),
4537 "The functions GridTools::exchange_cell_data_to_ghosts() and "
4538 "GridTools::exchange_cell_data_to_level_ghosts() can only "
4539 "operate on a single layer of ghost cells. However, you have "
4540 "given a Triangulation object of type "
4541 "parallel::shared::Triangulation without artificial cells "
4542 "resulting in arbitrary numbers of ghost layers."));
4546 std::set<::types::subdomain_id> ghost_owners =
4547 compute_ghost_owners(*
tria);
4549 std::vector<typename CellId::binary_type>>
4552 for (
const auto ghost_owner : ghost_owners)
4553 neighbor_cell_list[ghost_owner] = {};
4555 process_cells([&](
const auto &cell,
const auto key) ->
void {
4556 if (cell_filter(cell))
4558 constexpr
int spacedim = MeshType::space_dimension;
4559 neighbor_cell_list[key].emplace_back(
4560 cell->id().template to_binary<spacedim>());
4564 Assert(ghost_owners.size() == neighbor_cell_list.size(),
4576 const int mpi_tag_reply =
4580 std::vector<MPI_Request> requests(ghost_owners.size());
4582 unsigned int idx = 0;
4583 for (
const auto &it : neighbor_cell_list)
4586 const int ierr = MPI_Isend(it.second.data(),
4587 it.second.size() *
sizeof(it.second[0]),
4599 std::vector<MPI_Request> reply_requests(ghost_owners.size());
4600 std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
4602 for (
unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
4605 int ierr = MPI_Probe(MPI_ANY_SOURCE,
4612 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4619 std::vector<typename CellId::binary_type> cells_with_requests(
4621 std::vector<DataType> data_to_send;
4622 data_to_send.reserve(
n_cells);
4623 std::vector<bool> cell_carries_data(
n_cells,
false);
4625 ierr = MPI_Recv(cells_with_requests.data(),
4635 for (
unsigned int c = 0; c < static_cast<unsigned int>(
n_cells); ++c)
4640 MeshCellIteratorType mesh_it(
tria,
4645 const std_cxx17::optional<DataType> data =
pack(mesh_it);
4648 data_to_send.emplace_back(*data);
4649 cell_carries_data[c] =
true;
4657 sendbuffers[idx].resize(
sizeof(std::size_t));
4662 std::size_t size_of_send =
4666 std::memcpy(sendbuffers[idx].data(),
4668 sizeof(std::size_t));
4672 if (data_to_send.size() <
n_cells)
4678 ierr = MPI_Isend(sendbuffers[idx].data(),
4679 sendbuffers[idx].size(),
4684 &reply_requests[idx]);
4689 std::vector<char> receive;
4690 for (
unsigned int id = 0;
id < neighbor_cell_list.size(); ++id)
4693 int ierr = MPI_Probe(MPI_ANY_SOURCE,
4700 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
4703 receive.resize(len);
4705 ierr = MPI_Recv(receive.data(),
4716 auto data_iterator = receive.begin();
4717 std::size_t size_of_received_data =
4718 Utilities::unpack<std::size_t>(data_iterator,
4719 data_iterator +
sizeof(std::size_t));
4720 data_iterator +=
sizeof(std::size_t);
4723 auto received_data = Utilities::unpack<std::vector<DataType>>(
4725 data_iterator + size_of_received_data,
4727 data_iterator += size_of_received_data;
4732 const std::vector<typename CellId::binary_type> &this_cell_list =
4733 neighbor_cell_list[status.MPI_SOURCE];
4735 std::vector<bool> cells_with_data;
4736 if (received_data.size() < this_cell_list.size())
4738 cells_with_data = Utilities::unpack<std::vector<bool>>(
4739 data_iterator, receive.end(),
false);
4745 auto received_data_iterator = received_data.begin();
4746 for (
unsigned int c = 0; c < this_cell_list.size(); ++c)
4747 if (cells_with_data.empty() || cells_with_data[c])
4752 MeshCellIteratorType cell(
tria,
4757 unpack(cell, *received_data_iterator);
4758 ++received_data_iterator;
4764 if (requests.size() > 0)
4767 MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
4770 if (reply_requests.size() > 0)
4772 const int ierr = MPI_Waitall(reply_requests.size(),
4773 reply_requests.data(),
4774 MPI_STATUSES_IGNORE);
4784 template <
typename DataType,
typename MeshType>
4787 const MeshType & mesh,
4788 const std::function<std_cxx17::optional<DataType>(
4789 const typename MeshType::active_cell_iterator &)> &
pack,
4790 const std::function<
void(
const typename MeshType::active_cell_iterator &,
4791 const DataType &)> &
unpack,
4792 const std::function<
bool(
const typename MeshType::active_cell_iterator &)>
4795 # ifndef DEAL_II_WITH_MPI
4802 internal::exchange_cell_data<DataType,
4804 typename MeshType::active_cell_iterator>(
4809 [&](
const auto &process) {
4810 for (
const auto &cell : mesh.active_cell_iterators())
4811 if (cell->is_ghost())
4812 process(cell, cell->subdomain_id());
4814 [](
const auto &
tria) {
return tria.ghost_owners(); });
4820 template <
typename DataType,
typename MeshType>
4823 const MeshType & mesh,
4824 const std::function<std_cxx17::optional<DataType>(
4825 const typename MeshType::level_cell_iterator &)> &
pack,
4826 const std::function<
void(
const typename MeshType::level_cell_iterator &,
4827 const DataType &)> &
unpack,
4828 const std::function<
bool(
const typename MeshType::level_cell_iterator &)>
4831 # ifndef DEAL_II_WITH_MPI
4838 internal::exchange_cell_data<DataType,
4840 typename MeshType::level_cell_iterator>(
4845 [&](
const auto &process) {
4846 for (
const auto &cell : mesh.cell_iterators())
4847 if (cell->is_ghost_on_level())
4848 process(cell, cell->level_subdomain_id());
4850 [](
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)
constexpr bool is_triangulation_or_dof_handler
concept is_triangulation_or_dof_handler
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