15#ifndef dealii_grid_tools_geometry_h
16#define dealii_grid_tools_geometry_h
52 template <
int dim,
int spacedim>
79 template <
int dim,
int spacedim>
110 template <
int dim,
int spacedim>
125 template <
int dim,
int spacedim>
130 (ReferenceCells::get_hypercube<dim>()
132 .
template get_default_linear_mapping<dim, spacedim>()
134 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
148 template <
int dim,
int spacedim>
153 (ReferenceCells::get_hypercube<dim>()
155 .
template get_default_linear_mapping<dim, spacedim>()
157 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
196 template <
int dim,
int spacedim>
197 std::pair<unsigned int, double>
223 template <
int dim,
int spacedim>
288 template <
int dim,
int spacedim>
309 template <
typename MeshType>
315 const std::function<
bool(
316 const typename MeshType::
317 active_cell_iterator &)>
337 template <
typename Iterator>
340 const Iterator &
object,
350 namespace ProjectToObject
364 struct CrossDerivative
366 const unsigned int direction_0;
367 const unsigned int direction_1;
369 CrossDerivative(
const unsigned int d0,
const unsigned int d1);
372 inline CrossDerivative::CrossDerivative(
const unsigned int d0,
373 const unsigned int d1)
384 template <
typename F>
386 centered_first_difference(
const double center,
388 const F &f) ->
decltype(f(center) - f(center))
390 return (f(center + step) - f(center - step)) / (2.0 * step);
399 template <
typename F>
401 centered_second_difference(
const double center,
403 const F &f) ->
decltype(f(center) - f(center))
405 return (f(center + step) - 2.0 * f(center) + f(center - step)) /
420 template <
int structdim,
typename F>
423 const CrossDerivative cross_derivative,
426 const F &f) ->
decltype(f(center) - f(center))
429 simplex_vector[cross_derivative.direction_0] = 0.5 * step;
430 simplex_vector[cross_derivative.direction_1] = -0.5 * step;
431 return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
432 1.0 / 3.0 * f(center - simplex_vector) +
433 16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
445 template <
int spacedim,
int structdim,
typename F>
448 const unsigned int row_n,
449 const unsigned int dependent_direction,
456 dependent_direction <
458 ExcMessage(
"This function assumes that the last weight is a "
459 "dependent variable (and hence we cannot take its "
460 "derivative directly)."));
461 Assert(row_n != dependent_direction,
463 "We cannot differentiate with respect to the variable "
464 "that is assumed to be dependent."));
468 {row_n, dependent_direction}, center, step, f);
470 for (
unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
472 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
481 template <
typename Iterator,
int spacedim,
int structdim>
483 project_to_d_linear_object(
const Iterator &
object,
519 for (
unsigned int d = 0;
d < structdim; ++
d)
524 x_k += object->vertex(i) *
525 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
530 for (
const unsigned int i :
533 (x_k - trial_point) * object->vertex(i) *
534 GeometryInfo<structdim>::d_linear_shape_function_gradient(xi,
538 for (
const unsigned int i :
540 for (const unsigned
int j :
548 H_k += (
object->vertex(i) *
object->vertex(j)) * tmp;
555 for (
const unsigned int i :
557 x_k += object->vertex(i) *
558 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
560 if (delta_xi.
norm() < 1e-7)
571 template <
int structdim>
578 static const std::size_t n_vertices_per_cell =
580 n_independent_components;
581 std::array<double, n_vertices_per_cell> copied_weights;
582 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
584 copied_weights[i] = v[i];
585 if (v[i] < 0.0 || v[i] > 1.0)
590 std::sort(copied_weights.begin(), copied_weights.end());
592 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
599 template <
typename Iterator>
602 const Iterator &
object,
605 const int spacedim = Iterator::AccessorType::space_dimension;
606 const int structdim = Iterator::AccessorType::structure_dimension;
610 if (structdim >= spacedim)
611 return projected_point;
612 else if (structdim == 1 || structdim == 2)
614 using namespace internal::ProjectToObject;
619 const int dim = Iterator::AccessorType::dimension;
622 &manifold) !=
nullptr)
625 project_to_d_linear_object<Iterator, spacedim, structdim>(
626 object, trial_point);
671 const double step_size =
object->diameter() / 64.0;
673 constexpr unsigned int n_vertices_per_cell =
676 std::array<Point<spacedim>, n_vertices_per_cell> vertices;
677 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
679 vertices[vertex_n] = object->vertex(vertex_n);
681 auto get_point_from_weights =
684 return object->get_manifold().get_new_point(
692 double guess_weights_sum = 0.0;
693 for (
unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
696 const double distance =
697 vertices[vertex_n].distance(trial_point);
701 guess_weights[vertex_n] = 1.0;
702 guess_weights_sum = 1.0;
707 guess_weights[vertex_n] = 1.0 / distance;
708 guess_weights_sum += guess_weights[vertex_n];
711 guess_weights /= guess_weights_sum;
712 Assert(internal::weights_are_ok<structdim>(guess_weights),
721 for (
unsigned int outer_n = 0; outer_n < 40; ++outer_n)
723 const unsigned int dependent_direction =
724 n_vertices_per_cell - 1;
726 for (
unsigned int row_n = 0; row_n < n_vertices_per_cell;
729 if (row_n != dependent_direction)
731 current_gradient[row_n] =
732 gradient_entry<spacedim, structdim>(
738 get_point_from_weights);
740 current_gradient[dependent_direction] -=
741 current_gradient[row_n];
759 double gradient_weight = -0.5;
760 auto gradient_weight_objective_function =
761 [&](
const double gradient_weight_guess) ->
double {
762 return (trial_point -
763 get_point_from_weights(guess_weights +
764 gradient_weight_guess *
769 for (
unsigned int inner_n = 0; inner_n < 10; ++inner_n)
771 const double update_numerator = centered_first_difference(
774 gradient_weight_objective_function);
775 const double update_denominator =
776 centered_second_difference(
779 gradient_weight_objective_function);
783 if (
std::abs(update_denominator) == 0.0)
786 gradient_weight - update_numerator / update_denominator;
793 gradient_weight = -10.0;
803 guess_weights + gradient_weight * current_gradient;
805 double new_gradient_weight = gradient_weight;
806 for (
unsigned int iteration_count = 0; iteration_count < 40;
809 if (internal::weights_are_ok<structdim>(tentative_weights))
812 for (
unsigned int i = 0; i < n_vertices_per_cell; ++i)
814 if (tentative_weights[i] < 0.0)
817 (tentative_weights[i] / current_gradient[i]) *
820 if (tentative_weights[i] < 0.0 ||
821 1.0 < tentative_weights[i])
823 new_gradient_weight /= 2.0;
826 new_gradient_weight * current_gradient;
833 if (!internal::weights_are_ok<structdim>(tentative_weights))
838 if (get_point_from_weights(tentative_weights)
839 .distance(trial_point) <
840 get_point_from_weights(guess_weights).distance(trial_point))
841 guess_weights = tentative_weights;
844 Assert(internal::weights_are_ok<structdim>(guess_weights),
847 Assert(internal::weights_are_ok<structdim>(guess_weights),
849 projected_point = get_point_from_weights(guess_weights);
859 for (
unsigned int line_n = 0;
860 line_n < GeometryInfo<structdim>::lines_per_cell;
863 line_projections[line_n] =
866 std::sort(line_projections.begin(),
867 line_projections.end(),
869 return a.distance(trial_point) <
870 b.distance(trial_point);
872 if (line_projections[0].distance(trial_point) <
873 projected_point.
distance(trial_point))
874 projected_point = line_projections[0];
880 return projected_point;
883 return projected_point;
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Abstract base class for mapping classes.
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)