Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2684-gc61376a70f 2025-02-22 15: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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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 <numeric>
34#include <utility>
35#include <vector>
36
38
39namespace GridTools
40{
52 template <int dim, int spacedim>
53 double
55
79 template <int dim, int spacedim>
80 double
82
110 template <int dim, int spacedim>
111 double
113 const Mapping<dim, spacedim> &mapping);
114
125 template <int dim, int spacedim>
126 double
129 const Mapping<dim, spacedim> &mapping =
130 (ReferenceCells::get_hypercube<dim>()
131#ifndef _MSC_VER
132 .template get_default_linear_mapping<dim, spacedim>()
133#else
134 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
135#endif
136 ));
137
148 template <int dim, int spacedim>
149 double
152 const Mapping<dim, spacedim> &mapping =
153 (ReferenceCells::get_hypercube<dim>()
154#ifndef _MSC_VER
155 .template get_default_linear_mapping<dim, spacedim>()
156#else
157 .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
158#endif
159 ));
160
180 template <int dim>
181 double
182 cell_measure(const std::vector<Point<dim>> &all_vertices,
184
196 template <int dim, int spacedim>
197 std::pair<unsigned int, double>
200
223 template <int dim, int spacedim>
224 std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
226
256 template <int dim>
260 const Quadrature<dim> &quadrature);
261
269 template <int dim>
270 double
273 const Quadrature<dim> &quadrature);
274
288 template <int dim, int spacedim>
291
309 template <typename MeshType>
311 std::pair<
313 Point<MeshType::
314 space_dimension>> compute_bounding_box(const MeshType &mesh,
315 const std::function<bool(
316 const typename MeshType::
317 active_cell_iterator &)>
318 &predicate);
319
337 template <typename Iterator>
340 const Iterator &object,
343} // namespace GridTools
344
345#ifndef DOXYGEN
346namespace GridTools
347{
348 namespace internal
349 {
350 namespace ProjectToObject
351 {
364 struct CrossDerivative
365 {
366 const unsigned int direction_0;
367 const unsigned int direction_1;
368
369 CrossDerivative(const unsigned int d0, const unsigned int d1);
370 };
371
372 inline CrossDerivative::CrossDerivative(const unsigned int d0,
373 const unsigned int d1)
374 : direction_0(d0)
375 , direction_1(d1)
376 {}
377
378
379
384 template <typename F>
385 inline auto
386 centered_first_difference(const double center,
387 const double step,
388 const F &f) -> decltype(f(center) - f(center))
389 {
390 return (f(center + step) - f(center - step)) / (2.0 * step);
391 }
392
393
394
399 template <typename F>
400 inline auto
401 centered_second_difference(const double center,
402 const double step,
403 const F &f) -> decltype(f(center) - f(center))
404 {
405 return (f(center + step) - 2.0 * f(center) + f(center - step)) /
406 (step * step);
407 }
408
409
410
420 template <int structdim, typename F>
421 inline auto
422 cross_stencil(
423 const CrossDerivative cross_derivative,
425 const double step,
426 const F &f) -> decltype(f(center) - f(center))
427 {
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)) /
434 step;
435 }
436
437
438
445 template <int spacedim, int structdim, typename F>
446 inline double
447 gradient_entry(
448 const unsigned int row_n,
449 const unsigned int dependent_direction,
450 const Point<spacedim> &p0,
452 const double step,
453 const F &f)
454 {
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."));
465
466 const Point<spacedim> manifold_point = f(center);
467 const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>(
468 {row_n, dependent_direction}, center, step, f);
469 double entry = 0.0;
470 for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
471 entry +=
472 -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
473 return entry;
474 }
475
481 template <typename Iterator, int spacedim, int structdim>
483 project_to_d_linear_object(const Iterator &object,
484 const Point<spacedim> &trial_point)
485 {
486 // let's look at this for simplicity for a quadrilateral
487 // (structdim==2) in a space with spacedim>2 (notate trial_point by
488 // y): all points on the surface are given by
489 // x(\xi) = sum_i v_i phi_x(\xi)
490 // where v_i are the vertices of the quadrilateral, and
491 // \xi=(\xi_1,\xi_2) are the reference coordinates of the
492 // quadrilateral. so what we are trying to do is find a point x on the
493 // surface that is closest to the point y. there are different ways to
494 // solve this problem, but in the end it's a nonlinear problem and we
495 // have to find reference coordinates \xi so that J(\xi) = 1/2 ||
496 // x(\xi)-y ||^2 is minimal. x(\xi) is a function that is
497 // structdim-linear in \xi, so J(\xi) is a polynomial of degree
498 // 2*structdim that we'd like to minimize. unless structdim==1, we'll
499 // have to use a Newton method to find the answer. This leads to the
500 // following formulation of Newton steps:
501 //
502 // Given \xi_k, find \delta\xi_k so that
503 // H_k \delta\xi_k = - F_k
504 // where H_k is an approximation to the second derivatives of J at
505 // \xi_k, and F_k is the first derivative of J. We'll iterate this a
506 // number of times until the right hand side is small enough. As a
507 // stopping criterion, we terminate if ||\delta\xi||<eps.
508 //
509 // As for the Hessian, the best choice would be
510 // H_k = J''(\xi_k)
511 // but we'll opt for the simpler Gauss-Newton form
512 // H_k = A^T A
513 // i.e.
514 // (H_k)_{nm} = \sum_{i,j} v_i*v_j *
515 // \partial_n phi_i *
516 // \partial_m phi_j
517 // we start at xi=(0.5, 0.5).
519 for (unsigned int d = 0; d < structdim; ++d)
520 xi[d] = 0.5;
521
522 Point<spacedim> x_k;
523 for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
524 x_k += object->vertex(i) *
525 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
526
527 do
528 {
530 for (const unsigned int i :
531 GeometryInfo<structdim>::vertex_indices())
532 F_k +=
533 (x_k - trial_point) * object->vertex(i) *
534 GeometryInfo<structdim>::d_linear_shape_function_gradient(xi,
535 i);
536
538 for (const unsigned int i :
539 GeometryInfo<structdim>::vertex_indices())
540 for (const unsigned int j :
541 GeometryInfo<structdim>::vertex_indices())
542 {
545 xi, i),
547 xi, j));
548 H_k += (object->vertex(i) * object->vertex(j)) * tmp;
549 }
550
551 const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
552 xi += delta_xi;
553
554 x_k = Point<spacedim>();
555 for (const unsigned int i :
556 GeometryInfo<structdim>::vertex_indices())
557 x_k += object->vertex(i) *
558 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
559
560 if (delta_xi.norm() < 1e-7)
561 break;
562 }
563 while (true);
564
565 return x_k;
566 }
567 } // namespace ProjectToObject
568
569 // We hit an internal compiler error in ICC 15 if we define this as a lambda
570 // inside the project_to_object function below.
571 template <int structdim>
572 inline bool
573 weights_are_ok(
575 {
576 // clang has trouble figuring out structdim here, so define it
577 // again:
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)
583 {
584 copied_weights[i] = v[i];
585 if (v[i] < 0.0 || v[i] > 1.0)
586 return false;
587 }
588
589 // check the sum: try to avoid some roundoff errors by summing in order
590 std::sort(copied_weights.begin(), copied_weights.end());
591 const double sum =
592 std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
593 return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
594 }
595 } // namespace internal
596
597
598
599 template <typename Iterator>
602 const Iterator &object,
604 {
605 const int spacedim = Iterator::AccessorType::space_dimension;
606 const int structdim = Iterator::AccessorType::structure_dimension;
607
608 Point<spacedim> projected_point = trial_point;
609
610 if (structdim >= spacedim)
611 return projected_point;
612 else if (structdim == 1 || structdim == 2)
613 {
614 using namespace internal::ProjectToObject;
615 // Try to use the special flat algorithm for quads (this is better
616 // than the general algorithm in 3d). This does not take into account
617 // whether projected_point is outside the quad, but we optimize along
618 // lines below anyway:
619 const int dim = Iterator::AccessorType::dimension;
620 const Manifold<dim, spacedim> &manifold = object->get_manifold();
621 if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
622 &manifold) != nullptr)
623 {
624 projected_point =
625 project_to_d_linear_object<Iterator, spacedim, structdim>(
626 object, trial_point);
627 }
628 else
629 {
630 // We want to find a point on the convex hull (defined by the
631 // vertices of the object and the manifold description) that is
632 // relatively close to the trial point. This has a few issues:
633 //
634 // 1. For a general convex hull we are not guaranteed that a unique
635 // minimum exists.
636 // 2. The independent variables in the optimization process are the
637 // weights given to Manifold::get_new_point, which must sum to 1,
638 // so we cannot use standard finite differences to approximate a
639 // gradient.
640 //
641 // There is not much we can do about 1., but for 2. we can derive
642 // finite difference stencils that work on a structdim-dimensional
643 // simplex and rewrite the optimization problem to use those
644 // instead. Consider the structdim 2 case and let
645 //
646 // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
647 // c2, c3})
648 //
649 // where {c0, c1, c2, c3} are the weights for the four vertices on
650 // the quadrilateral. We seek to minimize the Euclidean distance
651 // between F(...) and trial_point. We can solve for c3 in terms of
652 // the other weights and get, for one coordinate direction
653 //
654 // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
655 // = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
656 //
657 // where we substitute back in for c3 after taking the
658 // derivative. We can compute a stencil for the cross derivative
659 // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
660 // (and gradient_entry computes the sum over the independent
661 // variables). Below, we somewhat arbitrarily pick the last
662 // component as the dependent one.
663 //
664 // Since we can now calculate derivatives of the objective
665 // function we can use gradient descent to minimize it.
666 //
667 // Of course, this is much simpler in the structdim = 1 case (we
668 // could rewrite the projection as a 1d optimization problem), but
669 // to reduce the potential for bugs we use the same code in both
670 // cases.
671 const double step_size = object->diameter() / 64.0;
672
673 constexpr unsigned int n_vertices_per_cell =
675
676 std::array<Point<spacedim>, n_vertices_per_cell> vertices;
677 for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
678 ++vertex_n)
679 vertices[vertex_n] = object->vertex(vertex_n);
680
681 auto get_point_from_weights =
682 [&](const Tensor<1, n_vertices_per_cell> &weights)
683 -> Point<spacedim> {
684 return object->get_manifold().get_new_point(
685 make_array_view(vertices.begin(), vertices.end()),
686 make_array_view(weights.begin_raw(), weights.end_raw()));
687 };
688
689 // pick the initial weights as (normalized) inverse distances from
690 // the trial point:
691 Tensor<1, n_vertices_per_cell> guess_weights;
692 double guess_weights_sum = 0.0;
693 for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
694 ++vertex_n)
695 {
696 const double distance =
697 vertices[vertex_n].distance(trial_point);
698 if (distance == 0.0)
699 {
700 guess_weights = 0.0;
701 guess_weights[vertex_n] = 1.0;
702 guess_weights_sum = 1.0;
703 break;
704 }
705 else
706 {
707 guess_weights[vertex_n] = 1.0 / distance;
708 guess_weights_sum += guess_weights[vertex_n];
709 }
710 }
711 guess_weights /= guess_weights_sum;
712 Assert(internal::weights_are_ok<structdim>(guess_weights),
714
715 // The optimization algorithm consists of two parts:
716 //
717 // 1. An outer loop where we apply the gradient descent algorithm.
718 // 2. An inner loop where we do a line search to find the optimal
719 // length of the step one should take in the gradient direction.
720 //
721 for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
722 {
723 const unsigned int dependent_direction =
724 n_vertices_per_cell - 1;
725 Tensor<1, n_vertices_per_cell> current_gradient;
726 for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
727 ++row_n)
728 {
729 if (row_n != dependent_direction)
730 {
731 current_gradient[row_n] =
732 gradient_entry<spacedim, structdim>(
733 row_n,
734 dependent_direction,
735 trial_point,
736 guess_weights,
737 step_size,
738 get_point_from_weights);
739
740 current_gradient[dependent_direction] -=
741 current_gradient[row_n];
742 }
743 }
744
745 // We need to travel in the -gradient direction, as noted
746 // above, but we may not want to take a full step in that
747 // direction; instead, guess that we will go -0.5*gradient and
748 // do quasi-Newton iteration to pick the best multiplier. The
749 // goal is to find a scalar alpha such that
750 //
751 // F(x - alpha g)
752 //
753 // is minimized, where g is the gradient and F is the
754 // objective function. To find the optimal value we find roots
755 // of the derivative of the objective function with respect to
756 // alpha by Newton iteration, where we approximate the first
757 // and second derivatives of F(x - alpha g) with centered
758 // finite differences.
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 *
765 current_gradient))
766 .norm_square();
767 };
768
769 for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
770 {
771 const double update_numerator = centered_first_difference(
772 gradient_weight,
773 step_size,
774 gradient_weight_objective_function);
775 const double update_denominator =
776 centered_second_difference(
777 gradient_weight,
778 step_size,
779 gradient_weight_objective_function);
780
781 // avoid division by zero. Note that we limit the gradient
782 // weight below
783 if (std::abs(update_denominator) == 0.0)
784 break;
785 gradient_weight =
786 gradient_weight - update_numerator / update_denominator;
787
788 // Put a fairly lenient bound on the largest possible
789 // gradient (things tend to be locally flat, so the gradient
790 // itself is usually small)
791 if (std::abs(gradient_weight) > 10)
792 {
793 gradient_weight = -10.0;
794 break;
795 }
796 }
797
798 // It only makes sense to take convex combinations with weights
799 // between zero and one. If the update takes us outside of this
800 // region then rescale the update to stay within the region and
801 // try again
802 Tensor<1, n_vertices_per_cell> tentative_weights =
803 guess_weights + gradient_weight * current_gradient;
804
805 double new_gradient_weight = gradient_weight;
806 for (unsigned int iteration_count = 0; iteration_count < 40;
807 ++iteration_count)
808 {
809 if (internal::weights_are_ok<structdim>(tentative_weights))
810 break;
811
812 for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
813 {
814 if (tentative_weights[i] < 0.0)
815 {
816 tentative_weights -=
817 (tentative_weights[i] / current_gradient[i]) *
818 current_gradient;
819 }
820 if (tentative_weights[i] < 0.0 ||
821 1.0 < tentative_weights[i])
822 {
823 new_gradient_weight /= 2.0;
824 tentative_weights =
825 guess_weights +
826 new_gradient_weight * current_gradient;
827 }
828 }
829 }
830
831 // the update might still send us outside the valid region, so
832 // check again and quit if the update is still not valid
833 if (!internal::weights_are_ok<structdim>(tentative_weights))
834 break;
835
836 // if we cannot get closer by traveling in the gradient
837 // direction then quit
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;
842 else
843 break;
844 Assert(internal::weights_are_ok<structdim>(guess_weights),
846 }
847 Assert(internal::weights_are_ok<structdim>(guess_weights),
849 projected_point = get_point_from_weights(guess_weights);
850 }
851
852 // if structdim == 2 and the optimal point is not on the interior then
853 // we may be able to get a more accurate result by projecting onto the
854 // lines.
855 if (structdim == 2)
856 {
857 std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
858 line_projections;
859 for (unsigned int line_n = 0;
860 line_n < GeometryInfo<structdim>::lines_per_cell;
861 ++line_n)
862 {
863 line_projections[line_n] =
864 project_to_object(object->line(line_n), trial_point);
865 }
866 std::sort(line_projections.begin(),
867 line_projections.end(),
868 [&](const Point<spacedim> &a, const Point<spacedim> &b) {
869 return a.distance(trial_point) <
870 b.distance(trial_point);
871 });
872 if (line_projections[0].distance(trial_point) <
873 projected_point.distance(trial_point))
874 projected_point = line_projections[0];
875 }
876 }
877 else
878 {
880 return projected_point;
881 }
882
883 return projected_point;
884 }
885} // namespace GridTools
886#endif // DOXYGEN
887
889
890#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:320
Definition point.h:113
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:522
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:194
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:523
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int vertex_indices[2]
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
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)