deal.II version GIT relicensing-2248-g6a8f632045 2024-12-12 07:30:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
grid_tools_geometry.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_grid_tools_geometry_h
16#define dealii_grid_tools_geometry_h
17
18#include <deal.II/base/config.h>
19
21#include <deal.II/base/point.h>
23#include <deal.II/base/tensor.h>
24
25#include <deal.II/fe/mapping.h>
26
28#include <deal.II/grid/tria.h>
29
30#include <algorithm>
31#include <array>
32#include <cmath>
33#include <utility>
34#include <vector>
35
37
38namespace GridTools
39{
51 template <int dim, int spacedim>
52 double
54
78 template <int dim, int spacedim>
79 double
81
109 template <int dim, int spacedim>
110 double
112 const Mapping<dim, spacedim> &mapping);
113
124 template <int dim, int spacedim>
125 double
128 const Mapping<dim, spacedim> &mapping =
129 (ReferenceCells::get_hypercube<dim>()
130#ifndef _MSC_VER
131 .template get_default_linear_mapping<dim, spacedim>()
132#else
133 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
134#endif
135 ));
136
147 template <int dim, int spacedim>
148 double
151 const Mapping<dim, spacedim> &mapping =
152 (ReferenceCells::get_hypercube<dim>()
153#ifndef _MSC_VER
154 .template get_default_linear_mapping<dim, spacedim>()
155#else
156 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
157#endif
158 ));
159
179 template <int dim>
180 double
181 cell_measure(const std::vector<Point<dim>> &all_vertices,
183
195 template <int dim, int spacedim>
196 std::pair<unsigned int, double>
199
222 template <int dim, int spacedim>
223 std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
225
255 template <int dim>
259 const Quadrature<dim> &quadrature);
260
268 template <int dim>
269 double
272 const Quadrature<dim> &quadrature);
273
287 template <int dim, int spacedim>
290
308 template <typename MeshType>
310 std::pair<
312 Point<MeshType::
313 space_dimension>> compute_bounding_box(const MeshType &mesh,
314 const std::function<bool(
315 const typename MeshType::
316 active_cell_iterator &)>
317 &predicate);
318
336 template <typename Iterator>
339 const Iterator &object,
342} // namespace GridTools
343
344#ifndef DOXYGEN
345namespace GridTools
346{
347 namespace internal
348 {
349 namespace ProjectToObject
350 {
363 struct CrossDerivative
364 {
365 const unsigned int direction_0;
366 const unsigned int direction_1;
367
368 CrossDerivative(const unsigned int d0, const unsigned int d1);
369 };
370
371 inline CrossDerivative::CrossDerivative(const unsigned int d0,
372 const unsigned int d1)
373 : direction_0(d0)
374 , direction_1(d1)
375 {}
376
377
378
383 template <typename F>
384 inline auto
385 centered_first_difference(const double center,
386 const double step,
387 const F &f) -> decltype(f(center) - f(center))
388 {
389 return (f(center + step) - f(center - step)) / (2.0 * step);
390 }
391
392
393
398 template <typename F>
399 inline auto
400 centered_second_difference(const double center,
401 const double step,
402 const F &f) -> decltype(f(center) - f(center))
403 {
404 return (f(center + step) - 2.0 * f(center) + f(center - step)) /
405 (step * step);
406 }
407
408
409
419 template <int structdim, typename F>
420 inline auto
421 cross_stencil(
422 const CrossDerivative cross_derivative,
424 const double step,
425 const F &f) -> decltype(f(center) - f(center))
426 {
428 simplex_vector[cross_derivative.direction_0] = 0.5 * step;
429 simplex_vector[cross_derivative.direction_1] = -0.5 * step;
430 return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
431 1.0 / 3.0 * f(center - simplex_vector) +
432 16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
433 step;
434 }
435
436
437
444 template <int spacedim, int structdim, typename F>
445 inline double
446 gradient_entry(
447 const unsigned int row_n,
448 const unsigned int dependent_direction,
449 const Point<spacedim> &p0,
451 const double step,
452 const F &f)
453 {
455 dependent_direction <
457 ExcMessage("This function assumes that the last weight is a "
458 "dependent variable (and hence we cannot take its "
459 "derivative directly)."));
460 Assert(row_n != dependent_direction,
462 "We cannot differentiate with respect to the variable "
463 "that is assumed to be dependent."));
464
465 const Point<spacedim> manifold_point = f(center);
466 const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
467 {row_n, dependent_direction}, center, step, f);
468 double entry = 0.0;
469 for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
470 entry +=
471 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
472 return entry;
473 }
474
480 template <typename Iterator, int spacedim, int structdim>
482 project_to_d_linear_object(const Iterator &object,
483 const Point<spacedim> &trial_point)
484 {
485 // let's look at this for simplicity for a quadrilateral
486 // (structdim==2) in a space with spacedim>2 (notate trial_point by
487 // y): all points on the surface are given by
488 // x(\xi) = sum_i v_i phi_x(\xi)
489 // where v_i are the vertices of the quadrilateral, and
490 // \xi=(\xi_1,\xi_2) are the reference coordinates of the
491 // quadrilateral. so what we are trying to do is find a point x on the
492 // surface that is closest to the point y. there are different ways to
493 // solve this problem, but in the end it's a nonlinear problem and we
494 // have to find reference coordinates \xi so that J(\xi) = 1/2 ||
495 // x(\xi)-y ||^2 is minimal. x(\xi) is a function that is
496 // structdim-linear in \xi, so J(\xi) is a polynomial of degree
497 // 2*structdim that we'd like to minimize. unless structdim==1, we'll
498 // have to use a Newton method to find the answer. This leads to the
499 // following formulation of Newton steps:
500 //
501 // Given \xi_k, find \delta\xi_k so that
502 // H_k \delta\xi_k = - F_k
503 // where H_k is an approximation to the second derivatives of J at
504 // \xi_k, and F_k is the first derivative of J. We'll iterate this a
505 // number of times until the right hand side is small enough. As a
506 // stopping criterion, we terminate if ||\delta\xi||<eps.
507 //
508 // As for the Hessian, the best choice would be
509 // H_k = J''(\xi_k)
510 // but we'll opt for the simpler Gauss-Newton form
511 // H_k = A^T A
512 // i.e.
513 // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
514 // \partial_n phi_i *
515 // \partial_m phi_j
516 // we start at xi=(0.5, 0.5).
518 for (unsigned int d = 0; d < structdim; ++d)
519 xi[d] = 0.5;
520
521 Point<spacedim> x_k;
522 for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
523 x_k += object->vertex(i) *
524 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
525
526 do
527 {
529 for (const unsigned int i :
530 GeometryInfo<structdim>::vertex_indices())
531 F_k +=
532 (x_k - trial_point) * object->vertex(i) *
533 GeometryInfo<structdim>::d_linear_shape_function_gradient(xi,
534 i);
535
537 for (const unsigned int i :
538 GeometryInfo<structdim>::vertex_indices())
539 for (const unsigned int j :
540 GeometryInfo<structdim>::vertex_indices())
541 {
544 xi, i),
546 xi, j));
547 H_k += (object->vertex(i) * object->vertex(j)) * tmp;
548 }
549
550 const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
551 xi += delta_xi;
552
553 x_k = Point<spacedim>();
554 for (const unsigned int i :
555 GeometryInfo<structdim>::vertex_indices())
556 x_k += object->vertex(i) *
557 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
558
559 if (delta_xi.norm() < 1e-7)
560 break;
561 }
562 while (true);
563
564 return x_k;
565 }
566 } // namespace ProjectToObject
567
568 // We hit an internal compiler error in ICC 15 if we define this as a lambda
569 // inside the project_to_object function below.
570 template <int structdim>
571 inline bool
572 weights_are_ok(
574 {
575 // clang has trouble figuring out structdim here, so define it
576 // again:
577 static const std::size_t n_vertices_per_cell =
579 n_independent_components;
580 std::array<double, n_vertices_per_cell> copied_weights;
581 for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
582 {
583 copied_weights[i] = v[i];
584 if (v[i] < 0.0 || v[i] > 1.0)
585 return false;
586 }
587
588 // check the sum: try to avoid some roundoff errors by summing in order
589 std::sort(copied_weights.begin(), copied_weights.end());
590 const double sum =
591 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
592 return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
593 }
594 } // namespace internal
595
596
597
598 template <typename Iterator>
601 const Iterator &object,
603 {
604 const int spacedim = Iterator::AccessorType::space_dimension;
605 const int structdim = Iterator::AccessorType::structure_dimension;
606
607 Point<spacedim> projected_point = trial_point;
608
609 if (structdim >= spacedim)
610 return projected_point;
611 else if (structdim == 1 || structdim == 2)
612 {
613 using namespace internal::ProjectToObject;
614 // Try to use the special flat algorithm for quads (this is better
615 // than the general algorithm in 3d). This does not take into account
616 // whether projected_point is outside the quad, but we optimize along
617 // lines below anyway:
618 const int dim = Iterator::AccessorType::dimension;
619 const Manifold<dim, spacedim> &manifold = object->get_manifold();
620 if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
621 &manifold) != nullptr)
622 {
623 projected_point =
624 project_to_d_linear_object<Iterator, spacedim, structdim>(
625 object, trial_point);
626 }
627 else
628 {
629 // We want to find a point on the convex hull (defined by the
630 // vertices of the object and the manifold description) that is
631 // relatively close to the trial point. This has a few issues:
632 //
633 // 1. For a general convex hull we are not guaranteed that a unique
634 // minimum exists.
635 // 2. The independent variables in the optimization process are the
636 // weights given to Manifold::get_new_point, which must sum to 1,
637 // so we cannot use standard finite differences to approximate a
638 // gradient.
639 //
640 // There is not much we can do about 1., but for 2. we can derive
641 // finite difference stencils that work on a structdim-dimensional
642 // simplex and rewrite the optimization problem to use those
643 // instead. Consider the structdim 2 case and let
644 //
645 // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
646 // c2, c3})
647 //
648 // where {c0, c1, c2, c3} are the weights for the four vertices on
649 // the quadrilateral. We seek to minimize the Euclidean distance
650 // between F(...) and trial_point. We can solve for c3 in terms of
651 // the other weights and get, for one coordinate direction
652 //
653 // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
654 // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
655 //
656 // where we substitute back in for c3 after taking the
657 // derivative. We can compute a stencil for the cross derivative
658 // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
659 // (and gradient_entry computes the sum over the independent
660 // variables). Below, we somewhat arbitrarily pick the last
661 // component as the dependent one.
662 //
663 // Since we can now calculate derivatives of the objective
664 // function we can use gradient descent to minimize it.
665 //
666 // Of course, this is much simpler in the structdim = 1 case (we
667 // could rewrite the projection as a 1d optimization problem), but
668 // to reduce the potential for bugs we use the same code in both
669 // cases.
670 const double step_size = object->diameter() / 64.0;
671
672 constexpr unsigned int n_vertices_per_cell =
674
675 std::array<Point<spacedim>, n_vertices_per_cell> vertices;
676 for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
677 ++vertex_n)
678 vertices[vertex_n] = object->vertex(vertex_n);
679
680 auto get_point_from_weights =
681 [&](const Tensor<1, n_vertices_per_cell> &weights)
682 -> Point<spacedim> {
683 return object->get_manifold().get_new_point(
684 make_array_view(vertices.begin(), vertices.end()),
685 make_array_view(weights.begin_raw(), weights.end_raw()));
686 };
687
688 // pick the initial weights as (normalized) inverse distances from
689 // the trial point:
690 Tensor<1, n_vertices_per_cell> guess_weights;
691 double guess_weights_sum = 0.0;
692 for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
693 ++vertex_n)
694 {
695 const double distance =
696 vertices[vertex_n].distance(trial_point);
697 if (distance == 0.0)
698 {
699 guess_weights = 0.0;
700 guess_weights[vertex_n] = 1.0;
701 guess_weights_sum = 1.0;
702 break;
703 }
704 else
705 {
706 guess_weights[vertex_n] = 1.0 / distance;
707 guess_weights_sum += guess_weights[vertex_n];
708 }
709 }
710 guess_weights /= guess_weights_sum;
711 Assert(internal::weights_are_ok<structdim>(guess_weights),
713
714 // The optimization algorithm consists of two parts:
715 //
716 // 1. An outer loop where we apply the gradient descent algorithm.
717 // 2. An inner loop where we do a line search to find the optimal
718 // length of the step one should take in the gradient direction.
719 //
720 for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
721 {
722 const unsigned int dependent_direction =
723 n_vertices_per_cell - 1;
724 Tensor<1, n_vertices_per_cell> current_gradient;
725 for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
726 ++row_n)
727 {
728 if (row_n != dependent_direction)
729 {
730 current_gradient[row_n] =
731 gradient_entry<spacedim, structdim>(
732 row_n,
733 dependent_direction,
734 trial_point,
735 guess_weights,
736 step_size,
737 get_point_from_weights);
738
739 current_gradient[dependent_direction] -=
740 current_gradient[row_n];
741 }
742 }
743
744 // We need to travel in the -gradient direction, as noted
745 // above, but we may not want to take a full step in that
746 // direction; instead, guess that we will go -0.5*gradient and
747 // do quasi-Newton iteration to pick the best multiplier. The
748 // goal is to find a scalar alpha such that
749 //
750 // F(x - alpha g)
751 //
752 // is minimized, where g is the gradient and F is the
753 // objective function. To find the optimal value we find roots
754 // of the derivative of the objective function with respect to
755 // alpha by Newton iteration, where we approximate the first
756 // and second derivatives of F(x - alpha g) with centered
757 // finite differences.
758 double gradient_weight = -0.5;
759 auto gradient_weight_objective_function =
760 [&](const double gradient_weight_guess) -> double {
761 return (trial_point -
762 get_point_from_weights(guess_weights +
763 gradient_weight_guess *
764 current_gradient))
765 .norm_square();
766 };
767
768 for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
769 {
770 const double update_numerator = centered_first_difference(
771 gradient_weight,
772 step_size,
773 gradient_weight_objective_function);
774 const double update_denominator =
775 centered_second_difference(
776 gradient_weight,
777 step_size,
778 gradient_weight_objective_function);
779
780 // avoid division by zero. Note that we limit the gradient
781 // weight below
782 if (std::abs(update_denominator) == 0.0)
783 break;
784 gradient_weight =
785 gradient_weight - update_numerator / update_denominator;
786
787 // Put a fairly lenient bound on the largest possible
788 // gradient (things tend to be locally flat, so the gradient
789 // itself is usually small)
790 if (std::abs(gradient_weight) > 10)
791 {
792 gradient_weight = -10.0;
793 break;
794 }
795 }
796
797 // It only makes sense to take convex combinations with weights
798 // between zero and one. If the update takes us outside of this
799 // region then rescale the update to stay within the region and
800 // try again
801 Tensor<1, n_vertices_per_cell> tentative_weights =
802 guess_weights + gradient_weight * current_gradient;
803
804 double new_gradient_weight = gradient_weight;
805 for (unsigned int iteration_count = 0; iteration_count < 40;
806 ++iteration_count)
807 {
808 if (internal::weights_are_ok<structdim>(tentative_weights))
809 break;
810
811 for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
812 {
813 if (tentative_weights[i] < 0.0)
814 {
815 tentative_weights -=
816 (tentative_weights[i] / current_gradient[i]) *
817 current_gradient;
818 }
819 if (tentative_weights[i] < 0.0 ||
820 1.0 < tentative_weights[i])
821 {
822 new_gradient_weight /= 2.0;
823 tentative_weights =
824 guess_weights +
825 new_gradient_weight * current_gradient;
826 }
827 }
828 }
829
830 // the update might still send us outside the valid region, so
831 // check again and quit if the update is still not valid
832 if (!internal::weights_are_ok<structdim>(tentative_weights))
833 break;
834
835 // if we cannot get closer by traveling in the gradient
836 // direction then quit
837 if (get_point_from_weights(tentative_weights)
838 .distance(trial_point) <
839 get_point_from_weights(guess_weights).distance(trial_point))
840 guess_weights = tentative_weights;
841 else
842 break;
843 Assert(internal::weights_are_ok<structdim>(guess_weights),
845 }
846 Assert(internal::weights_are_ok<structdim>(guess_weights),
848 projected_point = get_point_from_weights(guess_weights);
849 }
850
851 // if structdim == 2 and the optimal point is not on the interior then
852 // we may be able to get a more accurate result by projecting onto the
853 // lines.
854 if (structdim == 2)
855 {
856 std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
857 line_projections;
858 for (unsigned int line_n = 0;
859 line_n < GeometryInfo<structdim>::lines_per_cell;
860 ++line_n)
861 {
862 line_projections[line_n] =
863 project_to_object(object->line(line_n), trial_point);
864 }
865 std::sort(line_projections.begin(),
866 line_projections.end(),
867 [&](const Point<spacedim> &a, const Point<spacedim> &b) {
868 return a.distance(trial_point) <
869 b.distance(trial_point);
870 });
871 if (line_projections[0].distance(trial_point) <
872 projected_point.distance(trial_point))
873 projected_point = line_projections[0];
874 }
875 }
876 else
877 {
879 return projected_point;
880 }
881
882 return projected_point;
883 }
884} // namespace GridTools
885#endif // DOXYGEN
886
888
889#endif
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:175
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
unsigned int vertex_indices[2]
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_NOT_IMPLEMENTED()
std::pair< DerivativeForm< 1, dim, spacedim >, Tensor< 1, spacedim > > affine_cell_approximation(const ArrayView< const Point< spacedim > > &vertices)
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Point< Iterator::AccessorType::space_dimension > project_to_object(const Iterator &object, const Point< Iterator::AccessorType::space_dimension > &trial_point)
Vector< double > compute_aspect_ratio_of_cells(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
double compute_maximum_aspect_ratio(const Mapping< dim > &mapping, const Triangulation< dim > &triangulation, const Quadrature< dim > &quadrature)
double volume(const Triangulation< dim, spacedim > &tria)
double diameter(const Triangulation< dim, spacedim > &tria)
BoundingBox< spacedim > compute_bounding_box(const Triangulation< dim, spacedim > &triangulation)
double maximal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
std::pair< unsigned int, double > get_longest_direction(typename Triangulation< dim, spacedim >::active_cell_iterator cell)
double cell_measure(const std::vector< Point< dim > > &all_vertices, const ArrayView< const unsigned int > &vertex_indices)
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)