]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move geometric GridTools functions to a new file. 16149/head
authorDavid Wells <drwells@email.unc.edu>
Sun, 15 Oct 2023 22:30:58 +0000 (18:30 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sun, 15 Oct 2023 23:05:01 +0000 (19:05 -0400)
grid_tools.cc is presently our most expensive file to compile. This is the first
commit in a sequence which will chop it up in the same way we chopped up
VectorTools.

include/deal.II/fe/mapping_q_internal.h
include/deal.II/grid/grid_tools.h
include/deal.II/grid/grid_tools_geometry.h [new file with mode: 0644]
source/grid/CMakeLists.txt
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
source/grid/grid_tools_geometry.cc [new file with mode: 0644]
source/grid/grid_tools_geometry.inst.in [new file with mode: 0644]
source/grid/tria_accessor.cc

index 80250bbdd76c63516a69b50803b6ec39e870978c..dd56d422f0a4d4e204a40e991b12d35fd707c3bc 100644 (file)
@@ -31,7 +31,7 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q.h>
 
-#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_geometry.h>
 
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
index 7faf8c8afc3b87349af2f87fc4a648ef3cbffd76..b42ec6043f20877d7b89b8dda5ea6d60cf11bc3d 100644 (file)
@@ -32,6 +32,7 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping.h>
 
+#include <deal.II/grid/grid_tools_geometry.h>
 #include <deal.II/grid/manifold.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_accessor.h>
@@ -132,286 +133,6 @@ namespace internal
  */
 namespace GridTools
 {
-  /**
-   * @name Information about meshes and cells
-   */
-  /** @{ */
-
-  /**
-   * Return the diameter of a triangulation. The diameter is computed using
-   * only the vertices, i.e. if the diameter should be larger than the maximal
-   * distance between boundary vertices due to a higher order mapping, then
-   * this function will not catch this.
-   */
-  template <int dim, int spacedim>
-  double
-  diameter(const Triangulation<dim, spacedim> &tria);
-
-  /**
-   * Compute the volume (i.e. the dim-dimensional measure) of the
-   * triangulation. We compute the measure using the integral $\sum_K \int_K 1
-   * \; dx$ where $K$ are the cells of the given triangulation. The integral
-   * is approximated via quadrature. This version of the function uses a
-   * linear mapping to compute the JxW values on each cell.
-   *
-   * If the triangulation is a dim-dimensional one embedded in a higher
-   * dimensional space of dimension spacedim, then the value returned is the
-   * dim-dimensional measure. For example, for a two-dimensional triangulation
-   * in three-dimensional space, the value returned is the area of the surface
-   * so described. (This obviously makes sense since the spacedim-dimensional
-   * measure of a dim-dimensional triangulation would always be zero if dim @<
-   * spacedim).
-   *
-   * This function also works for objects of type
-   * parallel::distributed::Triangulation, in which case the function is a
-   * collective operation.
-   *
-   * @param tria The triangulation.
-   * @return The dim-dimensional measure of the domain described by the
-   * triangulation, as discussed above.
-   */
-  template <int dim, int spacedim>
-  double
-  volume(const Triangulation<dim, spacedim> &tria);
-
-  /**
-   * Compute the volume (i.e. the dim-dimensional measure) of the
-   * triangulation. We compute the measure using the integral $\sum_K \int_K 1
-   * \; dx$ where $K$ are the cells of the given triangulation. The integral
-   * is approximated via quadrature for which we use the mapping argument.
-   *
-   * If the triangulation is a dim-dimensional one embedded in a higher
-   * dimensional space of dimension spacedim, then the value returned is the
-   * dim-dimensional measure. For example, for a two-dimensional triangulation
-   * in three-dimensional space, the value returned is the area of the surface
-   * so described. (This obviously makes sense since the spacedim-dimensional
-   * measure of a dim-dimensional triangulation would always be zero if dim @<
-   * spacedim.
-   *
-   * This function also works for objects of type
-   * parallel::distributed::Triangulation, in which case the function is a
-   * collective operation.
-   *
-   * @param tria The triangulation.
-   * @param mapping The Mapping which computes the Jacobians used to
-   * approximate the volume via quadrature. Explicitly using a higher-order
-   * Mapping (i.e., instead of using the other version of this function) will
-   * result in a more accurate approximation of the volume on Triangulations
-   * with curvature described by Manifold objects.
-   * @return The dim-dimensional measure of the domain described by the
-   * triangulation, as discussed above.
-   */
-  template <int dim, int spacedim>
-  double
-  volume(const Triangulation<dim, spacedim> &tria,
-         const Mapping<dim, spacedim>       &mapping);
-
-  /**
-   * Return an approximation of the diameter of the smallest active cell of a
-   * triangulation. See step-24 for an example of use of this function.
-   *
-   * Notice that, even if you pass a non-trivial mapping, the returned value is
-   * computed only using information on the vertices of the triangulation,
-   * possibly transformed by the mapping. While this is accurate most of the
-   * times, it may fail to give the correct result when the triangulation
-   * contains very distorted cells.
-   */
-  template <int dim, int spacedim>
-  double
-  minimal_cell_diameter(
-    const Triangulation<dim, spacedim> &triangulation,
-    const Mapping<dim, spacedim>       &mapping =
-      (ReferenceCells::get_hypercube<dim>()
-#ifndef _MSC_VER
-         .template get_default_linear_mapping<dim, spacedim>()
-#else
-         .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
-#endif
-         ));
-
-  /**
-   * Return an approximation of the diameter of the largest active cell of a
-   * triangulation.
-   *
-   * Notice that, even if you pass a non-trivial mapping to this function, the
-   * returned value is computed only using information on the vertices of the
-   * triangulation, possibly transformed by the mapping. While this is accurate
-   * most of the times, it may fail to give the correct result when the
-   * triangulation contains very distorted cells.
-   */
-  template <int dim, int spacedim>
-  double
-  maximal_cell_diameter(
-    const Triangulation<dim, spacedim> &triangulation,
-    const Mapping<dim, spacedim>       &mapping =
-      (ReferenceCells::get_hypercube<dim>()
-#ifndef _MSC_VER
-         .template get_default_linear_mapping<dim, spacedim>()
-#else
-         .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
-#endif
-         ));
-
-  /**
-   * Given a list of vertices (typically obtained using
-   * Triangulation::get_vertices()) as the first, and a list of vertex indices
-   * that characterize a single cell as the second argument, return the
-   * measure (area, volume) of this cell. If this is a real cell, then you can
-   * get the same result using <code>cell-@>measure()</code>, but this
-   * function also works for cells that do not exist except that you make it
-   * up by naming its vertices from the list.
-   *
-   * The size of @p vertex_indices, combined with `dim`, implicitly encodes
-   * the ReferenceCell type of the provided cell. For example, if `dim == 2` and
-   * `vertex_indices.size() == 3` then the cell is a triangle, but if
-   * `dim == 2` and `vertex_indices.size() == 4` then the cell is a
-   * quadrilateral. A std::vector is implicitly convertible to an ArrayView, so
-   * it can be passed directly to this function. See the ArrayView class for
-   * more information.
-   *
-   * @note This function is only implemented for codimension zero objects.
-   */
-  template <int dim>
-  double
-  cell_measure(const std::vector<Point<dim>>       &all_vertices,
-               const ArrayView<const unsigned int> &vertex_indices);
-
-  /**
-   * This function computes an affine approximation of the map from the unit
-   * coordinates to the real coordinates of the form $p_\text{real} = A
-   * p_\text{unit} + b $ by a least squares fit of this affine function to the
-   * $2^\text{dim}$ vertices representing a quadrilateral or hexahedral cell
-   * in `spacedim` dimensions. The result is returned as a pair with the
-   * matrix <i>A</i> as the first argument and the vector <i>b</i> describing
-   * distance of the plane to the origin.
-   *
-   * For any valid mesh cell whose geometry is not degenerate, this operation
-   * results in a unique affine mapping, even in cases where the actual
-   * transformation by a bi-/trilinear or higher order mapping might be
-   * singular. The result is exact in case the transformation from the unit to
-   * the real cell is indeed affine, such as in one dimension or for Cartesian
-   * and affine (parallelogram) meshes in 2d/3d.
-   *
-   * This approximation is underlying the function
-   * TriaAccessor::real_to_unit_cell_affine_approximation() function.
-   *
-   * For exact transformations to the unit cell, use
-   * Mapping::transform_real_to_unit_cell().
-   */
-  template <int dim, int spacedim>
-  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
-  affine_cell_approximation(const ArrayView<const Point<spacedim>> &vertices);
-
-  /**
-   * Computes an aspect ratio measure for all locally-owned active cells and
-   * fills a vector with one entry per cell, given a @p triangulation and
-   * @p mapping. The size of the vector that is returned equals the number of
-   * active cells. The vector contains zero for non locally-owned cells. The
-   * aspect ratio of a cell is defined as the ratio of the maximum to minimum
-   * singular value of the Jacobian, taking the maximum over all quadrature
-   * points of a quadrature rule specified via @p quadrature. For example, for
-   * the special case of rectangular elements in 2d with dimensions $a$ and $b$
-   * ($a \geq b$), this function returns the usual aspect ratio definition
-   * $a/b$. The above definition using singular values is a generalization to
-   * arbitrarily deformed elements. This function is intended to be used for
-   * $d=2,3$ space dimensions, but it can also be used for $d=1$ returning a
-   * value of 1.
-   *
-   * @note Inverted elements do not throw an exception. Instead, a value of inf
-   * is written into the vector in case of inverted elements.
-   *
-   * @note Make sure to use enough quadrature points for a precise calculation
-   * of the aspect ratio in case of deformed elements.
-   *
-   * @note In parallel computations the return value will have the length
-   * n_active_cells but the aspect ratio is only computed for the cells that
-   * are locally owned and placed at index CellAccessor::active_cell_index(),
-   * respectively. All other values are set to 0.
-   *
-   * @note This function can only be used if deal.II was configured with
-   * support for LAPACK.
-   */
-  template <int dim>
-  Vector<double>
-  compute_aspect_ratio_of_cells(const Mapping<dim>       &mapping,
-                                const Triangulation<dim> &triangulation,
-                                const Quadrature<dim>    &quadrature);
-
-  /**
-   * Computes the maximum aspect ratio by taking the maximum over all cells.
-   *
-   * @note When running in parallel with a Triangulation that supports MPI,
-   * this is a collective call and the return value is the maximum over all
-   * processors.
-   */
-  template <int dim>
-  double
-  compute_maximum_aspect_ratio(const Mapping<dim>       &mapping,
-                               const Triangulation<dim> &triangulation,
-                               const Quadrature<dim>    &quadrature);
-
-  /**
-   * Compute the smallest box containing the entire triangulation.
-   *
-   * If the input triangulation is a `parallel::distributed::Triangulation`,
-   * then each processor will compute a bounding box enclosing all locally
-   * owned, ghost, and artificial cells. In the case of a domain without curved
-   * boundaries, these bounding boxes will all agree between processors because
-   * the union of the areas occupied by artificial and ghost cells equals the
-   * union of the areas occupied by the cells that other processors own.
-   * However, if the domain has curved boundaries, this is no longer the case.
-   * The bounding box returned may be appropriate for the current processor,
-   * but different from the bounding boxes computed on other processors.
-   */
-  template <int dim, int spacedim>
-  BoundingBox<spacedim>
-  compute_bounding_box(const Triangulation<dim, spacedim> &triangulation);
-
-  /**
-   * Return the point on the geometrical object @p object closest to the given
-   * point @p trial_point. For example, if @p object is a one-dimensional line
-   * or edge, then the returned point will be a point on the geodesic that
-   * connects the vertices as the manifold associated with the object sees it
-   * (i.e., the geometric line may be curved if it lives in a higher
-   * dimensional space). If the iterator points to a quadrilateral in a higher
-   * dimensional space, then the returned point lies within the convex hull of
-   * the vertices of the quad as seen by the associated manifold.
-   *
-   * @note This projection is usually not well-posed since there may be
-   * multiple points on the object that minimize the distance. The algorithm
-   * used in this function is robust (and the output is guaranteed to be on
-   * the given @p object) but may only provide a few correct digits if the
-   * object has high curvature. If your manifold supports it then the
-   * specialized function Manifold::project_to_manifold() may perform better.
-   */
-  template <typename Iterator>
-  Point<Iterator::AccessorType::space_dimension>
-  project_to_object(
-    const Iterator                                       &object,
-    const Point<Iterator::AccessorType::space_dimension> &trial_point);
-
-  /**
-   * Return the arrays that define the coarse mesh of a Triangulation. This
-   * function is the inverse of Triangulation::create_triangulation() in the
-   * sense that if one called this function on a triangulation, then that
-   * triangulation could be recreated by some kind of refinement from the
-   * results of this function.
-   *
-   * The return value is a tuple with the vector of vertices, the vector of
-   * cells, and a SubCellData structure. The latter contains additional
-   * information about faces and lines. These three objects are exactly
-   * the arguments to Triangulation::create_triangulation().
-   *
-   * This function is useful in cases where one needs to deconstruct a
-   * Triangulation or manipulate the numbering of the vertices in some way: an
-   * example is GridGenerator::merge_triangulations().
-   */
-  template <int dim, int spacedim>
-  std::
-    tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
-    get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria);
-
-  /** @} */
   /**
    * @name Functions supporting the creation of meshes
    */
@@ -536,6 +257,26 @@ namespace GridTools
    * @name Rotating, stretching and otherwise transforming meshes
    */
   /** @{ */
+  /**
+   * Return the arrays that define the coarse mesh of a Triangulation. This
+   * function is the inverse of Triangulation::create_triangulation() in the
+   * sense that if one called this function on a triangulation, then that
+   * triangulation could be recreated by some kind of refinement from the
+   * results of this function.
+   *
+   * The return value is a tuple with the vector of vertices, the vector of
+   * cells, and a SubCellData structure. The latter contains additional
+   * information about faces and lines. These three objects are exactly
+   * the arguments to Triangulation::create_triangulation().
+   *
+   * This function is useful in cases where one needs to deconstruct a
+   * Triangulation or manipulate the numbering of the vertices in some way: an
+   * example is GridGenerator::merge_triangulations().
+   */
+  template <int dim, int spacedim>
+  std::
+    tuple<std::vector<Point<spacedim>>, std::vector<CellData<dim>>, SubCellData>
+    get_coarse_mesh_description(const Triangulation<dim, spacedim> &tria);
 
   /**
    * Transform the vertices of the given triangulation by applying the
@@ -2130,34 +1871,6 @@ namespace GridTools
       const MeshType &mesh,
       const double    layer_thickness);
 
-  /**
-   * Compute and return a bounding box, defined through a pair of points
-   * bottom left and top right, that surrounds a subdomain of the @p mesh.
-   * Here, the "subdomain" consists of exactly all of those
-   * active cells for which the @p predicate returns @p true.
-   *
-   * For a description of how @p predicate works,
-   * see compute_active_cell_halo_layer().
-   *
-   * @note This function was written before the BoundingBox class was invented.
-   *   Consequently, it returns a pair of points, rather than a BoundingBox
-   * object as one may expect. However, BoundingBox has a conversion constructor
-   * from pairs of points, so the result of this function can still be assigned
-   * to a BoundingBox object.
-   *
-   * @dealiiConceptRequires{concepts::is_triangulation_or_dof_handler<MeshType>}
-   */
-  template <typename MeshType>
-  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
-  std::pair<
-    Point<MeshType::space_dimension>,
-    Point<MeshType::
-            space_dimension>> compute_bounding_box(const MeshType &mesh,
-                                                   const std::function<bool(
-                                                     const typename MeshType::
-                                                       active_cell_iterator &)>
-                                                     &predicate);
-
   /**
    * Compute a collection of bounding boxes so that all active cells for which
    * the given predicate is true, are completely enclosed in at least one of the
@@ -4054,549 +3767,6 @@ namespace GridTools
 
 
 
-  namespace internal
-  {
-    namespace ProjectToObject
-    {
-      /**
-       * The method GridTools::project_to_object requires taking derivatives
-       * along the surface of a simplex. In general these cannot be
-       * approximated with finite differences but special differences of the
-       * form
-       *
-       *     df/dx_i - df/dx_j
-       *
-       * <em>can</em> be approximated. This <code>struct</code> just stores
-       * the two derivatives approximated by the stencil (in the case of the
-       * example above <code>i</code> and <code>j</code>).
-       */
-      struct CrossDerivative
-      {
-        const unsigned int direction_0;
-        const unsigned int direction_1;
-
-        CrossDerivative(const unsigned int d0, const unsigned int d1);
-      };
-
-      inline CrossDerivative::CrossDerivative(const unsigned int d0,
-                                              const unsigned int d1)
-        : direction_0(d0)
-        , direction_1(d1)
-      {}
-
-
-
-      /**
-       * Standard second-order approximation to the first derivative with a
-       * two-point centered scheme. This is used below in a 1d Newton method.
-       */
-      template <typename F>
-      inline auto
-      centered_first_difference(const double center,
-                                const double step,
-                                const F &f) -> decltype(f(center) - f(center))
-      {
-        return (f(center + step) - f(center - step)) / (2.0 * step);
-      }
-
-
-
-      /**
-       * Standard second-order approximation to the second derivative with a
-       * three-point centered scheme. This is used below in a 1d Newton method.
-       */
-      template <typename F>
-      inline auto
-      centered_second_difference(const double center,
-                                 const double step,
-                                 const F &f) -> decltype(f(center) - f(center))
-      {
-        return (f(center + step) - 2.0 * f(center) + f(center - step)) /
-               (step * step);
-      }
-
-
-
-      /**
-       * Fourth order approximation of the derivative
-       *
-       *     df/dx_i - df/dx_j
-       *
-       * where <code>i</code> and <code>j</code> are specified by @p
-       * cross_derivative. The derivative approximation is at @p center with a
-       * step size of @p step and function @p f.
-       */
-      template <int structdim, typename F>
-      inline auto
-      cross_stencil(
-        const CrossDerivative cross_derivative,
-        const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
-        const double                                                 step,
-        const F &f) -> decltype(f(center) - f(center))
-      {
-        Tensor<1, GeometryInfo<structdim>::vertices_per_cell> simplex_vector;
-        simplex_vector[cross_derivative.direction_0] = 0.5 * step;
-        simplex_vector[cross_derivative.direction_1] = -0.5 * step;
-        return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
-                1.0 / 3.0 * f(center - simplex_vector) +
-                16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
-               step;
-      }
-
-
-
-      /**
-       * The optimization algorithm used in GridTools::project_to_object is
-       * essentially a gradient descent method. This function computes entries
-       * in the gradient of the objective function; see the description in the
-       * comments inside GridTools::project_to_object for more information.
-       */
-      template <int spacedim, int structdim, typename F>
-      inline double
-      gradient_entry(
-        const unsigned int     row_n,
-        const unsigned int     dependent_direction,
-        const Point<spacedim> &p0,
-        const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
-        const double                                                 step,
-        const F                                                     &f)
-      {
-        Assert(row_n < GeometryInfo<structdim>::vertices_per_cell &&
-                 dependent_direction <
-                   GeometryInfo<structdim>::vertices_per_cell,
-               ExcMessage("This function assumes that the last weight is a "
-                          "dependent variable (and hence we cannot take its "
-                          "derivative directly)."));
-        Assert(row_n != dependent_direction,
-               ExcMessage(
-                 "We cannot differentiate with respect to the variable "
-                 "that is assumed to be dependent."));
-
-        const Point<spacedim>     manifold_point = f(center);
-        const Tensor<1, spacedim> stencil_value  = cross_stencil<structdim>(
-          {row_n, dependent_direction}, center, step, f);
-        double entry = 0.0;
-        for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
-          entry +=
-            -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
-        return entry;
-      }
-
-      /**
-       * Project onto a d-linear object. This is more accurate than the
-       * general algorithm in project_to_object but only works for geometries
-       * described by linear, bilinear, or trilinear mappings.
-       */
-      template <typename Iterator, int spacedim, int structdim>
-      Point<spacedim>
-      project_to_d_linear_object(const Iterator        &object,
-                                 const Point<spacedim> &trial_point)
-      {
-        // let's look at this for simplicity for a quadrilateral
-        // (structdim==2) in a space with spacedim>2 (notate trial_point by
-        // y): all points on the surface are given by
-        //   x(\xi) = sum_i v_i phi_x(\xi)
-        // where v_i are the vertices of the quadrilateral, and
-        // \xi=(\xi_1,\xi_2) are the reference coordinates of the
-        // quadrilateral. so what we are trying to do is find a point x on the
-        // surface that is closest to the point y. there are different ways to
-        // solve this problem, but in the end it's a nonlinear problem and we
-        // have to find reference coordinates \xi so that J(\xi) = 1/2 ||
-        // x(\xi)-y ||^2 is minimal. x(\xi) is a function that is
-        // structdim-linear in \xi, so J(\xi) is a polynomial of degree
-        // 2*structdim that we'd like to minimize. unless structdim==1, we'll
-        // have to use a Newton method to find the answer. This leads to the
-        // following formulation of Newton steps:
-        //
-        // Given \xi_k, find \delta\xi_k so that
-        //   H_k \delta\xi_k = - F_k
-        // where H_k is an approximation to the second derivatives of J at
-        // \xi_k, and F_k is the first derivative of J.  We'll iterate this a
-        // number of times until the right hand side is small enough. As a
-        // stopping criterion, we terminate if ||\delta\xi||<eps.
-        //
-        // As for the Hessian, the best choice would be
-        //   H_k = J''(\xi_k)
-        // but we'll opt for the simpler Gauss-Newton form
-        //   H_k = A^T A
-        // i.e.
-        //   (H_k)_{nm} = \sum_{i,j} v_i*v_j *
-        //                   \partial_n phi_i *
-        //                   \partial_m phi_j
-        // we start at xi=(0.5, 0.5).
-        Point<structdim> xi;
-        for (unsigned int d = 0; d < structdim; ++d)
-          xi[d] = 0.5;
-
-        Point<spacedim> x_k;
-        for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
-          x_k += object->vertex(i) *
-                 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
-
-        do
-          {
-            Tensor<1, structdim> F_k;
-            for (const unsigned int i :
-                 GeometryInfo<structdim>::vertex_indices())
-              F_k +=
-                (x_k - trial_point) * object->vertex(i) *
-                GeometryInfo<structdim>::d_linear_shape_function_gradient(xi,
-                                                                          i);
-
-            Tensor<2, structdim> H_k;
-            for (const unsigned int i :
-                 GeometryInfo<structdim>::vertex_indices())
-              for (const unsigned int j :
-                   GeometryInfo<structdim>::vertex_indices())
-                {
-                  Tensor<2, structdim> tmp = outer_product(
-                    GeometryInfo<structdim>::d_linear_shape_function_gradient(
-                      xi, i),
-                    GeometryInfo<structdim>::d_linear_shape_function_gradient(
-                      xi, j));
-                  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
-                }
-
-            const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
-            xi += delta_xi;
-
-            x_k = Point<spacedim>();
-            for (const unsigned int i :
-                 GeometryInfo<structdim>::vertex_indices())
-              x_k += object->vertex(i) *
-                     GeometryInfo<structdim>::d_linear_shape_function(xi, i);
-
-            if (delta_xi.norm() < 1e-7)
-              break;
-          }
-        while (true);
-
-        return x_k;
-      }
-    } // namespace ProjectToObject
-  }   // namespace internal
-
-
-
-  namespace internal
-  {
-    // We hit an internal compiler error in ICC 15 if we define this as a lambda
-    // inside the project_to_object function below.
-    template <int structdim>
-    inline bool
-    weights_are_ok(
-      const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &v)
-    {
-      // clang has trouble figuring out structdim here, so define it
-      // again:
-      static const std::size_t n_vertices_per_cell =
-        Tensor<1, GeometryInfo<structdim>::vertices_per_cell>::
-          n_independent_components;
-      std::array<double, n_vertices_per_cell> copied_weights;
-      for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
-        {
-          copied_weights[i] = v[i];
-          if (v[i] < 0.0 || v[i] > 1.0)
-            return false;
-        }
-
-      // check the sum: try to avoid some roundoff errors by summing in order
-      std::sort(copied_weights.begin(), copied_weights.end());
-      const double sum =
-        std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
-      return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
-    }
-  } // namespace internal
-
-  template <typename Iterator>
-  Point<Iterator::AccessorType::space_dimension>
-  project_to_object(
-    const Iterator                                       &object,
-    const Point<Iterator::AccessorType::space_dimension> &trial_point)
-  {
-    const int spacedim  = Iterator::AccessorType::space_dimension;
-    const int structdim = Iterator::AccessorType::structure_dimension;
-
-    Point<spacedim> projected_point = trial_point;
-
-    if (structdim >= spacedim)
-      return projected_point;
-    else if (structdim == 1 || structdim == 2)
-      {
-        using namespace internal::ProjectToObject;
-        // Try to use the special flat algorithm for quads (this is better
-        // than the general algorithm in 3d). This does not take into account
-        // whether projected_point is outside the quad, but we optimize along
-        // lines below anyway:
-        const int                      dim = Iterator::AccessorType::dimension;
-        const Manifold<dim, spacedim> &manifold = object->get_manifold();
-        if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
-                                &manifold) != nullptr)
-          {
-            projected_point =
-              project_to_d_linear_object<Iterator, spacedim, structdim>(
-                object, trial_point);
-          }
-        else
-          {
-            // We want to find a point on the convex hull (defined by the
-            // vertices of the object and the manifold description) that is
-            // relatively close to the trial point. This has a few issues:
-            //
-            // 1. For a general convex hull we are not guaranteed that a unique
-            //    minimum exists.
-            // 2. The independent variables in the optimization process are the
-            //    weights given to Manifold::get_new_point, which must sum to 1,
-            //    so we cannot use standard finite differences to approximate a
-            //    gradient.
-            //
-            // There is not much we can do about 1., but for 2. we can derive
-            // finite difference stencils that work on a structdim-dimensional
-            // simplex and rewrite the optimization problem to use those
-            // instead. Consider the structdim 2 case and let
-            //
-            // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
-            // c2, c3})
-            //
-            // where {c0, c1, c2, c3} are the weights for the four vertices on
-            // the quadrilateral. We seek to minimize the Euclidean distance
-            // between F(...) and trial_point. We can solve for c3 in terms of
-            // the other weights and get, for one coordinate direction
-            //
-            // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
-            //      = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
-            //
-            // where we substitute back in for c3 after taking the
-            // derivative. We can compute a stencil for the cross derivative
-            // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
-            // (and gradient_entry computes the sum over the independent
-            // variables). Below, we somewhat arbitrarily pick the last
-            // component as the dependent one.
-            //
-            // Since we can now calculate derivatives of the objective
-            // function we can use gradient descent to minimize it.
-            //
-            // Of course, this is much simpler in the structdim = 1 case (we
-            // could rewrite the projection as a 1d optimization problem), but
-            // to reduce the potential for bugs we use the same code in both
-            // cases.
-            const double step_size = object->diameter() / 64.0;
-
-            constexpr unsigned int n_vertices_per_cell =
-              GeometryInfo<structdim>::vertices_per_cell;
-
-            std::array<Point<spacedim>, n_vertices_per_cell> vertices;
-            for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
-                 ++vertex_n)
-              vertices[vertex_n] = object->vertex(vertex_n);
-
-            auto get_point_from_weights =
-              [&](const Tensor<1, n_vertices_per_cell> &weights)
-              -> Point<spacedim> {
-              return object->get_manifold().get_new_point(
-                make_array_view(vertices.begin(), vertices.end()),
-                make_array_view(weights.begin_raw(), weights.end_raw()));
-            };
-
-            // pick the initial weights as (normalized) inverse distances from
-            // the trial point:
-            Tensor<1, n_vertices_per_cell> guess_weights;
-            double                         guess_weights_sum = 0.0;
-            for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
-                 ++vertex_n)
-              {
-                const double distance =
-                  vertices[vertex_n].distance(trial_point);
-                if (distance == 0.0)
-                  {
-                    guess_weights           = 0.0;
-                    guess_weights[vertex_n] = 1.0;
-                    guess_weights_sum       = 1.0;
-                    break;
-                  }
-                else
-                  {
-                    guess_weights[vertex_n] = 1.0 / distance;
-                    guess_weights_sum += guess_weights[vertex_n];
-                  }
-              }
-            guess_weights /= guess_weights_sum;
-            Assert(internal::weights_are_ok<structdim>(guess_weights),
-                   ExcInternalError());
-
-            // The optimization algorithm consists of two parts:
-            //
-            // 1. An outer loop where we apply the gradient descent algorithm.
-            // 2. An inner loop where we do a line search to find the optimal
-            //    length of the step one should take in the gradient direction.
-            //
-            for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
-              {
-                const unsigned int dependent_direction =
-                  n_vertices_per_cell - 1;
-                Tensor<1, n_vertices_per_cell> current_gradient;
-                for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
-                     ++row_n)
-                  {
-                    if (row_n != dependent_direction)
-                      {
-                        current_gradient[row_n] =
-                          gradient_entry<spacedim, structdim>(
-                            row_n,
-                            dependent_direction,
-                            trial_point,
-                            guess_weights,
-                            step_size,
-                            get_point_from_weights);
-
-                        current_gradient[dependent_direction] -=
-                          current_gradient[row_n];
-                      }
-                  }
-
-                // We need to travel in the -gradient direction, as noted
-                // above, but we may not want to take a full step in that
-                // direction; instead, guess that we will go -0.5*gradient and
-                // do quasi-Newton iteration to pick the best multiplier. The
-                // goal is to find a scalar alpha such that
-                //
-                // F(x - alpha g)
-                //
-                // is minimized, where g is the gradient and F is the
-                // objective function. To find the optimal value we find roots
-                // of the derivative of the objective function with respect to
-                // alpha by Newton iteration, where we approximate the first
-                // and second derivatives of F(x - alpha g) with centered
-                // finite differences.
-                double gradient_weight = -0.5;
-                auto   gradient_weight_objective_function =
-                  [&](const double gradient_weight_guess) -> double {
-                  return (trial_point -
-                          get_point_from_weights(guess_weights +
-                                                 gradient_weight_guess *
-                                                   current_gradient))
-                    .norm_square();
-                };
-
-                for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
-                  {
-                    const double update_numerator = centered_first_difference(
-                      gradient_weight,
-                      step_size,
-                      gradient_weight_objective_function);
-                    const double update_denominator =
-                      centered_second_difference(
-                        gradient_weight,
-                        step_size,
-                        gradient_weight_objective_function);
-
-                    // avoid division by zero. Note that we limit the gradient
-                    // weight below
-                    if (std::abs(update_denominator) == 0.0)
-                      break;
-                    gradient_weight =
-                      gradient_weight - update_numerator / update_denominator;
-
-                    // Put a fairly lenient bound on the largest possible
-                    // gradient (things tend to be locally flat, so the gradient
-                    // itself is usually small)
-                    if (std::abs(gradient_weight) > 10)
-                      {
-                        gradient_weight = -10.0;
-                        break;
-                      }
-                  }
-
-                // It only makes sense to take convex combinations with weights
-                // between zero and one. If the update takes us outside of this
-                // region then rescale the update to stay within the region and
-                // try again
-                Tensor<1, n_vertices_per_cell> tentative_weights =
-                  guess_weights + gradient_weight * current_gradient;
-
-                double new_gradient_weight = gradient_weight;
-                for (unsigned int iteration_count = 0; iteration_count < 40;
-                     ++iteration_count)
-                  {
-                    if (internal::weights_are_ok<structdim>(tentative_weights))
-                      break;
-
-                    for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
-                      {
-                        if (tentative_weights[i] < 0.0)
-                          {
-                            tentative_weights -=
-                              (tentative_weights[i] / current_gradient[i]) *
-                              current_gradient;
-                          }
-                        if (tentative_weights[i] < 0.0 ||
-                            1.0 < tentative_weights[i])
-                          {
-                            new_gradient_weight /= 2.0;
-                            tentative_weights =
-                              guess_weights +
-                              new_gradient_weight * current_gradient;
-                          }
-                      }
-                  }
-
-                // the update might still send us outside the valid region, so
-                // check again and quit if the update is still not valid
-                if (!internal::weights_are_ok<structdim>(tentative_weights))
-                  break;
-
-                // if we cannot get closer by traveling in the gradient
-                // direction then quit
-                if (get_point_from_weights(tentative_weights)
-                      .distance(trial_point) <
-                    get_point_from_weights(guess_weights).distance(trial_point))
-                  guess_weights = tentative_weights;
-                else
-                  break;
-                Assert(internal::weights_are_ok<structdim>(guess_weights),
-                       ExcInternalError());
-              }
-            Assert(internal::weights_are_ok<structdim>(guess_weights),
-                   ExcInternalError());
-            projected_point = get_point_from_weights(guess_weights);
-          }
-
-        // if structdim == 2 and the optimal point is not on the interior then
-        // we may be able to get a more accurate result by projecting onto the
-        // lines.
-        if (structdim == 2)
-          {
-            std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
-              line_projections;
-            for (unsigned int line_n = 0;
-                 line_n < GeometryInfo<structdim>::lines_per_cell;
-                 ++line_n)
-              {
-                line_projections[line_n] =
-                  project_to_object(object->line(line_n), trial_point);
-              }
-            std::sort(line_projections.begin(),
-                      line_projections.end(),
-                      [&](const Point<spacedim> &a, const Point<spacedim> &b) {
-                        return a.distance(trial_point) <
-                               b.distance(trial_point);
-                      });
-            if (line_projections[0].distance(trial_point) <
-                projected_point.distance(trial_point))
-              projected_point = line_projections[0];
-          }
-      }
-    else
-      {
-        Assert(false, ExcNotImplemented());
-        return projected_point;
-      }
-
-    return projected_point;
-  }
-
-
-
   namespace internal
   {
     template <typename DataType,
diff --git a/include/deal.II/grid/grid_tools_geometry.h b/include/deal.II/grid/grid_tools_geometry.h
new file mode 100644 (file)
index 0000000..efe0cc8
--- /dev/null
@@ -0,0 +1,879 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_grid_tools_geometry_h
+#define dealii_grid_tools_geometry_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/fe/mapping.h>
+
+#include <deal.II/grid/manifold.h>
+#include <deal.II/grid/tria.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * This namespace is a collection of algorithms working on triangulations,
+ * such as shifting or rotating triangulations, but also finding a cell that
+ * contains a given point. See the descriptions of the individual functions
+ * for more information.
+ *
+ * @ingroup grid
+ */
+namespace GridTools
+{
+  /**
+   * @name Information about meshes and cells
+   */
+  /** @{ */
+
+  /**
+   * Return the diameter of a triangulation. The diameter is computed using
+   * only the vertices, i.e. if the diameter should be larger than the maximal
+   * distance between boundary vertices due to a higher order mapping, then
+   * this function will not catch this.
+   */
+  template <int dim, int spacedim>
+  double
+  diameter(const Triangulation<dim, spacedim> &tria);
+
+  /**
+   * Compute the volume (i.e. the dim-dimensional measure) of the
+   * triangulation. We compute the measure using the integral $\sum_K \int_K 1
+   * \; dx$ where $K$ are the cells of the given triangulation. The integral
+   * is approximated via quadrature. This version of the function uses a
+   * linear mapping to compute the JxW values on each cell.
+   *
+   * If the triangulation is a dim-dimensional one embedded in a higher
+   * dimensional space of dimension spacedim, then the value returned is the
+   * dim-dimensional measure. For example, for a two-dimensional triangulation
+   * in three-dimensional space, the value returned is the area of the surface
+   * so described. (This obviously makes sense since the spacedim-dimensional
+   * measure of a dim-dimensional triangulation would always be zero if dim @<
+   * spacedim).
+   *
+   * This function also works for objects of type
+   * parallel::distributed::Triangulation, in which case the function is a
+   * collective operation.
+   *
+   * @param tria The triangulation.
+   * @return The dim-dimensional measure of the domain described by the
+   * triangulation, as discussed above.
+   */
+  template <int dim, int spacedim>
+  double
+  volume(const Triangulation<dim, spacedim> &tria);
+
+  /**
+   * Compute the volume (i.e. the dim-dimensional measure) of the
+   * triangulation. We compute the measure using the integral $\sum_K \int_K 1
+   * \; dx$ where $K$ are the cells of the given triangulation. The integral
+   * is approximated via quadrature for which we use the mapping argument.
+   *
+   * If the triangulation is a dim-dimensional one embedded in a higher
+   * dimensional space of dimension spacedim, then the value returned is the
+   * dim-dimensional measure. For example, for a two-dimensional triangulation
+   * in three-dimensional space, the value returned is the area of the surface
+   * so described. (This obviously makes sense since the spacedim-dimensional
+   * measure of a dim-dimensional triangulation would always be zero if dim @<
+   * spacedim.
+   *
+   * This function also works for objects of type
+   * parallel::distributed::Triangulation, in which case the function is a
+   * collective operation.
+   *
+   * @param tria The triangulation.
+   * @param mapping The Mapping which computes the Jacobians used to
+   * approximate the volume via quadrature. Explicitly using a higher-order
+   * Mapping (i.e., instead of using the other version of this function) will
+   * result in a more accurate approximation of the volume on Triangulations
+   * with curvature described by Manifold objects.
+   * @return The dim-dimensional measure of the domain described by the
+   * triangulation, as discussed above.
+   */
+  template <int dim, int spacedim>
+  double
+  volume(const Triangulation<dim, spacedim> &tria,
+         const Mapping<dim, spacedim>       &mapping);
+
+  /**
+   * Return an approximation of the diameter of the smallest active cell of a
+   * triangulation. See step-24 for an example of use of this function.
+   *
+   * Notice that, even if you pass a non-trivial mapping, the returned value is
+   * computed only using information on the vertices of the triangulation,
+   * possibly transformed by the mapping. While this is accurate most of the
+   * times, it may fail to give the correct result when the triangulation
+   * contains very distorted cells.
+   */
+  template <int dim, int spacedim>
+  double
+  minimal_cell_diameter(
+    const Triangulation<dim, spacedim> &triangulation,
+    const Mapping<dim, spacedim>       &mapping =
+      (ReferenceCells::get_hypercube<dim>()
+#ifndef _MSC_VER
+         .template get_default_linear_mapping<dim, spacedim>()
+#else
+         .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
+#endif
+         ));
+
+  /**
+   * Return an approximation of the diameter of the largest active cell of a
+   * triangulation.
+   *
+   * Notice that, even if you pass a non-trivial mapping to this function, the
+   * returned value is computed only using information on the vertices of the
+   * triangulation, possibly transformed by the mapping. While this is accurate
+   * most of the times, it may fail to give the correct result when the
+   * triangulation contains very distorted cells.
+   */
+  template <int dim, int spacedim>
+  double
+  maximal_cell_diameter(
+    const Triangulation<dim, spacedim> &triangulation,
+    const Mapping<dim, spacedim>       &mapping =
+      (ReferenceCells::get_hypercube<dim>()
+#ifndef _MSC_VER
+         .template get_default_linear_mapping<dim, spacedim>()
+#else
+         .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
+#endif
+         ));
+
+  /**
+   * Given a list of vertices (typically obtained using
+   * Triangulation::get_vertices()) as the first, and a list of vertex indices
+   * that characterize a single cell as the second argument, return the
+   * measure (area, volume) of this cell. If this is a real cell, then you can
+   * get the same result using <code>cell-@>measure()</code>, but this
+   * function also works for cells that do not exist except that you make it
+   * up by naming its vertices from the list.
+   *
+   * The size of @p vertex_indices, combined with `dim`, implicitly encodes
+   * the ReferenceCell type of the provided cell. For example, if `dim == 2` and
+   * `vertex_indices.size() == 3` then the cell is a triangle, but if
+   * `dim == 2` and `vertex_indices.size() == 4` then the cell is a
+   * quadrilateral. A std::vector is implicitly convertible to an ArrayView, so
+   * it can be passed directly to this function. See the ArrayView class for
+   * more information.
+   *
+   * @note This function is only implemented for codimension zero objects.
+   */
+  template <int dim>
+  double
+  cell_measure(const std::vector<Point<dim>>       &all_vertices,
+               const ArrayView<const unsigned int> &vertex_indices);
+
+  /**
+   * This function computes an affine approximation of the map from the unit
+   * coordinates to the real coordinates of the form $p_\text{real} = A
+   * p_\text{unit} + b $ by a least squares fit of this affine function to the
+   * $2^\text{dim}$ vertices representing a quadrilateral or hexahedral cell
+   * in `spacedim` dimensions. The result is returned as a pair with the
+   * matrix <i>A</i> as the first argument and the vector <i>b</i> describing
+   * distance of the plane to the origin.
+   *
+   * For any valid mesh cell whose geometry is not degenerate, this operation
+   * results in a unique affine mapping, even in cases where the actual
+   * transformation by a bi-/trilinear or higher order mapping might be
+   * singular. The result is exact in case the transformation from the unit to
+   * the real cell is indeed affine, such as in one dimension or for Cartesian
+   * and affine (parallelogram) meshes in 2d/3d.
+   *
+   * This approximation is underlying the function
+   * TriaAccessor::real_to_unit_cell_affine_approximation() function.
+   *
+   * For exact transformations to the unit cell, use
+   * Mapping::transform_real_to_unit_cell().
+   */
+  template <int dim, int spacedim>
+  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
+  affine_cell_approximation(const ArrayView<const Point<spacedim>> &vertices);
+
+  /**
+   * Computes an aspect ratio measure for all locally-owned active cells and
+   * fills a vector with one entry per cell, given a @p triangulation and
+   * @p mapping. The size of the vector that is returned equals the number of
+   * active cells. The vector contains zero for non locally-owned cells. The
+   * aspect ratio of a cell is defined as the ratio of the maximum to minimum
+   * singular value of the Jacobian, taking the maximum over all quadrature
+   * points of a quadrature rule specified via @p quadrature. For example, for
+   * the special case of rectangular elements in 2d with dimensions $a$ and $b$
+   * ($a \geq b$), this function returns the usual aspect ratio definition
+   * $a/b$. The above definition using singular values is a generalization to
+   * arbitrarily deformed elements. This function is intended to be used for
+   * $d=2,3$ space dimensions, but it can also be used for $d=1$ returning a
+   * value of 1.
+   *
+   * @note Inverted elements do not throw an exception. Instead, a value of inf
+   * is written into the vector in case of inverted elements.
+   *
+   * @note Make sure to use enough quadrature points for a precise calculation
+   * of the aspect ratio in case of deformed elements.
+   *
+   * @note In parallel computations the return value will have the length
+   * n_active_cells but the aspect ratio is only computed for the cells that
+   * are locally owned and placed at index CellAccessor::active_cell_index(),
+   * respectively. All other values are set to 0.
+   *
+   * @note This function can only be used if deal.II was configured with
+   * support for LAPACK.
+   */
+  template <int dim>
+  Vector<double>
+  compute_aspect_ratio_of_cells(const Mapping<dim>       &mapping,
+                                const Triangulation<dim> &triangulation,
+                                const Quadrature<dim>    &quadrature);
+
+  /**
+   * Computes the maximum aspect ratio by taking the maximum over all cells.
+   *
+   * @note When running in parallel with a Triangulation that supports MPI,
+   * this is a collective call and the return value is the maximum over all
+   * processors.
+   */
+  template <int dim>
+  double
+  compute_maximum_aspect_ratio(const Mapping<dim>       &mapping,
+                               const Triangulation<dim> &triangulation,
+                               const Quadrature<dim>    &quadrature);
+
+  /**
+   * Compute the smallest box containing the entire triangulation.
+   *
+   * If the input triangulation is a `parallel::distributed::Triangulation`,
+   * then each processor will compute a bounding box enclosing all locally
+   * owned, ghost, and artificial cells. In the case of a domain without curved
+   * boundaries, these bounding boxes will all agree between processors because
+   * the union of the areas occupied by artificial and ghost cells equals the
+   * union of the areas occupied by the cells that other processors own.
+   * However, if the domain has curved boundaries, this is no longer the case.
+   * The bounding box returned may be appropriate for the current processor,
+   * but different from the bounding boxes computed on other processors.
+   */
+  template <int dim, int spacedim>
+  BoundingBox<spacedim>
+  compute_bounding_box(const Triangulation<dim, spacedim> &triangulation);
+
+  /**
+   * Compute and return a bounding box, defined through a pair of points
+   * bottom left and top right, that surrounds a subdomain of the @p mesh.
+   * Here, the "subdomain" consists of exactly all of those
+   * active cells for which the @p predicate returns @p true.
+   *
+   * For a description of how @p predicate works,
+   * see compute_active_cell_halo_layer().
+   *
+   * @note This function was written before the BoundingBox class was invented.
+   *   Consequently, it returns a pair of points, rather than a BoundingBox
+   * object as one may expect. However, BoundingBox has a conversion constructor
+   * from pairs of points, so the result of this function can still be assigned
+   * to a BoundingBox object.
+   *
+   * @dealiiConceptRequires{concepts::is_triangulation_or_dof_handler<MeshType>}
+   */
+  template <typename MeshType>
+  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
+  std::pair<
+    Point<MeshType::space_dimension>,
+    Point<MeshType::
+            space_dimension>> compute_bounding_box(const MeshType &mesh,
+                                                   const std::function<bool(
+                                                     const typename MeshType::
+                                                       active_cell_iterator &)>
+                                                     &predicate);
+
+  /**
+   * Return the point on the geometrical object @p object closest to the given
+   * point @p trial_point. For example, if @p object is a one-dimensional line
+   * or edge, then the returned point will be a point on the geodesic that
+   * connects the vertices as the manifold associated with the object sees it
+   * (i.e., the geometric line may be curved if it lives in a higher
+   * dimensional space). If the iterator points to a quadrilateral in a higher
+   * dimensional space, then the returned point lies within the convex hull of
+   * the vertices of the quad as seen by the associated manifold.
+   *
+   * @note This projection is usually not well-posed since there may be
+   * multiple points on the object that minimize the distance. The algorithm
+   * used in this function is robust (and the output is guaranteed to be on
+   * the given @p object) but may only provide a few correct digits if the
+   * object has high curvature. If your manifold supports it then the
+   * specialized function Manifold::project_to_manifold() may perform better.
+   */
+  template <typename Iterator>
+  Point<Iterator::AccessorType::space_dimension>
+  project_to_object(
+    const Iterator                                       &object,
+    const Point<Iterator::AccessorType::space_dimension> &trial_point);
+  /** @} */
+} // namespace GridTools
+
+#ifndef DOXYGEN
+namespace GridTools
+{
+  namespace internal
+  {
+    namespace ProjectToObject
+    {
+      /**
+       * The method GridTools::project_to_object requires taking derivatives
+       * along the surface of a simplex. In general these cannot be
+       * approximated with finite differences but special differences of the
+       * form
+       *
+       *     df/dx_i - df/dx_j
+       *
+       * <em>can</em> be approximated. This <code>struct</code> just stores
+       * the two derivatives approximated by the stencil (in the case of the
+       * example above <code>i</code> and <code>j</code>).
+       */
+      struct CrossDerivative
+      {
+        const unsigned int direction_0;
+        const unsigned int direction_1;
+
+        CrossDerivative(const unsigned int d0, const unsigned int d1);
+      };
+
+      inline CrossDerivative::CrossDerivative(const unsigned int d0,
+                                              const unsigned int d1)
+        : direction_0(d0)
+        , direction_1(d1)
+      {}
+
+
+
+      /**
+       * Standard second-order approximation to the first derivative with a
+       * two-point centered scheme. This is used below in a 1d Newton method.
+       */
+      template <typename F>
+      inline auto
+      centered_first_difference(const double center,
+                                const double step,
+                                const F &f) -> decltype(f(center) - f(center))
+      {
+        return (f(center + step) - f(center - step)) / (2.0 * step);
+      }
+
+
+
+      /**
+       * Standard second-order approximation to the second derivative with a
+       * three-point centered scheme. This is used below in a 1d Newton method.
+       */
+      template <typename F>
+      inline auto
+      centered_second_difference(const double center,
+                                 const double step,
+                                 const F &f) -> decltype(f(center) - f(center))
+      {
+        return (f(center + step) - 2.0 * f(center) + f(center - step)) /
+               (step * step);
+      }
+
+
+
+      /**
+       * Fourth order approximation of the derivative
+       *
+       *     df/dx_i - df/dx_j
+       *
+       * where <code>i</code> and <code>j</code> are specified by @p
+       * cross_derivative. The derivative approximation is at @p center with a
+       * step size of @p step and function @p f.
+       */
+      template <int structdim, typename F>
+      inline auto
+      cross_stencil(
+        const CrossDerivative cross_derivative,
+        const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
+        const double                                                 step,
+        const F &f) -> decltype(f(center) - f(center))
+      {
+        Tensor<1, GeometryInfo<structdim>::vertices_per_cell> simplex_vector;
+        simplex_vector[cross_derivative.direction_0] = 0.5 * step;
+        simplex_vector[cross_derivative.direction_1] = -0.5 * step;
+        return (-4.0 * f(center) - 1.0 * f(center + simplex_vector) -
+                1.0 / 3.0 * f(center - simplex_vector) +
+                16.0 / 3.0 * f(center + 0.5 * simplex_vector)) /
+               step;
+      }
+
+
+
+      /**
+       * The optimization algorithm used in GridTools::project_to_object is
+       * essentially a gradient descent method. This function computes entries
+       * in the gradient of the objective function; see the description in the
+       * comments inside GridTools::project_to_object for more information.
+       */
+      template <int spacedim, int structdim, typename F>
+      inline double
+      gradient_entry(
+        const unsigned int     row_n,
+        const unsigned int     dependent_direction,
+        const Point<spacedim> &p0,
+        const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
+        const double                                                 step,
+        const F                                                     &f)
+      {
+        Assert(row_n < GeometryInfo<structdim>::vertices_per_cell &&
+                 dependent_direction <
+                   GeometryInfo<structdim>::vertices_per_cell,
+               ExcMessage("This function assumes that the last weight is a "
+                          "dependent variable (and hence we cannot take its "
+                          "derivative directly)."));
+        Assert(row_n != dependent_direction,
+               ExcMessage(
+                 "We cannot differentiate with respect to the variable "
+                 "that is assumed to be dependent."));
+
+        const Point<spacedim>     manifold_point = f(center);
+        const Tensor<1, spacedim> stencil_value  = cross_stencil<structdim>(
+          {row_n, dependent_direction}, center, step, f);
+        double entry = 0.0;
+        for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
+          entry +=
+            -2.0 * (p0[dim_n] - manifold_point[dim_n]) * stencil_value[dim_n];
+        return entry;
+      }
+
+      /**
+       * Project onto a d-linear object. This is more accurate than the
+       * general algorithm in project_to_object but only works for geometries
+       * described by linear, bilinear, or trilinear mappings.
+       */
+      template <typename Iterator, int spacedim, int structdim>
+      Point<spacedim>
+      project_to_d_linear_object(const Iterator        &object,
+                                 const Point<spacedim> &trial_point)
+      {
+        // let's look at this for simplicity for a quadrilateral
+        // (structdim==2) in a space with spacedim>2 (notate trial_point by
+        // y): all points on the surface are given by
+        //   x(\xi) = sum_i v_i phi_x(\xi)
+        // where v_i are the vertices of the quadrilateral, and
+        // \xi=(\xi_1,\xi_2) are the reference coordinates of the
+        // quadrilateral. so what we are trying to do is find a point x on the
+        // surface that is closest to the point y. there are different ways to
+        // solve this problem, but in the end it's a nonlinear problem and we
+        // have to find reference coordinates \xi so that J(\xi) = 1/2 ||
+        // x(\xi)-y ||^2 is minimal. x(\xi) is a function that is
+        // structdim-linear in \xi, so J(\xi) is a polynomial of degree
+        // 2*structdim that we'd like to minimize. unless structdim==1, we'll
+        // have to use a Newton method to find the answer. This leads to the
+        // following formulation of Newton steps:
+        //
+        // Given \xi_k, find \delta\xi_k so that
+        //   H_k \delta\xi_k = - F_k
+        // where H_k is an approximation to the second derivatives of J at
+        // \xi_k, and F_k is the first derivative of J.  We'll iterate this a
+        // number of times until the right hand side is small enough. As a
+        // stopping criterion, we terminate if ||\delta\xi||<eps.
+        //
+        // As for the Hessian, the best choice would be
+        //   H_k = J''(\xi_k)
+        // but we'll opt for the simpler Gauss-Newton form
+        //   H_k = A^T A
+        // i.e.
+        //   (H_k)_{nm} = \sum_{i,j} v_i*v_j *
+        //                   \partial_n phi_i *
+        //                   \partial_m phi_j
+        // we start at xi=(0.5, 0.5).
+        Point<structdim> xi;
+        for (unsigned int d = 0; d < structdim; ++d)
+          xi[d] = 0.5;
+
+        Point<spacedim> x_k;
+        for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
+          x_k += object->vertex(i) *
+                 GeometryInfo<structdim>::d_linear_shape_function(xi, i);
+
+        do
+          {
+            Tensor<1, structdim> F_k;
+            for (const unsigned int i :
+                 GeometryInfo<structdim>::vertex_indices())
+              F_k +=
+                (x_k - trial_point) * object->vertex(i) *
+                GeometryInfo<structdim>::d_linear_shape_function_gradient(xi,
+                                                                          i);
+
+            Tensor<2, structdim> H_k;
+            for (const unsigned int i :
+                 GeometryInfo<structdim>::vertex_indices())
+              for (const unsigned int j :
+                   GeometryInfo<structdim>::vertex_indices())
+                {
+                  Tensor<2, structdim> tmp = outer_product(
+                    GeometryInfo<structdim>::d_linear_shape_function_gradient(
+                      xi, i),
+                    GeometryInfo<structdim>::d_linear_shape_function_gradient(
+                      xi, j));
+                  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
+                }
+
+            const Tensor<1, structdim> delta_xi = -invert(H_k) * F_k;
+            xi += delta_xi;
+
+            x_k = Point<spacedim>();
+            for (const unsigned int i :
+                 GeometryInfo<structdim>::vertex_indices())
+              x_k += object->vertex(i) *
+                     GeometryInfo<structdim>::d_linear_shape_function(xi, i);
+
+            if (delta_xi.norm() < 1e-7)
+              break;
+          }
+        while (true);
+
+        return x_k;
+      }
+    } // namespace ProjectToObject
+
+    // We hit an internal compiler error in ICC 15 if we define this as a lambda
+    // inside the project_to_object function below.
+    template <int structdim>
+    inline bool
+    weights_are_ok(
+      const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &v)
+    {
+      // clang has trouble figuring out structdim here, so define it
+      // again:
+      static const std::size_t n_vertices_per_cell =
+        Tensor<1, GeometryInfo<structdim>::vertices_per_cell>::
+          n_independent_components;
+      std::array<double, n_vertices_per_cell> copied_weights;
+      for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
+        {
+          copied_weights[i] = v[i];
+          if (v[i] < 0.0 || v[i] > 1.0)
+            return false;
+        }
+
+      // check the sum: try to avoid some roundoff errors by summing in order
+      std::sort(copied_weights.begin(), copied_weights.end());
+      const double sum =
+        std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
+      return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
+    }
+  } // namespace internal
+
+  template <typename Iterator>
+  Point<Iterator::AccessorType::space_dimension>
+  project_to_object(
+    const Iterator                                       &object,
+    const Point<Iterator::AccessorType::space_dimension> &trial_point)
+  {
+    const int spacedim  = Iterator::AccessorType::space_dimension;
+    const int structdim = Iterator::AccessorType::structure_dimension;
+
+    Point<spacedim> projected_point = trial_point;
+
+    if (structdim >= spacedim)
+      return projected_point;
+    else if (structdim == 1 || structdim == 2)
+      {
+        using namespace internal::ProjectToObject;
+        // Try to use the special flat algorithm for quads (this is better
+        // than the general algorithm in 3d). This does not take into account
+        // whether projected_point is outside the quad, but we optimize along
+        // lines below anyway:
+        const int                      dim = Iterator::AccessorType::dimension;
+        const Manifold<dim, spacedim> &manifold = object->get_manifold();
+        if (structdim == 2 && dynamic_cast<const FlatManifold<dim, spacedim> *>(
+                                &manifold) != nullptr)
+          {
+            projected_point =
+              project_to_d_linear_object<Iterator, spacedim, structdim>(
+                object, trial_point);
+          }
+        else
+          {
+            // We want to find a point on the convex hull (defined by the
+            // vertices of the object and the manifold description) that is
+            // relatively close to the trial point. This has a few issues:
+            //
+            // 1. For a general convex hull we are not guaranteed that a unique
+            //    minimum exists.
+            // 2. The independent variables in the optimization process are the
+            //    weights given to Manifold::get_new_point, which must sum to 1,
+            //    so we cannot use standard finite differences to approximate a
+            //    gradient.
+            //
+            // There is not much we can do about 1., but for 2. we can derive
+            // finite difference stencils that work on a structdim-dimensional
+            // simplex and rewrite the optimization problem to use those
+            // instead. Consider the structdim 2 case and let
+            //
+            // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1,
+            // c2, c3})
+            //
+            // where {c0, c1, c2, c3} are the weights for the four vertices on
+            // the quadrilateral. We seek to minimize the Euclidean distance
+            // between F(...) and trial_point. We can solve for c3 in terms of
+            // the other weights and get, for one coordinate direction
+            //
+            // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
+            //      = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
+            //
+            // where we substitute back in for c3 after taking the
+            // derivative. We can compute a stencil for the cross derivative
+            // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
+            // (and gradient_entry computes the sum over the independent
+            // variables). Below, we somewhat arbitrarily pick the last
+            // component as the dependent one.
+            //
+            // Since we can now calculate derivatives of the objective
+            // function we can use gradient descent to minimize it.
+            //
+            // Of course, this is much simpler in the structdim = 1 case (we
+            // could rewrite the projection as a 1d optimization problem), but
+            // to reduce the potential for bugs we use the same code in both
+            // cases.
+            const double step_size = object->diameter() / 64.0;
+
+            constexpr unsigned int n_vertices_per_cell =
+              GeometryInfo<structdim>::vertices_per_cell;
+
+            std::array<Point<spacedim>, n_vertices_per_cell> vertices;
+            for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
+                 ++vertex_n)
+              vertices[vertex_n] = object->vertex(vertex_n);
+
+            auto get_point_from_weights =
+              [&](const Tensor<1, n_vertices_per_cell> &weights)
+              -> Point<spacedim> {
+              return object->get_manifold().get_new_point(
+                make_array_view(vertices.begin(), vertices.end()),
+                make_array_view(weights.begin_raw(), weights.end_raw()));
+            };
+
+            // pick the initial weights as (normalized) inverse distances from
+            // the trial point:
+            Tensor<1, n_vertices_per_cell> guess_weights;
+            double                         guess_weights_sum = 0.0;
+            for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
+                 ++vertex_n)
+              {
+                const double distance =
+                  vertices[vertex_n].distance(trial_point);
+                if (distance == 0.0)
+                  {
+                    guess_weights           = 0.0;
+                    guess_weights[vertex_n] = 1.0;
+                    guess_weights_sum       = 1.0;
+                    break;
+                  }
+                else
+                  {
+                    guess_weights[vertex_n] = 1.0 / distance;
+                    guess_weights_sum += guess_weights[vertex_n];
+                  }
+              }
+            guess_weights /= guess_weights_sum;
+            Assert(internal::weights_are_ok<structdim>(guess_weights),
+                   ExcInternalError());
+
+            // The optimization algorithm consists of two parts:
+            //
+            // 1. An outer loop where we apply the gradient descent algorithm.
+            // 2. An inner loop where we do a line search to find the optimal
+            //    length of the step one should take in the gradient direction.
+            //
+            for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
+              {
+                const unsigned int dependent_direction =
+                  n_vertices_per_cell - 1;
+                Tensor<1, n_vertices_per_cell> current_gradient;
+                for (unsigned int row_n = 0; row_n < n_vertices_per_cell;
+                     ++row_n)
+                  {
+                    if (row_n != dependent_direction)
+                      {
+                        current_gradient[row_n] =
+                          gradient_entry<spacedim, structdim>(
+                            row_n,
+                            dependent_direction,
+                            trial_point,
+                            guess_weights,
+                            step_size,
+                            get_point_from_weights);
+
+                        current_gradient[dependent_direction] -=
+                          current_gradient[row_n];
+                      }
+                  }
+
+                // We need to travel in the -gradient direction, as noted
+                // above, but we may not want to take a full step in that
+                // direction; instead, guess that we will go -0.5*gradient and
+                // do quasi-Newton iteration to pick the best multiplier. The
+                // goal is to find a scalar alpha such that
+                //
+                // F(x - alpha g)
+                //
+                // is minimized, where g is the gradient and F is the
+                // objective function. To find the optimal value we find roots
+                // of the derivative of the objective function with respect to
+                // alpha by Newton iteration, where we approximate the first
+                // and second derivatives of F(x - alpha g) with centered
+                // finite differences.
+                double gradient_weight = -0.5;
+                auto   gradient_weight_objective_function =
+                  [&](const double gradient_weight_guess) -> double {
+                  return (trial_point -
+                          get_point_from_weights(guess_weights +
+                                                 gradient_weight_guess *
+                                                   current_gradient))
+                    .norm_square();
+                };
+
+                for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
+                  {
+                    const double update_numerator = centered_first_difference(
+                      gradient_weight,
+                      step_size,
+                      gradient_weight_objective_function);
+                    const double update_denominator =
+                      centered_second_difference(
+                        gradient_weight,
+                        step_size,
+                        gradient_weight_objective_function);
+
+                    // avoid division by zero. Note that we limit the gradient
+                    // weight below
+                    if (std::abs(update_denominator) == 0.0)
+                      break;
+                    gradient_weight =
+                      gradient_weight - update_numerator / update_denominator;
+
+                    // Put a fairly lenient bound on the largest possible
+                    // gradient (things tend to be locally flat, so the gradient
+                    // itself is usually small)
+                    if (std::abs(gradient_weight) > 10)
+                      {
+                        gradient_weight = -10.0;
+                        break;
+                      }
+                  }
+
+                // It only makes sense to take convex combinations with weights
+                // between zero and one. If the update takes us outside of this
+                // region then rescale the update to stay within the region and
+                // try again
+                Tensor<1, n_vertices_per_cell> tentative_weights =
+                  guess_weights + gradient_weight * current_gradient;
+
+                double new_gradient_weight = gradient_weight;
+                for (unsigned int iteration_count = 0; iteration_count < 40;
+                     ++iteration_count)
+                  {
+                    if (internal::weights_are_ok<structdim>(tentative_weights))
+                      break;
+
+                    for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
+                      {
+                        if (tentative_weights[i] < 0.0)
+                          {
+                            tentative_weights -=
+                              (tentative_weights[i] / current_gradient[i]) *
+                              current_gradient;
+                          }
+                        if (tentative_weights[i] < 0.0 ||
+                            1.0 < tentative_weights[i])
+                          {
+                            new_gradient_weight /= 2.0;
+                            tentative_weights =
+                              guess_weights +
+                              new_gradient_weight * current_gradient;
+                          }
+                      }
+                  }
+
+                // the update might still send us outside the valid region, so
+                // check again and quit if the update is still not valid
+                if (!internal::weights_are_ok<structdim>(tentative_weights))
+                  break;
+
+                // if we cannot get closer by traveling in the gradient
+                // direction then quit
+                if (get_point_from_weights(tentative_weights)
+                      .distance(trial_point) <
+                    get_point_from_weights(guess_weights).distance(trial_point))
+                  guess_weights = tentative_weights;
+                else
+                  break;
+                Assert(internal::weights_are_ok<structdim>(guess_weights),
+                       ExcInternalError());
+              }
+            Assert(internal::weights_are_ok<structdim>(guess_weights),
+                   ExcInternalError());
+            projected_point = get_point_from_weights(guess_weights);
+          }
+
+        // if structdim == 2 and the optimal point is not on the interior then
+        // we may be able to get a more accurate result by projecting onto the
+        // lines.
+        if (structdim == 2)
+          {
+            std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
+              line_projections;
+            for (unsigned int line_n = 0;
+                 line_n < GeometryInfo<structdim>::lines_per_cell;
+                 ++line_n)
+              {
+                line_projections[line_n] =
+                  project_to_object(object->line(line_n), trial_point);
+              }
+            std::sort(line_projections.begin(),
+                      line_projections.end(),
+                      [&](const Point<spacedim> &a, const Point<spacedim> &b) {
+                        return a.distance(trial_point) <
+                               b.distance(trial_point);
+                      });
+            if (line_projections[0].distance(trial_point) <
+                projected_point.distance(trial_point))
+              projected_point = line_projections[0];
+          }
+      }
+    else
+      {
+        Assert(false, ExcNotImplemented());
+        return projected_point;
+      }
+
+    return projected_point;
+  }
+} // namespace GridTools
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 3d70ba8095710f43f9019e8eb33da461829b2e92..4c31e658aaaa38ee416429ce81eeb125a177b73e 100644 (file)
@@ -52,9 +52,10 @@ set(_separate_src
   grid_generator_pipe_junction.cc
   grid_in.cc
   grid_out.cc
-  grid_tools_cache.cc
   grid_tools.cc
+  grid_tools_cache.cc
   grid_tools_dof_handlers.cc
+  grid_tools_geometry.cc
   grid_tools_nontemplates.cc
   tria.cc
   )
@@ -78,8 +79,9 @@ set(_inst
   grid_out.inst.in
   grid_refinement.inst.in
   grid_tools.inst.in
-  grid_tools_dof_handlers.inst.in
   grid_tools_cache.inst.in
+  grid_tools_dof_handlers.inst.in
+  grid_tools_geometry.inst.in
   intergrid_map.inst.in
   manifold.inst.in
   manifold_lib.inst.in
index 6abda211d979adf9d27935bd01065536bd0d5d5b..ba458f8fa2184c6170c202a5eb0df9d083d79440 100644 (file)
@@ -106,338 +106,6 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace GridTools
 {
-  template <int dim, int spacedim>
-  double
-  diameter(const Triangulation<dim, spacedim> &tria)
-  {
-    // we can't deal with distributed meshes since we don't have all
-    // vertices locally. there is one exception, however: if the mesh has
-    // never been refined. the way to test this is not to ask
-    // tria.n_levels()==1, since this is something that can happen on one
-    // processor without being true on all. however, we can ask for the
-    // global number of active cells and use that
-#if defined(DEAL_II_WITH_P4EST) && defined(DEBUG)
-    if (const parallel::distributed::Triangulation<dim, spacedim> *p_tria =
-          dynamic_cast<
-            const parallel::distributed::Triangulation<dim, spacedim> *>(&tria))
-      Assert(p_tria->n_global_active_cells() == tria.n_cells(0),
-             ExcNotImplemented());
-#endif
-
-    // the algorithm used simply traverses all cells and picks out the
-    // boundary vertices. it may or may not be faster to simply get all
-    // vectors, don't mark boundary vertices, and compute the distances
-    // thereof, but at least as the mesh is refined, it seems better to
-    // first mark boundary nodes, as marking is O(N) in the number of
-    // cells/vertices, while computing the maximal distance is O(N*N)
-    const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
-    std::vector<bool> boundary_vertices(vertices.size(), false);
-
-    typename Triangulation<dim, spacedim>::active_cell_iterator cell =
-      tria.begin_active();
-    const typename Triangulation<dim, spacedim>::active_cell_iterator endc =
-      tria.end();
-    for (; cell != endc; ++cell)
-      for (const unsigned int face : cell->face_indices())
-        if (cell->face(face)->at_boundary())
-          for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
-            boundary_vertices[cell->face(face)->vertex_index(i)] = true;
-
-    // now traverse the list of boundary vertices and check distances.
-    // since distances are symmetric, we only have to check one half
-    double                            max_distance_sqr = 0;
-    std::vector<bool>::const_iterator pi = boundary_vertices.begin();
-    const unsigned int                N  = boundary_vertices.size();
-    for (unsigned int i = 0; i < N; ++i, ++pi)
-      {
-        std::vector<bool>::const_iterator pj = pi + 1;
-        for (unsigned int j = i + 1; j < N; ++j, ++pj)
-          if ((*pi == true) && (*pj == true) &&
-              ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
-            max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
-      }
-
-    return std::sqrt(max_distance_sqr);
-  }
-
-
-
-  template <int dim, int spacedim>
-  double
-  volume(const Triangulation<dim, spacedim> &triangulation)
-  {
-    Assert(triangulation.get_reference_cells().size() == 1,
-           ExcNotImplemented());
-    const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
-    return volume(
-      triangulation,
-      reference_cell.template get_default_linear_mapping<dim, spacedim>());
-  }
-
-
-
-  template <int dim, int spacedim>
-  double
-  volume(const Triangulation<dim, spacedim> &triangulation,
-         const Mapping<dim, spacedim>       &mapping)
-  {
-    // get the degree of the mapping if possible. if not, just assume 1
-    unsigned int mapping_degree = 1;
-    if (const auto *p = dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
-      mapping_degree = p->get_degree();
-    else if (const auto *p =
-               dynamic_cast<const MappingFE<dim, spacedim> *>(&mapping))
-      mapping_degree = p->get_degree();
-
-    // then initialize an appropriate quadrature formula
-    Assert(triangulation.get_reference_cells().size() == 1,
-           ExcNotImplemented());
-    const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
-    const Quadrature<dim> quadrature_formula =
-      reference_cell.template get_gauss_type_quadrature<dim>(mapping_degree +
-                                                             1);
-    const unsigned int n_q_points = quadrature_formula.size();
-
-    // we really want the JxW values from the FEValues object, but it
-    // wants a finite element. create a cheap element as a dummy
-    // element
-    FE_Nothing<dim, spacedim> dummy_fe(reference_cell);
-    FEValues<dim, spacedim>   fe_values(mapping,
-                                      dummy_fe,
-                                      quadrature_formula,
-                                      update_JxW_values);
-
-    double local_volume = 0;
-
-    // compute the integral quantities by quadrature
-    for (const auto &cell : triangulation.active_cell_iterators())
-      if (cell->is_locally_owned())
-        {
-          fe_values.reinit(cell);
-          for (unsigned int q = 0; q < n_q_points; ++q)
-            local_volume += fe_values.JxW(q);
-        }
-
-    const double global_volume =
-      Utilities::MPI::sum(local_volume, triangulation.get_communicator());
-
-    return global_volume;
-  }
-
-
-
-  namespace
-  {
-    /**
-     * The algorithm to compute the affine approximation to a cell finds an
-     * affine map A x_hat + b from the reference cell to the real space.
-     *
-     * Some details about how we compute the least square plane. We look
-     * for a spacedim x (dim + 1) matrix X such that X * M = Y where M is
-     * a (dim+1) x n_vertices matrix and Y a spacedim x n_vertices.  And:
-     * The i-th column of M is unit_vertex[i] and the last row all
-     * 1's. The i-th column of Y is real_vertex[i].  If we split X=[A|b],
-     * the least square approx is A x_hat+b Classically X = Y * (M^t (M
-     * M^t)^{-1}) Let K = M^t * (M M^t)^{-1} = [KA Kb] this can be
-     * precomputed, and that is exactly what we do.  Finally A = Y*KA and
-     * b = Y*Kb.
-     */
-    template <int dim>
-    struct TransformR2UAffine
-    {
-      static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
-      static const double Kb[GeometryInfo<dim>::vertices_per_cell];
-    };
-
-
-    /*
-      Octave code:
-      M=[0 1; 1 1];
-      K1 = transpose(M) * inverse (M*transpose(M));
-      printf ("{%f, %f},\n", K1' );
-    */
-    template <>
-    const double TransformR2UAffine<1>::KA[GeometryInfo<1>::vertices_per_cell]
-                                          [1] = {{-1.000000}, {1.000000}};
-
-    template <>
-    const double TransformR2UAffine<1>::Kb[GeometryInfo<1>::vertices_per_cell] =
-      {1.000000, 0.000000};
-
-
-    /*
-      Octave code:
-      M=[0 1 0 1;0 0 1 1;1 1 1 1];
-      K2 = transpose(M) * inverse (M*transpose(M));
-      printf ("{%f, %f, %f},\n", K2' );
-    */
-    template <>
-    const double TransformR2UAffine<2>::KA[GeometryInfo<2>::vertices_per_cell]
-                                          [2] = {{-0.500000, -0.500000},
-                                                 {0.500000, -0.500000},
-                                                 {-0.500000, 0.500000},
-                                                 {0.500000, 0.500000}};
-
-    /*
-      Octave code:
-      M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
-      K3 = transpose(M) * inverse (M*transpose(M))
-      printf ("{%f, %f, %f, %f},\n", K3' );
-    */
-    template <>
-    const double TransformR2UAffine<2>::Kb[GeometryInfo<2>::vertices_per_cell] =
-      {0.750000, 0.250000, 0.250000, -0.250000};
-
-
-    template <>
-    const double TransformR2UAffine<3>::KA[GeometryInfo<3>::vertices_per_cell]
-                                          [3] = {
-                                            {-0.250000, -0.250000, -0.250000},
-                                            {0.250000, -0.250000, -0.250000},
-                                            {-0.250000, 0.250000, -0.250000},
-                                            {0.250000, 0.250000, -0.250000},
-                                            {-0.250000, -0.250000, 0.250000},
-                                            {0.250000, -0.250000, 0.250000},
-                                            {-0.250000, 0.250000, 0.250000},
-                                            {0.250000, 0.250000, 0.250000}
-
-    };
-
-
-    template <>
-    const double TransformR2UAffine<3>::Kb[GeometryInfo<3>::vertices_per_cell] =
-      {0.500000,
-       0.250000,
-       0.250000,
-       0.000000,
-       0.250000,
-       0.000000,
-       0.000000,
-       -0.250000};
-  } // namespace
-
-
-
-  template <int dim, int spacedim>
-  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
-  affine_cell_approximation(const ArrayView<const Point<spacedim>> &vertices)
-  {
-    AssertDimension(vertices.size(), GeometryInfo<dim>::vertices_per_cell);
-
-    // A = vertex * KA
-    DerivativeForm<1, dim, spacedim> A;
-
-    for (unsigned int d = 0; d < spacedim; ++d)
-      for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
-        for (unsigned int e = 0; e < dim; ++e)
-          A[d][e] += vertices[v][d] * TransformR2UAffine<dim>::KA[v][e];
-
-    // b = vertex * Kb
-    Tensor<1, spacedim> b;
-    for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
-      b += vertices[v] * TransformR2UAffine<dim>::Kb[v];
-
-    return std::make_pair(A, b);
-  }
-
-
-
-  template <int dim>
-  Vector<double>
-  compute_aspect_ratio_of_cells(const Mapping<dim>       &mapping,
-                                const Triangulation<dim> &triangulation,
-                                const Quadrature<dim>    &quadrature)
-  {
-    FE_Nothing<dim> fe;
-    FEValues<dim>   fe_values(mapping, fe, quadrature, update_jacobians);
-
-    Vector<double> aspect_ratio_vector(triangulation.n_active_cells());
-
-    // loop over cells of processor
-    for (const auto &cell : triangulation.active_cell_iterators())
-      {
-        if (cell->is_locally_owned())
-          {
-            double aspect_ratio_cell = 0.0;
-
-            fe_values.reinit(cell);
-
-            // loop over quadrature points
-            for (unsigned int q = 0; q < quadrature.size(); ++q)
-              {
-                const Tensor<2, dim, double> jacobian =
-                  Tensor<2, dim, double>(fe_values.jacobian(q));
-
-                // We intentionally do not want to throw an exception in case of
-                // inverted elements since this is not the task of this
-                // function. Instead, inf is written into the vector in case of
-                // inverted elements.
-                if (determinant(jacobian) <= 0)
-                  {
-                    aspect_ratio_cell = std::numeric_limits<double>::infinity();
-                  }
-                else
-                  {
-                    LAPACKFullMatrix<double> J = LAPACKFullMatrix<double>(dim);
-                    for (unsigned int i = 0; i < dim; ++i)
-                      for (unsigned int j = 0; j < dim; ++j)
-                        J(i, j) = jacobian[i][j];
-
-                    J.compute_svd();
-
-                    const double max_sv = J.singular_value(0);
-                    const double min_sv = J.singular_value(dim - 1);
-                    const double ar     = max_sv / min_sv;
-
-                    // Take the max between the previous and the current
-                    // aspect ratio value; if we had previously encountered
-                    // an inverted cell, we will have placed an infinity
-                    // in the aspect_ratio_cell variable, and that value
-                    // will survive this max operation.
-                    aspect_ratio_cell = std::max(aspect_ratio_cell, ar);
-                  }
-              }
-
-            // fill vector
-            aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
-          }
-      }
-
-    return aspect_ratio_vector;
-  }
-
-
-
-  template <int dim>
-  double
-  compute_maximum_aspect_ratio(const Mapping<dim>       &mapping,
-                               const Triangulation<dim> &triangulation,
-                               const Quadrature<dim>    &quadrature)
-  {
-    Vector<double> aspect_ratio_vector =
-      compute_aspect_ratio_of_cells(mapping, triangulation, quadrature);
-
-    return VectorTools::compute_global_error(triangulation,
-                                             aspect_ratio_vector,
-                                             VectorTools::Linfty_norm);
-  }
-
-
-
-  template <int dim, int spacedim>
-  BoundingBox<spacedim>
-  compute_bounding_box(const Triangulation<dim, spacedim> &tria)
-  {
-    using iterator =
-      typename Triangulation<dim, spacedim>::active_cell_iterator;
-    const auto predicate = [](const iterator &) { return true; };
-
-    return compute_bounding_box(
-      tria, std::function<bool(const iterator &)>(predicate));
-  }
-
-
-
   // Generic functions for appending face data in 2d or 3d. TODO: we can
   // remove these once we have 'if constexpr'.
   namespace internal
@@ -4442,40 +4110,6 @@ namespace GridTools
 
 
 
-  template <int dim, int spacedim>
-  double
-  minimal_cell_diameter(const Triangulation<dim, spacedim> &triangulation,
-                        const Mapping<dim, spacedim>       &mapping)
-  {
-    double min_diameter = std::numeric_limits<double>::max();
-    for (const auto &cell : triangulation.active_cell_iterators())
-      if (!cell->is_artificial())
-        min_diameter = std::min(min_diameter, cell->diameter(mapping));
-
-    const double global_min_diameter =
-      Utilities::MPI::min(min_diameter, triangulation.get_communicator());
-    return global_min_diameter;
-  }
-
-
-
-  template <int dim, int spacedim>
-  double
-  maximal_cell_diameter(const Triangulation<dim, spacedim> &triangulation,
-                        const Mapping<dim, spacedim>       &mapping)
-  {
-    double max_diameter = 0.;
-    for (const auto &cell : triangulation.active_cell_iterators())
-      if (!cell->is_artificial())
-        max_diameter = std::max(max_diameter, cell->diameter(mapping));
-
-    const double global_max_diameter =
-      Utilities::MPI::max(max_diameter, triangulation.get_communicator());
-    return global_max_diameter;
-  }
-
-
-
   namespace internal
   {
     namespace FixUpDistortedChildCells
index 446fd33c27788492f658ca879f63434e999ec5df..4824f738883d84d027a21bf6a36bcb53d718c6b2 100644 (file)
@@ -255,16 +255,6 @@ for (deal_II_space_dimension : SPACE_DIMENSIONS)
     GridTools::build_global_description_tree(
       const std::vector<BoundingBox<deal_II_space_dimension>> &,
       const MPI_Comm);
-
-    template Vector<double> GridTools::compute_aspect_ratio_of_cells(
-      const Mapping<deal_II_space_dimension> &,
-      const Triangulation<deal_II_space_dimension> &,
-      const Quadrature<deal_II_space_dimension> &);
-
-    template double GridTools::compute_maximum_aspect_ratio(
-      const Mapping<deal_II_space_dimension> &,
-      const Triangulation<deal_II_space_dimension> &,
-      const Quadrature<deal_II_space_dimension> &);
   }
 
 
@@ -273,27 +263,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 #if deal_II_dimension <= deal_II_space_dimension
     namespace GridTools
     \{
-      template double
-      diameter(
-        const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
-
-      template double
-      volume(const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
-
-      template double
-      volume(const Triangulation<deal_II_dimension, deal_II_space_dimension> &,
-             const Mapping<deal_II_dimension, deal_II_space_dimension> &);
-
-      template std::pair<
-        DerivativeForm<1, deal_II_dimension, deal_II_space_dimension>,
-        Tensor<1, deal_II_space_dimension>>
-      affine_cell_approximation<deal_II_dimension, deal_II_space_dimension>(
-        const ArrayView<const Point<deal_II_space_dimension>> &);
-
-      template BoundingBox<deal_II_space_dimension>
-      compute_bounding_box(
-        const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
-
       template std::tuple<std::vector<Point<deal_II_space_dimension>>,
                           std::vector<CellData<deal_II_dimension>>,
                           SubCellData>
@@ -441,18 +410,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
       get_locally_owned_vertices(
         const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
 
-      template double
-      minimal_cell_diameter(
-        const Triangulation<deal_II_dimension, deal_II_space_dimension>
-          &triangulation,
-        const Mapping<deal_II_dimension, deal_II_space_dimension> &);
-
-      template double
-      maximal_cell_diameter(
-        const Triangulation<deal_II_dimension, deal_II_space_dimension>
-          &triangulation,
-        const Mapping<deal_II_dimension, deal_II_space_dimension> &);
-
       template std::map<unsigned int, Point<deal_II_space_dimension>>
       get_all_vertices_at_boundary(
         const Triangulation<deal_II_dimension, deal_II_space_dimension> &tria);
diff --git a/source/grid/grid_tools_geometry.cc b/source/grid/grid_tools_geometry.cc
new file mode 100644 (file)
index 0000000..fd20d88
--- /dev/null
@@ -0,0 +1,410 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/derivative_form.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria_base.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_tools_geometry.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/lapack_full_matrix.h>
+
+#include <deal.II/numerics/vector_tools_integrate_difference.h>
+
+#include <functional>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace GridTools
+{
+  template <int dim, int spacedim>
+  double
+  diameter(const Triangulation<dim, spacedim> &tria)
+  {
+    // we can't deal with distributed meshes since we don't have all
+    // vertices locally. there is one exception, however: if the mesh has
+    // never been refined. the way to test this is not to ask
+    // tria.n_levels()==1, since this is something that can happen on one
+    // processor without being true on all. however, we can ask for the
+    // global number of active cells and use that
+#ifdef DEBUG
+    if (const auto *p_tria = dynamic_cast<
+          const parallel::DistributedTriangulationBase<dim, spacedim> *>(&tria))
+      Assert(p_tria->n_global_active_cells() == tria.n_cells(0),
+             ExcNotImplemented());
+#endif
+
+    // the algorithm used simply traverses all cells and picks out the
+    // boundary vertices. it may or may not be faster to simply get all
+    // vectors, don't mark boundary vertices, and compute the distances
+    // thereof, but at least as the mesh is refined, it seems better to
+    // first mark boundary nodes, as marking is O(N) in the number of
+    // cells/vertices, while computing the maximal distance is O(N*N)
+    const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
+    std::vector<bool> boundary_vertices(vertices.size(), false);
+
+    typename Triangulation<dim, spacedim>::active_cell_iterator cell =
+      tria.begin_active();
+    const typename Triangulation<dim, spacedim>::active_cell_iterator endc =
+      tria.end();
+    for (; cell != endc; ++cell)
+      for (const unsigned int face : cell->face_indices())
+        if (cell->face(face)->at_boundary())
+          for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
+            boundary_vertices[cell->face(face)->vertex_index(i)] = true;
+
+    // now traverse the list of boundary vertices and check distances.
+    // since distances are symmetric, we only have to check one half
+    double                            max_distance_sqr = 0;
+    std::vector<bool>::const_iterator pi = boundary_vertices.begin();
+    const unsigned int                N  = boundary_vertices.size();
+    for (unsigned int i = 0; i < N; ++i, ++pi)
+      {
+        std::vector<bool>::const_iterator pj = pi + 1;
+        for (unsigned int j = i + 1; j < N; ++j, ++pj)
+          if ((*pi == true) && (*pj == true) &&
+              ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
+            max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
+      }
+
+    return std::sqrt(max_distance_sqr);
+  }
+
+
+
+  template <int dim, int spacedim>
+  double
+  volume(const Triangulation<dim, spacedim> &triangulation)
+  {
+    Assert(triangulation.get_reference_cells().size() == 1,
+           ExcNotImplemented());
+    const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
+    return volume(
+      triangulation,
+      reference_cell.template get_default_linear_mapping<dim, spacedim>());
+  }
+
+
+
+  template <int dim, int spacedim>
+  double
+  volume(const Triangulation<dim, spacedim> &triangulation,
+         const Mapping<dim, spacedim>       &mapping)
+  {
+    // get the degree of the mapping if possible. if not, just assume 1
+    unsigned int mapping_degree = 1;
+    if (const auto *p = dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
+      mapping_degree = p->get_degree();
+    else if (const auto *p =
+               dynamic_cast<const MappingFE<dim, spacedim> *>(&mapping))
+      mapping_degree = p->get_degree();
+
+    // then initialize an appropriate quadrature formula
+    Assert(triangulation.get_reference_cells().size() == 1,
+           ExcNotImplemented());
+    const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
+    const Quadrature<dim> quadrature_formula =
+      reference_cell.template get_gauss_type_quadrature<dim>(mapping_degree +
+                                                             1);
+    const unsigned int n_q_points = quadrature_formula.size();
+
+    // we really want the JxW values from the FEValues object, but it
+    // wants a finite element. create a cheap element as a dummy
+    // element
+    FE_Nothing<dim, spacedim> dummy_fe(reference_cell);
+    FEValues<dim, spacedim>   fe_values(mapping,
+                                      dummy_fe,
+                                      quadrature_formula,
+                                      update_JxW_values);
+
+    double local_volume = 0;
+
+    // compute the integral quantities by quadrature
+    for (const auto &cell : triangulation.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          fe_values.reinit(cell);
+          for (unsigned int q = 0; q < n_q_points; ++q)
+            local_volume += fe_values.JxW(q);
+        }
+
+    const double global_volume =
+      Utilities::MPI::sum(local_volume, triangulation.get_communicator());
+
+    return global_volume;
+  }
+
+
+
+  namespace
+  {
+    /**
+     * The algorithm to compute the affine approximation to a cell finds an
+     * affine map A x_hat + b from the reference cell to the real space.
+     *
+     * Some details about how we compute the least square plane. We look
+     * for a spacedim x (dim + 1) matrix X such that X * M = Y where M is
+     * a (dim+1) x n_vertices matrix and Y a spacedim x n_vertices.  And:
+     * The i-th column of M is unit_vertex[i] and the last row all
+     * 1's. The i-th column of Y is real_vertex[i].  If we split X=[A|b],
+     * the least square approx is A x_hat+b Classically X = Y * (M^t (M
+     * M^t)^{-1}) Let K = M^t * (M M^t)^{-1} = [KA Kb] this can be
+     * precomputed, and that is exactly what we do.  Finally A = Y*KA and
+     * b = Y*Kb.
+     */
+    template <int dim>
+    struct TransformR2UAffine
+    {
+      static const double KA[GeometryInfo<dim>::vertices_per_cell][dim];
+      static const double Kb[GeometryInfo<dim>::vertices_per_cell];
+    };
+
+
+    /*
+      Octave code:
+      M=[0 1; 1 1];
+      K1 = transpose(M) * inverse (M*transpose(M));
+      printf ("{%f, %f},\n", K1' );
+    */
+    template <>
+    const double TransformR2UAffine<1>::KA[GeometryInfo<1>::vertices_per_cell]
+                                          [1] = {{-1.000000}, {1.000000}};
+
+    template <>
+    const double TransformR2UAffine<1>::Kb[GeometryInfo<1>::vertices_per_cell] =
+      {1.000000, 0.000000};
+
+
+    /*
+      Octave code:
+      M=[0 1 0 1;0 0 1 1;1 1 1 1];
+      K2 = transpose(M) * inverse (M*transpose(M));
+      printf ("{%f, %f, %f},\n", K2' );
+    */
+    template <>
+    const double TransformR2UAffine<2>::KA[GeometryInfo<2>::vertices_per_cell]
+                                          [2] = {{-0.500000, -0.500000},
+                                                 {0.500000, -0.500000},
+                                                 {-0.500000, 0.500000},
+                                                 {0.500000, 0.500000}};
+
+    /*
+      Octave code:
+      M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1];
+      K3 = transpose(M) * inverse (M*transpose(M))
+      printf ("{%f, %f, %f, %f},\n", K3' );
+    */
+    template <>
+    const double TransformR2UAffine<2>::Kb[GeometryInfo<2>::vertices_per_cell] =
+      {0.750000, 0.250000, 0.250000, -0.250000};
+
+
+    template <>
+    const double TransformR2UAffine<3>::KA[GeometryInfo<3>::vertices_per_cell]
+                                          [3] = {
+                                            {-0.250000, -0.250000, -0.250000},
+                                            {0.250000, -0.250000, -0.250000},
+                                            {-0.250000, 0.250000, -0.250000},
+                                            {0.250000, 0.250000, -0.250000},
+                                            {-0.250000, -0.250000, 0.250000},
+                                            {0.250000, -0.250000, 0.250000},
+                                            {-0.250000, 0.250000, 0.250000},
+                                            {0.250000, 0.250000, 0.250000}
+
+    };
+
+
+    template <>
+    const double TransformR2UAffine<3>::Kb[GeometryInfo<3>::vertices_per_cell] =
+      {0.500000,
+       0.250000,
+       0.250000,
+       0.000000,
+       0.250000,
+       0.000000,
+       0.000000,
+       -0.250000};
+  } // namespace
+
+
+
+  template <int dim, int spacedim>
+  std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
+  affine_cell_approximation(const ArrayView<const Point<spacedim>> &vertices)
+  {
+    AssertDimension(vertices.size(), GeometryInfo<dim>::vertices_per_cell);
+
+    // A = vertex * KA
+    DerivativeForm<1, dim, spacedim> A;
+
+    for (unsigned int d = 0; d < spacedim; ++d)
+      for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
+        for (unsigned int e = 0; e < dim; ++e)
+          A[d][e] += vertices[v][d] * TransformR2UAffine<dim>::KA[v][e];
+
+    // b = vertex * Kb
+    Tensor<1, spacedim> b;
+    for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
+      b += vertices[v] * TransformR2UAffine<dim>::Kb[v];
+
+    return std::make_pair(A, b);
+  }
+
+
+
+  template <int dim>
+  Vector<double>
+  compute_aspect_ratio_of_cells(const Mapping<dim>       &mapping,
+                                const Triangulation<dim> &triangulation,
+                                const Quadrature<dim>    &quadrature)
+  {
+    FE_Nothing<dim> fe;
+    FEValues<dim>   fe_values(mapping, fe, quadrature, update_jacobians);
+
+    Vector<double> aspect_ratio_vector(triangulation.n_active_cells());
+
+    // loop over cells of processor
+    for (const auto &cell : triangulation.active_cell_iterators())
+      {
+        if (cell->is_locally_owned())
+          {
+            double aspect_ratio_cell = 0.0;
+
+            fe_values.reinit(cell);
+
+            // loop over quadrature points
+            for (unsigned int q = 0; q < quadrature.size(); ++q)
+              {
+                const Tensor<2, dim, double> jacobian =
+                  Tensor<2, dim, double>(fe_values.jacobian(q));
+
+                // We intentionally do not want to throw an exception in case of
+                // inverted elements since this is not the task of this
+                // function. Instead, inf is written into the vector in case of
+                // inverted elements.
+                if (determinant(jacobian) <= 0)
+                  {
+                    aspect_ratio_cell = std::numeric_limits<double>::infinity();
+                  }
+                else
+                  {
+                    LAPACKFullMatrix<double> J = LAPACKFullMatrix<double>(dim);
+                    for (unsigned int i = 0; i < dim; ++i)
+                      for (unsigned int j = 0; j < dim; ++j)
+                        J(i, j) = jacobian[i][j];
+
+                    J.compute_svd();
+
+                    const double max_sv = J.singular_value(0);
+                    const double min_sv = J.singular_value(dim - 1);
+                    const double ar     = max_sv / min_sv;
+
+                    // Take the max between the previous and the current
+                    // aspect ratio value; if we had previously encountered
+                    // an inverted cell, we will have placed an infinity
+                    // in the aspect_ratio_cell variable, and that value
+                    // will survive this max operation.
+                    aspect_ratio_cell = std::max(aspect_ratio_cell, ar);
+                  }
+              }
+
+            // fill vector
+            aspect_ratio_vector(cell->active_cell_index()) = aspect_ratio_cell;
+          }
+      }
+
+    return aspect_ratio_vector;
+  }
+
+
+
+  template <int dim>
+  double
+  compute_maximum_aspect_ratio(const Mapping<dim>       &mapping,
+                               const Triangulation<dim> &triangulation,
+                               const Quadrature<dim>    &quadrature)
+  {
+    Vector<double> aspect_ratio_vector =
+      compute_aspect_ratio_of_cells(mapping, triangulation, quadrature);
+
+    return VectorTools::compute_global_error(triangulation,
+                                             aspect_ratio_vector,
+                                             VectorTools::Linfty_norm);
+  }
+
+
+
+  template <int dim, int spacedim>
+  BoundingBox<spacedim>
+  compute_bounding_box(const Triangulation<dim, spacedim> &tria)
+  {
+    using iterator =
+      typename Triangulation<dim, spacedim>::active_cell_iterator;
+    const auto predicate = [](const iterator &) { return true; };
+
+    return compute_bounding_box(
+      tria, std::function<bool(const iterator &)>(predicate));
+  }
+
+
+
+  template <int dim, int spacedim>
+  double
+  minimal_cell_diameter(const Triangulation<dim, spacedim> &triangulation,
+                        const Mapping<dim, spacedim>       &mapping)
+  {
+    double min_diameter = std::numeric_limits<double>::max();
+    for (const auto &cell : triangulation.active_cell_iterators())
+      if (!cell->is_artificial())
+        min_diameter = std::min(min_diameter, cell->diameter(mapping));
+
+    const double global_min_diameter =
+      Utilities::MPI::min(min_diameter, triangulation.get_communicator());
+    return global_min_diameter;
+  }
+
+
+
+  template <int dim, int spacedim>
+  double
+  maximal_cell_diameter(const Triangulation<dim, spacedim> &triangulation,
+                        const Mapping<dim, spacedim>       &mapping)
+  {
+    double max_diameter = 0.;
+    for (const auto &cell : triangulation.active_cell_iterators())
+      if (!cell->is_artificial())
+        max_diameter = std::max(max_diameter, cell->diameter(mapping));
+
+    const double global_max_diameter =
+      Utilities::MPI::max(max_diameter, triangulation.get_communicator());
+    return global_max_diameter;
+  }
+} /* namespace GridTools */
+
+
+// explicit instantiations
+#include "grid_tools_geometry.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/grid/grid_tools_geometry.inst.in b/source/grid/grid_tools_geometry.inst.in
new file mode 100644 (file)
index 0000000..5695304
--- /dev/null
@@ -0,0 +1,73 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
+  {
+#if deal_II_dimension <= deal_II_space_dimension
+    namespace GridTools
+    \{
+      template double
+      diameter(
+        const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
+
+      template double
+      volume(const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
+
+      template double
+      volume(const Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+             const Mapping<deal_II_dimension, deal_II_space_dimension> &);
+
+      template double
+      minimal_cell_diameter(
+        const Triangulation<deal_II_dimension, deal_II_space_dimension>
+          &triangulation,
+        const Mapping<deal_II_dimension, deal_II_space_dimension> &);
+
+      template double
+      maximal_cell_diameter(
+        const Triangulation<deal_II_dimension, deal_II_space_dimension>
+          &triangulation,
+        const Mapping<deal_II_dimension, deal_II_space_dimension> &);
+
+      template std::pair<
+        DerivativeForm<1, deal_II_dimension, deal_II_space_dimension>,
+        Tensor<1, deal_II_space_dimension>>
+      affine_cell_approximation<deal_II_dimension, deal_II_space_dimension>(
+        const ArrayView<const Point<deal_II_space_dimension>> &);
+
+      template BoundingBox<deal_II_space_dimension>
+      compute_bounding_box(
+        const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
+    \}
+#endif
+  }
+
+for (deal_II_space_dimension : SPACE_DIMENSIONS)
+  {
+    namespace GridTools
+    \{
+      template Vector<double>
+      GridTools::compute_aspect_ratio_of_cells(
+        const Mapping<deal_II_space_dimension> &,
+        const Triangulation<deal_II_space_dimension> &,
+        const Quadrature<deal_II_space_dimension> &);
+
+      template double
+      GridTools::compute_maximum_aspect_ratio(
+        const Mapping<deal_II_space_dimension> &,
+        const Triangulation<deal_II_space_dimension> &,
+        const Quadrature<deal_II_space_dimension> &);
+    \}
+  }
index 8723a402affee08077ac7229c169ae0ef679d750..0836ac904c70521c0e20569c3dd4ce0c82169715 100644 (file)
@@ -21,7 +21,7 @@
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/mapping.h>
 
-#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_geometry.h>
 #include <deal.II/grid/manifold.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_accessor.h>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.