#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>
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
}
}
}
+
+
+
+ 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> ¢er,
+ 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> ¢er,
+ 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
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
/**
- * 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)
{
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)
||
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
// 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;
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)
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);
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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