]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Overhaul GridTools::internal::fix_up_object.
authorDavid Wells <wellsd2@rpi.edu>
Sat, 12 Aug 2017 00:13:32 +0000 (20:13 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Sun, 8 Oct 2017 23:22:35 +0000 (19:22 -0400)
The primary issue with excising Boundary from the library is the lack of
an equivalent project_to_surface function in the Manifold code (since
Manifolds don't know anything about faces or lines, just points and
geodesics). Fortunately, this function is only called in one place, and
the default implementation in StraightManifold is the only
implementation in the library: hence, we get around this issue by
copying and pasting StraightManifold::project_to_surface into the new
internal project_to_object function and marking this change as
incompatible.

To the best of the author's knowledge, no users have implemented their
own custom project_to_surface nor are they using the one function that
calls this (GridTools::fix_up_distorted_child_cells) so this small
incompatibility is acceptable.

include/deal.II/grid/grid_tools.h
source/fe/mapping_c1.cc
source/grid/grid_tools.cc
tests/grid/project_to_object_01.cc [new file with mode: 0644]
tests/grid/project_to_object_01.output [new file with mode: 0644]

index 7472f190804ef8257075139e529e4a61d2c70a52..a5772ac992d5310552b30fcccd4648011d1d6c42 100644 (file)
@@ -23,6 +23,7 @@
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/fe/mapping.h>
 #include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/manifold.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
@@ -187,6 +188,30 @@ namespace GridTools
   template <int dim, int spacedim>
   BoundingBox<spacedim> compute_bounding_box(const Triangulation<dim, spacedim> &triangulation);
 
+  /**
+   * Return the point on the geometrical object @object closest to the given
+   * point @p trial_point. For example, if @p object is a one-dimensional line
+   * or edge, then the 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.
+   *
+   * @author Luca Heltai, David Wells, 2017.
+   */
+  template <typename Iterator>
+  Point<Iterator::AccessorType::space_dimension>
+  project_to_object(const Iterator &object,
+                    const Point<Iterator::AccessorType::space_dimension> &trial_point);
+
   /*@}*/
   /**
    * @name Functions supporting the creation of meshes
@@ -2253,6 +2278,512 @@ 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 quad (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 quad, and \xi=(\xi_1,\xi_2) are the
+        // reference coordinates of the quad. 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 (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
+          x_k += object->vertex(i) *
+                 GeometryInfo<structdim>::d_linear_shape_function (xi, i);
+
+        do
+          {
+            Tensor<1,structdim> F_k;
+            for (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
+              F_k += (x_k-trial_point)*object->vertex(i) *
+                     GeometryInfo<structdim>::d_linear_shape_function_gradient (xi, i);
+
+            Tensor<2,structdim> H_k;
+            for (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
+              for (unsigned int j=0; j<GeometryInfo<structdim>::vertices_per_cell; ++j)
+                {
+                  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 (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
+              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;
+      }
+    }
+  }
+
+
+
+  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.
+            auto weights_are_ok = [](const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &v)
+                                  -> bool
+            {
+              // 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
+            };
+            const double step_size = object->diameter()/64.0;
+
+
+            std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> vertices;
+            for (unsigned int vertex_n = 0; vertex_n < GeometryInfo<structdim>::vertices_per_cell;
+                 ++vertex_n)
+              vertices[vertex_n] = object->vertex(vertex_n);
+
+            auto get_point_from_weights =
+              [&](const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &weights)
+              -> Point<spacedim>
+            {
+              return object->get_manifold().get_new_point
+              (make_array_view(vertices.begin(), vertices.end()),
+              make_array_view(&weights[0],
+              &weights[GeometryInfo<structdim>::vertices_per_cell - 1] + 1));
+            };
+
+            // pick the initial weights as (normalized) inverse distances from
+            // the trial point:
+            Tensor<1, GeometryInfo<structdim>::vertices_per_cell> guess_weights;
+            double guess_weights_sum = 0.0;
+            for (unsigned int vertex_n = 0; vertex_n < GeometryInfo<structdim>::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(weights_are_ok(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 = GeometryInfo<structdim>::vertices_per_cell - 1;
+                Tensor<1, GeometryInfo<structdim>::vertices_per_cell> current_gradient;
+                for (unsigned int row_n = 0;
+                     row_n < GeometryInfo<structdim>::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, GeometryInfo<structdim>::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 (weights_are_ok(tentative_weights))
+                      break;
+
+                    for (unsigned int i = 0; i < GeometryInfo<structdim>::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 (!weights_are_ok(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(weights_are_ok(guess_weights), ExcInternalError());
+              }
+            Assert(weights_are_ok(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;
+  }
 }
 
 #endif
index f20482269e71977e1e72f8e4624a012b97afb5d2..cff8eb50a54ce2812bc3f867b34c3acea60835fe 100644 (file)
@@ -69,7 +69,6 @@ MappingC1<2>::MappingC1Generic::add_line_support_points (const Triangulation<2>:
   const std::array<double, 2> interior_gl_points
   {{0.5 - 0.5*std::sqrt(1.0/5.0), 0.5 + 0.5*std::sqrt(1.0/5.0)}};
 
-
   // loop over each of the lines, and if it is at the boundary, then first get
   // the boundary description and second compute the points on it. if not at
   // the boundary, get the respective points from another function
index 16c19984e65d67e200ff17df3a4fd9f24b403e99..c469d0e734e2c34211d54838c2f5f276ff6a42c7 100644 (file)
@@ -3269,75 +3269,45 @@ next_cell:
 
 
       /**
-       * Try to fix up a single cell. Return
-       * whether we succeeded with this.
-       *
-       * The second argument indicates
-       * whether we need to respect the
-       * manifold/boundary on which this
-       * object lies when moving around its
-       * mid-point.
+       * Try to fix up a single cell by moving around its midpoint. Return whether we succeeded with this.
        */
       template <typename Iterator>
       bool
-      fix_up_object (const Iterator &object,
-                     const bool respect_manifold)
+      fix_up_object (const Iterator &object)
       {
-        const Boundary<Iterator::AccessorType::dimension,
-              Iterator::AccessorType::space_dimension>
-              *manifold = (respect_manifold ?
-                           &object->get_boundary() :
-                           nullptr);
-
         const unsigned int structdim = Iterator::AccessorType::structure_dimension;
         const unsigned int spacedim  = Iterator::AccessorType::space_dimension;
 
-        // right now we can only deal
-        // with cells that have been
-        // refined isotropically
-        // because that is the only
-        // case where we have a cell
-        // mid-point that can be moved
-        // around without having to
-        // consider boundary
-        // information
+        // right now we can only deal with cells that have been refined
+        // isotropically because that is the only case where we have a cell
+        // mid-point that can be moved around without having to consider
+        // boundary information
         Assert (object->has_children(), ExcInternalError());
         Assert (object->refinement_case() == RefinementCase<structdim>::isotropic_refinement,
                 ExcNotImplemented());
 
-        // get the current location of
-        // the object mid-vertex:
+        // get the current location of the object mid-vertex:
         Point<spacedim> object_mid_point
           = object->child(0)->vertex (GeometryInfo<structdim>::max_children_per_cell-1);
 
-        // now do a few steepest descent
-        // steps to reduce the objective
-        // function. compute the diameter in
-        // the helper function above
+        // now do a few steepest descent steps to reduce the objective
+        // function. compute the diameter in the helper function above
         unsigned int iteration = 0;
         const double diameter = minimal_diameter (object);
 
-        // current value of objective
-        // function and initial delta
+        // current value of objective function and initial delta
         double current_value = objective_function (object, object_mid_point);
         double initial_delta = 0;
 
         do
           {
-            // choose a step length
-            // that is initially 1/4
-            // of the child objects'
-            // diameter, and a sequence
-            // whose sum does not
-            // converge (to avoid
-            // premature termination of
-            // the iteration)
+            // choose a step length that is initially 1/4 of the child
+            // objects' diameter, and a sequence whose sum does not converge
+            // (to avoid premature termination of the iteration)
             const double step_length = diameter / 4 / (iteration + 1);
 
-            // compute the objective
-            // function's derivative using a
-            // two-sided difference formula
-            // with eps=step_length/10
+            // compute the objective function's derivative using a two-sided
+            // difference formula with eps=step_length/10
             Tensor<1,spacedim> gradient;
             for (unsigned int d=0; d<spacedim; ++d)
               {
@@ -3346,68 +3316,35 @@ next_cell:
                 Tensor<1,spacedim> h;
                 h[d] = eps/2;
 
-                if (respect_manifold == false)
-                  gradient[d]
-                    = ((objective_function (object, object_mid_point + h)
-                        -
-                        objective_function (object, object_mid_point - h))
-                       /
-                       eps);
-                else
-                  gradient[d]
-                    = ((objective_function (object,
-                                            manifold->project_to_surface(object,
-                                                                         object_mid_point + h))
-                        -
-                        objective_function (object,
-                                            manifold->project_to_surface(object,
-                                                                         object_mid_point - h)))
-                       /
-                       eps);
+                gradient[d] = (objective_function (object,
+                                                   project_to_object(object, object_mid_point + h))
+                               -
+                               objective_function (object,
+                                                   project_to_object(object, object_mid_point - h)))
+                              /eps;
               }
 
-            // sometimes, the
-            // (unprojected) gradient
-            // is perpendicular to
-            // the manifold, but we
-            // can't go there if
-            // respect_manifold==true. in
-            // that case, gradient=0,
-            // and we simply need to
-            // quite the loop here
+            // there is nowhere to go
             if (gradient.norm() == 0)
               break;
 
-            // so we need to go in
-            // direction -gradient. the
-            // optimal value of the
-            // objective function is
-            // zero, so assuming that
-            // the model is quadratic
-            // we would have to go
-            // -2*val/||gradient|| in
-            // this direction, make
-            // sure we go at most
-            // step_length into this
+            // We need to go in direction -gradient. the optimal value of the
+            // objective function is zero, so assuming that the model is
+            // quadratic we would have to go -2*val/||gradient|| in this
+            // direction, make sure we go at most step_length into this
             // direction
             object_mid_point -= std::min(2 * current_value / (gradient*gradient),
-                                         step_length / gradient.norm()) *
-                                gradient;
+                                         step_length / gradient.norm()) * gradient;
+            object_mid_point = project_to_object(object, object_mid_point);
 
-            if (respect_manifold == true)
-              object_mid_point = manifold->project_to_surface(object,
-                                                              object_mid_point);
-
-            // compute current value of the
-            // objective function
+            // compute current value of the objective function
             const double previous_value = current_value;
             current_value = objective_function (object, object_mid_point);
 
             if (iteration == 0)
               initial_delta = (previous_value - current_value);
 
-            // stop if we aren't moving much
-            // any more
+            // stop if we aren't moving much any more
             if ((iteration >= 1) &&
                 ((previous_value - current_value < 0)
                  ||
@@ -3516,18 +3453,15 @@ next_cell:
                          std::integral_constant<int, 1>,
                          std::integral_constant<int, 1>)
       {
-        // nothing to do for the faces of
-        // cells in 1d
+        // nothing to do for the faces of cells in 1d
       }
 
 
 
-      // possibly fix up the faces of
-      // a cell by moving around its
-      // mid-points
-      template <int structdim, int spacedim>
-      void fix_up_faces (const typename dealii::Triangulation<structdim,spacedim>::cell_iterator &cell,
-                         std::integral_constant<int, structdim>,
+      // possibly fix up the faces of a cell by moving around its mid-points
+      template <int dim, int spacedim>
+      void fix_up_faces (const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
+                         std::integral_constant<int, dim>,
                          std::integral_constant<int, spacedim>)
       {
         // see if we first can fix up
@@ -3555,15 +3489,15 @@ next_cell:
         // ignored at the time; we
         // should then also be able to
         // ignore it this time as well
-        for (unsigned int f=0; f<GeometryInfo<structdim>::faces_per_cell; ++f)
+        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
           {
             Assert (cell->face(f)->has_children(), ExcInternalError());
             Assert (cell->face(f)->refinement_case() ==
-                    RefinementCase<structdim-1>::isotropic_refinement,
+                    RefinementCase<dim - 1>::isotropic_refinement,
                     ExcInternalError());
 
             bool subface_is_more_refined = false;
-            for (unsigned int g=0; g<GeometryInfo<structdim>::max_children_per_face; ++g)
+            for (unsigned int g=0; g<GeometryInfo<dim>::max_children_per_face; ++g)
               if (cell->face(f)->child(g)->has_children())
                 {
                   subface_is_more_refined = true;
@@ -3573,28 +3507,23 @@ next_cell:
             if (subface_is_more_refined == true)
               continue;
 
-            // so, now we finally know
-            // that we can do something
-            // about this face
-            fix_up_object (cell->face(f), cell->at_boundary(f));
+            // we finally know that we can do something about this face
+            fix_up_object (cell->face(f));
           }
       }
-
-
     } /* namespace FixUpDistortedChildCells */
   } /* namespace internal */
 
 
   template <int dim, int spacedim>
   typename Triangulation<dim,spacedim>::DistortedCellList
-
-  fix_up_distorted_child_cells (const typename Triangulation<dim,spacedim>::DistortedCellList &distorted_cells,
-                                Triangulation<dim,spacedim> &/*triangulation*/)
+  fix_up_distorted_child_cells
+  (const typename Triangulation<dim,spacedim>::DistortedCellList &distorted_cells,
+   Triangulation<dim,spacedim> &/*triangulation*/)
   {
     typename Triangulation<dim,spacedim>::DistortedCellList unfixable_subset;
 
-    // loop over all cells that we have
-    // to fix up
+    // loop over all cells that we have to fix up
     for (typename std::list<typename Triangulation<dim,spacedim>::cell_iterator>::const_iterator
          cell_ptr = distorted_cells.distorted_cells.begin();
          cell_ptr != distorted_cells.distorted_cells.end(); ++cell_ptr)
@@ -3611,14 +3540,8 @@ next_cell:
                         std::integral_constant<int, dim>(),
                         std::integral_constant<int, spacedim>());
 
-        // fix up the object. we need to
-        // respect the manifold if the cell is
-        // embedded in a higher dimensional
-        // space; otherwise, like a hex in 3d,
-        // every point within the cell interior
-        // is fair game
-        if (! internal::FixUpDistortedChildCells::fix_up_object (cell,
-                                                                 (dim < spacedim)))
+        // If possible, fix up the object.
+        if (!internal::FixUpDistortedChildCells::fix_up_object (cell))
           unfixable_subset.distorted_cells.push_back (cell);
       }
 
diff --git a/tests/grid/project_to_object_01.cc b/tests/grid/project_to_object_01.cc
new file mode 100644 (file)
index 0000000..33d3f64
--- /dev/null
@@ -0,0 +1,466 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <array>
+#include <cmath>
+
+#include "../tests.h"
+
+int main()
+{
+  using namespace dealii;
+  initlog();
+
+  // test the internal cross derivative function: the stencil should be order 2.
+  {
+    using namespace dealii::GridTools::internal::ProjectToObject;
+
+    constexpr int structdim = 2;
+    auto objective = [](const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &weights)
+    {
+      return std::sin(2.0*weights[0])*std::cos(3.0*weights[1]) + std::cos(weights[2]);
+    };
+    Tensor<1, GeometryInfo<structdim>::vertices_per_cell> c0;
+    for (unsigned int row_n = 0;
+         row_n < GeometryInfo<structdim>::vertices_per_cell;
+         ++row_n)
+      c0[row_n] = 1.0;
+
+    const std::array<double, 3> exact
+    {
+      {
+        2.0*std::cos(3.0) *std::cos(2.0) + 3.0*std::sin(3.0) *std::sin(2.0),
+        -2.0*std::cos(3.0) *std::cos(2.0) - 3.0*std::sin(3.0) *std::sin(2.0),
+        -3.0*std::sin(3.0) *std::sin(2.0) + std::sin(1.0)
+      }
+    };
+    const std::array<CrossDerivative, 3> cross_derivatives {{{0, 1}, {1, 0}, {1, 2}}};
+
+    for (unsigned int cross_derivative_n = 0; cross_derivative_n < cross_derivatives.size();
+         ++cross_derivative_n)
+      {
+        deallog << "testing cross derivative "
+                << cross_derivative_n
+                << std::endl;
+        for (unsigned int step_n = 5; step_n < 10; ++step_n)
+          {
+            deallog << cross_stencil<structdim>(cross_derivatives[cross_derivative_n],
+                                                c0,
+                                                std::pow(0.5, double(step_n)),
+                                                objective) - exact[cross_derivative_n]
+                    << std::endl;
+          }
+      }
+  }
+
+  // Test project_to_object in 2D with an annulus
+  {
+    PolarManifold<2> polar_manifold;
+    Triangulation<2> triangulation;
+    constexpr types::manifold_id polar_id = 42;
+    GridGenerator::hyper_shell(triangulation, Point<2>(), 1.0, 2.0);
+    triangulation.set_manifold(polar_id, polar_manifold);
+    triangulation.set_all_manifold_ids(42);
+
+    triangulation.refine_global(2);
+
+    auto cell = triangulation.begin_active();
+    while (!cell->at_boundary())
+      ++cell;
+    unsigned int face_n = 0;
+    while (!cell->face(face_n)->at_boundary())
+      ++face_n;
+    const auto face = cell->face(face_n);
+
+    // easy test: project the arithmetic mean of the vertices onto the
+    // geodesic
+    {
+      const Point<2> trial_point = 0.9*face->center();
+      const Point<2> projected_point = GridTools::project_to_object(face, trial_point);
+
+      const std::array<Point<2>, 2> vertices {{face->vertex(0), face->vertex(1)}};
+      const std::array<double, 2> weights {{0.5, 0.5}};
+      const auto vertices_view = make_array_view(vertices.begin(), vertices.end());
+      const auto weights_view = make_array_view(weights.begin(), weights.end());
+
+      deallog << std::endl
+              << "Project the arithmetic mean of two vertices on a circular"
+              << std::endl
+              << "arc. The result should be the geodesic midpoint:"
+              << std::endl
+              << std::endl
+              << "vertex 0 distance from origin: "
+              << face->vertex(0).distance(Point<2>()) << std::endl
+              << "vertex 1 distance from origin: "
+              << face->vertex(1).distance(Point<2>()) << std::endl
+              << "trial point distance from origin: "
+              << trial_point.distance(Point<2>()) << std::endl
+              << "projected point distance from origin: "
+              << projected_point.distance(Point<2>()) << std::endl
+              << "projected point:        "
+              << projected_point << std::endl
+              << "true geodesic midpoint: "
+              << face->get_manifold().get_new_point(vertices_view, weights_view)
+              << std::endl;
+    }
+
+    // test that we can project a convex combination point
+    {
+      const std::array<Point<2>, 2> vertices {{face->vertex(0), face->vertex(1)}};
+      const std::array<double, 2> weights {{0.125, 0.875}};
+      const auto vertices_view = make_array_view(vertices.begin(), vertices.end());
+      const auto weights_view = make_array_view(weights.begin(), weights.end());
+
+      const Point<2> trial_point = face->get_manifold().get_new_point(vertices_view, weights_view);
+      const Point<2> projected_point = GridTools::project_to_object(face, trial_point);
+      deallog << std::endl
+              << "Project the convex combination of two vertices:"
+              << std::endl
+              << "projected point:          "
+              << projected_point
+              << std::endl
+              << "Convex combination point: "
+              << trial_point
+              << std::endl
+              << "Distance:                 "
+              << trial_point.distance(projected_point)
+              << std::endl;
+    }
+
+    // easy test: project with spacedim == structdim (should be the identity)
+    {
+      const Point<2> trial_point {0.417941, 4242424242.4242};
+      const Point<2> projected_point = GridTools::project_to_object(cell,
+                                       trial_point);
+      deallog << std::endl
+              << "Check that projecting with spacedim == structdim"
+              << std::endl
+              << "is the identity map:"
+              << std::endl
+              << std::endl
+              << "trial point:     "
+              << trial_point
+              << std::endl
+              << "projected point: "
+              << projected_point
+              << std::endl;
+    }
+
+    // harder test: project a vertex onto the geodesic
+    {
+      const Point<2> trial_point = face->vertex(0);
+      const Point<2> projected_point = GridTools::project_to_object(cell->face(face_n),
+                                       trial_point);
+
+      deallog << std::endl
+              << "Project a vertex onto the geodesic:"
+              << std::endl
+              << std::endl
+              << "vertex distance from origin:          "
+              << face->vertex(1).distance(Point<2>())
+              << std::endl
+              << "trial point distance from origin:     "
+              << trial_point.distance(Point<2>())
+              << std::endl
+              << "projected point distance from origin: "
+              << projected_point.distance(Point<2>())
+              << std::endl
+              << "projected point: "
+              << projected_point << std::endl
+              << "actual vertex:   "
+              << face->vertex(0) << std::endl;
+    }
+  }
+
+  // Test project_to_object with a 2D surface in 3D
+  {
+    TorusManifold<2> torus_manifold(2.0, 1.0);
+    Triangulation<2, 3> triangulation;
+    GridGenerator::torus(triangulation, 2.0, 1.0);
+
+    constexpr types::manifold_id torus_id = 42;
+    triangulation.set_manifold(torus_id, torus_manifold);
+    triangulation.set_all_manifold_ids(torus_id);
+
+    // an ill-posed problem: project a point along the axis of symmetry onto a
+    // cell face. Make sure that we end up with something that is on the
+    // manifold and is near the equator (for the bottom cells) or on the top
+    // (for the top cells).
+    const Point<3> trial_point(0.0, 100.0, 0.0);
+    deallog << "Test for robustness by projecting points with nonunique"
+            << std::endl
+            << "minimizers. The output here has been eyeballed as decent."
+            << std::endl;
+    for (auto cell : triangulation.active_cell_iterators())
+      {
+        const Point<3> projected_point = GridTools::project_to_object(cell,
+                                         trial_point);
+        deallog << projected_point[0]
+                << ", "
+                << projected_point[1]
+                << ", "
+                << projected_point[2]
+                << std::endl;
+
+        Assert((torus_manifold.push_forward(torus_manifold.pull_back(projected_point)) -
+                projected_point).norm() < 1e-14, ExcInternalError());
+      }
+  }
+
+  // refine in 3D a few times so that we can observe that the projection error
+  // drops proportionally as the grid is refined
+  for (unsigned int n_refinements = 4; n_refinements < 7; ++ n_refinements)
+    {
+
+      deallog << "====================================================="
+              << std::endl
+              << "Number of global refinements: " << n_refinements
+              << std::endl
+              << "====================================================="
+              << std::endl;
+
+      SphericalManifold<3> spherical_manifold;
+      Triangulation<3> triangulation;
+      constexpr types::manifold_id spherical_id = 42;
+      GridGenerator::hyper_shell(triangulation, Point<3>(), 1.0, 2.0);
+      triangulation.set_manifold(spherical_id, spherical_manifold);
+      triangulation.set_all_manifold_ids(42);
+
+      triangulation.refine_global(n_refinements);
+
+      auto cell = triangulation.begin_active();
+      while (!cell->at_boundary())
+        ++cell;
+      unsigned int face_n = 0;
+      while (!cell->face(face_n)->at_boundary())
+        ++face_n;
+      const auto face = cell->face(face_n);
+
+      // easy test: project the arithmetic mean of the vertices onto the face
+      {
+        const Point<3> trial_point = 1.1*face->center();
+        const Point<3> projected_point = GridTools::project_to_object(face, trial_point);
+
+        deallog << std::endl
+                << "Project a reweighed center onto a face in 3D:"
+                << std::endl
+                << std::endl
+                << "vertex 0 distance from origin: "
+                << face->vertex(0).distance(Point<3>()) << std::endl
+                << "vertex 1 distance from origin: "
+                << face->vertex(1).distance(Point<3>()) << std::endl
+                << "trial point distance from origin: "
+                << trial_point.distance(Point<3>()) << std::endl
+                << "projected point distance from origin: "
+                << projected_point.distance(Point<3>()) << std::endl;
+
+        const std::array<Point<3>, 4> vertices
+        {{face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
+        const std::array<double, 4> weights {{0.25, 0.25, 0.25, 0.25}};
+        const auto vertices_view = make_array_view(vertices.begin(), vertices.end());
+        const auto weights_view = make_array_view(weights.begin(), weights.end());
+
+        const Point<3> true_midpoint = face->get_manifold().get_new_point(vertices_view,
+                                       weights_view);
+        deallog << "projected point:    "
+                << projected_point << std::endl
+                << "true line midpoint: "
+                << true_midpoint
+                << std::endl
+                << "distance: "
+                << projected_point.distance(true_midpoint)
+                << std::endl;
+      }
+
+      // easy test: project with spacedim == structdim
+      {
+        const Point<3> trial_point {0.417941, 4242424242.4242, 1};
+        const Point<3> projected_point = GridTools::project_to_object(cell,
+                                         trial_point);
+        deallog << std::endl
+                << "Check that projecting with spacedim == structdim "
+                << "is the identity map:"
+                << std::endl
+                << std::endl
+                << "trial point:     "
+                << trial_point
+                << std::endl
+                << "projected point: "
+                << projected_point
+                << std::endl;
+      }
+
+      // project a vertex onto the surface
+      {
+        const Point<3> trial_point = face->vertex(0);
+        const Point<3> projected_point = GridTools::project_to_object(face, trial_point);
+
+        deallog << std::endl
+                << "Project a vertex:"
+                << std::endl
+                << std::endl
+                << "vertex 0 distance from origin:        "
+                << face->vertex(0).distance(Point<3>())
+                << std::endl
+                << "trial point distance from origin:     "
+                << trial_point.distance(Point<3>())
+                << std::endl
+                << "projected point distance from origin: "
+                << projected_point.distance(Point<3>())
+                << std::endl
+                << "projected point: "
+                << projected_point
+                << std::endl
+                << "actual vertex:   "
+                << trial_point
+                << std::endl
+                << "distance:        "
+                << trial_point.distance(projected_point)
+                << std::endl;
+      }
+
+      // project (nearly) a vertex onto the surface
+      {
+        const std::array<Point<3>, 2> vertices {{face->vertex(0), face->vertex(1)}};
+        const std::array<double, 2> weights {{1.0 - 1.0/2048.0, 1.0/2048.0}};
+        const auto vertices_view = make_array_view(vertices.begin(), vertices.end());
+        const auto weights_view = make_array_view(weights.begin(), weights.end());
+
+        const Point<3> trial_point = face->get_manifold().get_new_point(vertices_view,
+                                     weights_view);
+        const Point<3> projected_point = GridTools::project_to_object(face, trial_point);
+
+        deallog << std::endl
+                << "Project a point near a vertex:"
+                << std::endl
+                << std::endl
+                << "vertex 0 distance from origin:        "
+                << face->vertex(0).distance(Point<3>())
+                << std::endl
+                << "trial point distance from origin:     "
+                << trial_point.distance(Point<3>())
+                << std::endl
+                << "projected point distance from origin: "
+                << projected_point.distance(Point<3>())
+                << std::endl
+                << "projected point: "
+                << projected_point
+                << std::endl
+                << "trial point:     "
+                << trial_point
+                << std::endl
+                << "distance:        "
+                << projected_point.distance(trial_point)
+                << std::endl;
+      }
+
+      // test that we can recover a point that is on the face
+      {
+        const std::array<Point<3>, 4> vertices
+        {{face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
+        const std::array<double, 4> weights {{0.0625, 0.5, 0.25, 0.1875}};
+        const auto vertices_view = make_array_view(vertices.begin(), vertices.end());
+        const auto weights_view = make_array_view(weights.begin(), weights.end());
+
+        const Point<3> trial_point = face->get_manifold().get_new_point(vertices_view,
+                                     weights_view);
+        const Point<3> projected_point = GridTools::project_to_object(face, trial_point);
+
+        deallog << std::endl
+                << "Project a weighed face point onto the same face in 3D:"
+                << std::endl
+                << std::endl
+                << "trial point distance from origin:     "
+                << trial_point.distance(Point<3>()) << std::endl
+                << "projected point distance from origin: "
+                << projected_point.distance(Point<3>())
+                << std::endl
+                << "Error:                                "
+                << projected_point.distance(trial_point)
+                << std::endl;
+
+      }
+
+      // test that we can recover a point that is moved off the surface along a
+      // normal vector
+      {
+        const std::array<Point<3>, 4> vertices
+        {{face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
+        const std::array<double, 4> weights {{0.0625, 0.5, 0.25, 0.1875}};
+        const auto vertices_view = make_array_view(vertices.begin(), vertices.end());
+        const auto weights_view = make_array_view(weights.begin(), weights.end());
+
+        Point<3> trial_point = face->get_manifold().get_new_point(vertices_view,
+                                                                  weights_view);
+        Tensor<1, 3> normal_vector = 0.1*(trial_point - Point<3>());
+        trial_point += normal_vector;
+        const Point<3> projected_point = GridTools::project_to_object(face, trial_point);
+
+        deallog << std::endl
+                << "Project a point offset from the face onto the face in 3D:"
+                << std::endl
+                << std::endl
+                << "trial point distance from origin:     "
+                << trial_point.distance(Point<3>()) << std::endl
+                << "projected point distance from origin: "
+                << projected_point.distance(Point<3>())
+                << std::endl
+                << "Error (vs. going along the normal):   "
+                << projected_point.distance(trial_point - normal_vector)
+                << std::endl;
+      }
+
+      // test that we can recover a point that is on a line in 3D
+      {
+        const auto line = face->line(0);
+        const std::array<Point<3>, 2> vertices {{line->vertex(0), line->vertex(1)}};
+        const std::array<double, 2> weights {{0.125, 1.0 - 0.125}};
+        const auto vertices_view = make_array_view(vertices.begin(), vertices.end());
+        const auto weights_view = make_array_view(weights.begin(), weights.end());
+
+        const Point<3> trial_point = face->get_manifold().get_new_point(vertices_view,
+                                     weights_view);
+        const Point<3> projected_point = GridTools::project_to_object(face->line(0), trial_point);
+
+        deallog << std::endl
+                << "Project a weighed line point onto the same line in 3D:"
+                << std::endl
+                << std::endl
+                << "trial point:                          "
+                << trial_point
+                << std::endl
+                << "projected point:                      "
+                << projected_point
+                << std::endl
+                << "trial point distance from origin:     "
+                << trial_point.distance(Point<3>()) << std::endl
+                << "projected point distance from origin: "
+                << projected_point.distance(Point<3>())
+                << std::endl
+                << "Error:                                "
+                << projected_point.distance(trial_point)
+                << std::endl;
+      }
+    }
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/grid/project_to_object_01.output b/tests/grid/project_to_object_01.output
new file mode 100644 (file)
index 0000000..26d19c3
--- /dev/null
@@ -0,0 +1,235 @@
+
+DEAL::testing cross derivative 0
+DEAL::2.08262e-05
+DEAL::2.61032e-06
+DEAL::3.26714e-07
+DEAL::4.08654e-08
+DEAL::5.10984e-09
+DEAL::testing cross derivative 1
+DEAL::2.10357e-05
+DEAL::2.62342e-06
+DEAL::3.27533e-07
+DEAL::4.09167e-08
+DEAL::5.11266e-09
+DEAL::testing cross derivative 2
+DEAL::5.75514e-06
+DEAL::7.19207e-07
+DEAL::8.98876e-08
+DEAL::1.12350e-08
+DEAL::1.40439e-09
+DEAL::
+DEAL::Project the arithmetic mean of two vertices on a circular
+DEAL::arc. The result should be the geodesic midpoint:
+DEAL::
+DEAL::vertex 0 distance from origin: 2.00000
+DEAL::vertex 1 distance from origin: 2.00000
+DEAL::trial point distance from origin: 1.79445
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point:        1.99383 0.156918
+DEAL::true geodesic midpoint: 1.99383 0.156918
+DEAL::
+DEAL::Project the convex combination of two vertices:
+DEAL::projected point:          1.98114 0.274025
+DEAL::Convex combination point: 1.98114 0.274025
+DEAL::Distance:                 5.55112e-17
+DEAL::
+DEAL::Check that projecting with spacedim == structdim
+DEAL::is the identity map:
+DEAL::
+DEAL::trial point:     0.417941 4.24242e+09
+DEAL::projected point: 0.417941 4.24242e+09
+DEAL::
+DEAL::Project a vertex onto the geodesic:
+DEAL::
+DEAL::vertex distance from origin:          2.00000
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point: 2.00000 0.00000
+DEAL::actual vertex:   2.00000 0.00000
+DEAL::Test for robustness by projecting points with nonunique
+DEAL::minimizers. The output here has been eyeballed as decent.
+DEAL::0.768956, 1.00000, 1.84627
+DEAL::1.00000, 1.22465e-16, 0.00000
+DEAL::3.00000, 0.00000, 0.00000
+DEAL::2.00000, 1.00000, 0.00000
+DEAL::-1.84627, 1.00000, 0.768956
+DEAL::6.12323e-17, 1.22465e-16, 1.00000
+DEAL::1.83697e-16, 0.00000, 3.00000
+DEAL::1.22465e-16, 1.00000, 2.00000
+DEAL::-2.00000, 1.00000, 2.44929e-16
+DEAL::-1.00000, 1.22465e-16, 1.22465e-16
+DEAL::-3.00000, 0.00000, 3.67394e-16
+DEAL::-2.00000, 1.00000, 2.44929e-16
+DEAL::1.84627, 1.00000, -0.768956
+DEAL::-1.83697e-16, 1.22465e-16, -1.00000
+DEAL::-5.51091e-16, 0.00000, -3.00000
+DEAL::-3.67394e-16, 1.00000, -2.00000
+DEAL::=====================================================
+DEAL::Number of global refinements: 4
+DEAL::=====================================================
+DEAL::
+DEAL::Project a reweighed center onto a face in 3D:
+DEAL::
+DEAL::vertex 0 distance from origin: 2.00000
+DEAL::vertex 1 distance from origin: 2.00000
+DEAL::trial point distance from origin: 2.19658
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point:    -1.12085 -1.12085 -1.21959
+DEAL::true line midpoint: -1.12085 -1.12085 -1.21959
+DEAL::distance: 8.01926e-07
+DEAL::
+DEAL::Check that projecting with spacedim == structdim is the identity map:
+DEAL::
+DEAL::trial point:     0.417941 4.24242e+09 1.00000
+DEAL::projected point: 0.417941 4.24242e+09 1.00000
+DEAL::
+DEAL::Project a vertex:
+DEAL::
+DEAL::vertex 0 distance from origin:        2.00000
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point: -1.15470 -1.15470 -1.15470
+DEAL::actual vertex:   -1.15470 -1.15470 -1.15470
+DEAL::distance:        0.00000
+DEAL::
+DEAL::Project a point near a vertex:
+DEAL::
+DEAL::vertex 0 distance from origin:        2.00000
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point: -1.15464 -1.15473 -1.15473
+DEAL::trial point:     -1.15464 -1.15473 -1.15473
+DEAL::distance:        3.63867e-13
+DEAL::
+DEAL::Project a weighed face point onto the same face in 3D:
+DEAL::
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error:                                4.74583e-11
+DEAL::
+DEAL::Project a point offset from the face onto the face in 3D:
+DEAL::
+DEAL::trial point distance from origin:     2.20000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error (vs. going along the normal):   3.01973e-06
+DEAL::
+DEAL::Project a weighed line point onto the same line in 3D:
+DEAL::
+DEAL::trial point:                          -1.20701 -1.04223 -1.20701
+DEAL::projected point:                      -1.20701 -1.04223 -1.20701
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error:                                3.28543e-10
+DEAL::=====================================================
+DEAL::Number of global refinements: 5
+DEAL::=====================================================
+DEAL::
+DEAL::Project a reweighed center onto a face in 3D:
+DEAL::
+DEAL::vertex 0 distance from origin: 2.00000
+DEAL::vertex 1 distance from origin: 2.00000
+DEAL::trial point distance from origin: 2.19917
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point:    -1.13839 -1.13839 -1.18664
+DEAL::true line midpoint: -1.13839 -1.13839 -1.18665
+DEAL::distance: 2.27995e-06
+DEAL::
+DEAL::Check that projecting with spacedim == structdim is the identity map:
+DEAL::
+DEAL::trial point:     0.417941 4.24242e+09 1.00000
+DEAL::projected point: 0.417941 4.24242e+09 1.00000
+DEAL::
+DEAL::Project a vertex:
+DEAL::
+DEAL::vertex 0 distance from origin:        2.00000
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point: -1.15470 -1.15470 -1.15470
+DEAL::actual vertex:   -1.15470 -1.15470 -1.15470
+DEAL::distance:        0.00000
+DEAL::
+DEAL::Project a point near a vertex:
+DEAL::
+DEAL::vertex 0 distance from origin:        2.00000
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point: -1.15467 -1.15472 -1.15472
+DEAL::trial point:     -1.15467 -1.15472 -1.15472
+DEAL::distance:        4.50346e-11
+DEAL::
+DEAL::Project a weighed face point onto the same face in 3D:
+DEAL::
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error:                                3.98773e-06
+DEAL::
+DEAL::Project a point offset from the face onto the face in 3D:
+DEAL::
+DEAL::trial point distance from origin:     2.20000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error (vs. going along the normal):   2.12114e-05
+DEAL::
+DEAL::Project a weighed line point onto the same line in 3D:
+DEAL::
+DEAL::trial point:                          -1.18152 -1.09909 -1.18152
+DEAL::projected point:                      -1.18152 -1.09909 -1.18152
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error:                                2.35082e-08
+DEAL::=====================================================
+DEAL::Number of global refinements: 6
+DEAL::=====================================================
+DEAL::
+DEAL::Project a reweighed center onto a face in 3D:
+DEAL::
+DEAL::vertex 0 distance from origin: 2.00000
+DEAL::vertex 1 distance from origin: 2.00000
+DEAL::trial point distance from origin: 2.19979
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point:    -1.14670 -1.14670 -1.17054
+DEAL::true line midpoint: -1.14670 -1.14670 -1.17054
+DEAL::distance: 1.03418e-07
+DEAL::
+DEAL::Check that projecting with spacedim == structdim is the identity map:
+DEAL::
+DEAL::trial point:     0.417941 4.24242e+09 1.00000
+DEAL::projected point: 0.417941 4.24242e+09 1.00000
+DEAL::
+DEAL::Project a vertex:
+DEAL::
+DEAL::vertex 0 distance from origin:        2.00000
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point: -1.15470 -1.15470 -1.15470
+DEAL::actual vertex:   -1.15470 -1.15470 -1.15470
+DEAL::distance:        0.00000
+DEAL::
+DEAL::Project a point near a vertex:
+DEAL::
+DEAL::vertex 0 distance from origin:        2.00000
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::projected point: -1.15469 -1.15471 -1.15471
+DEAL::trial point:     -1.15469 -1.15471 -1.15471
+DEAL::distance:        3.78080e-10
+DEAL::
+DEAL::Project a weighed face point onto the same face in 3D:
+DEAL::
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error:                                0.000830302
+DEAL::
+DEAL::Project a point offset from the face onto the face in 3D:
+DEAL::
+DEAL::trial point distance from origin:     2.20000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error (vs. going along the normal):   0.00114132
+DEAL::
+DEAL::Project a weighed line point onto the same line in 3D:
+DEAL::
+DEAL::trial point:                          -1.16828 -1.12706 -1.16828
+DEAL::projected point:                      -1.16828 -1.12706 -1.16828
+DEAL::trial point distance from origin:     2.00000
+DEAL::projected point distance from origin: 2.00000
+DEAL::Error:                                1.14617e-07
+DEAL::OK

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.