From: David Wells Date: Mon, 15 May 2023 16:36:25 +0000 (-0400) Subject: Add a function for projecting points onto reference cells. X-Git-Tag: v9.5.0-rc1~161^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0a82faa4f005e528fb9350cf52ee2e63bf511e95;p=dealii.git Add a function for projecting points onto reference cells. This can replace GeometryInfo::distance_to_unit_cell() and GeometryInfo::project_to_unit_cell() in most circumstances. There is no clear way to define the infinity norm for sloped surfaces (like the top face of a tetrahedron), so this new function uses the Euclidean norm. --- diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index db61250178..dccf46541d 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -588,7 +588,15 @@ public: bool contains_point(const Point &p, const double tolerance = 0) const; - /* + /** + * Return the point on the surface of the reference cell closest (in the + * Euclidean norm) to @p p. + */ + template + Point + closest_point(const Point &p) const; + + /** * Return $i$-th unit tangential vector of a face of the reference cell. * The vectors are arranged such that the * cross product between the two vectors returns the unit normal vector. diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc index 2eb2aa2d51..c9ce0979df 100644 --- a/source/grid/reference_cell.cc +++ b/source/grid/reference_cell.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -807,6 +808,329 @@ ReferenceCell::gmsh_element_type() const +namespace +{ + // Compute the nearest point to @p on the line segment and the square of its + // distance to @p. + template + std::pair, double> + project_to_line(const Point &x0, + const Point &x1, + const Point &p) + { + Assert(x0 != x1, ExcInternalError()); + // t is the convex combination coefficient (x = (1 - t) * x0 + t * x1) + // defining the position of the closest point on the line (not line segment) + // to p passing through x0 and x1. This formula is equivalent to the + // standard 'project a vector onto another vector', where each vector is + // shifted to start at x0. + const double t = ((x1 - x0) * (p - x0)) / ((x1 - x0).norm_square()); + + // Only consider points between x0 and x1 + if (t <= 0) + return std::make_pair(x0, x0.distance_square(p)); + else if (t <= 1) + { + const auto p2 = x0 + t * (x1 - x0); + return std::make_pair(p2, p2.distance_square(p)); + } + else + return std::make_pair(x1, x1.distance_square(p)); + } + + // template base-case + template + std::pair, double> + project_to_quad(const std::array, 3> & /*vertices*/, + const Point & /*p*/, + const ReferenceCell /*reference_cell*/) + { + Assert(false, ExcInternalError()); + return std::make_pair(Point(), + std::numeric_limits::signaling_NaN()); + } + + /** + * Compute the nearest point on a quad (in the deal.II sense: i.e., something + * with structdim = 2 and spacedim = 3) and the square of that point's + * distance to @p p. Here, the quad is described with three vertices: either + * the three vertices of a Triangle or the first three of a Quadrilateral (as + * the fourth one can be computed, in that case, from the first three). + * + * If the given point cannot be projected via a normal vector (i.e., if the + * line parallel to the normal vector intersecting @p p does not intersect the + * quad) then this function returns the origin and the largest double + * precision number. distance_to_line_square() is in charge of computing the + * distance to lines. + * + * @note This function is for Quadrilaterals and Triangles because it is + * intended to find the shortest distance to the face of a reference cell + * (which must be a Triangle or Quadrilateral) in 3d. + */ + template <> + std::pair, double> + project_to_quad(const std::array, 3> &vertices, + const Point<3> & p, + const ReferenceCell face_reference_cell) + { + Assert(face_reference_cell == ReferenceCells::Triangle || + face_reference_cell == ReferenceCells::Quadrilateral, + ExcNotImplemented()); + + // Make the problem slightly easier by shifting everything to avoid a point + // at the origin (this way we can invert the matrix of vertices). Use 2.0 so + // that the bottom left vertex of a Pyramid is now at x = 1. + std::array, 3> shifted_vertices = vertices; + const Tensor<1, 3> shift{{2.0, 2.0, 2.0}}; + for (Point<3> &shifted_vertex : shifted_vertices) + shifted_vertex += shift; + const Point<3> shifted_p = p + shift; + + // As we are projecting onto a face of a reference cell, the vectors + // describing its local coordinate system should be orthogonal. We don't + // know which of the three vectors computed from p are mutually orthogonal + // for triangles so that case requires an extra check. + Tensor<1, 3> e0; + Tensor<1, 3> e1; + const Point<3> vertex = shifted_vertices[0]; + // Triangles are difficult because of two cases: + // 1. the top face of a Tetrahedron, which does not have a right angle + // 2. wedges and pyramids, whose faces do not lie on the reference simplex + // + // Deal with both by creating a locally orthogonal (but not necessarily + // orthonormal) coordinate system and testing if the projected point is in + // the triangle by expressing it as a convex combination of the vertices. + if (face_reference_cell == ReferenceCells::Triangle) + { + e0 = shifted_vertices[1] - shifted_vertices[0]; + e1 = shifted_vertices[2] - shifted_vertices[0]; + e1 -= (e0 * e1) * e0 / (e0.norm_square()); + } + else + { + e0 = shifted_vertices[1] - shifted_vertices[0]; + e1 = shifted_vertices[2] - shifted_vertices[0]; + } + Assert(std::abs(e0 * e1) <= 1e-14, ExcInternalError()); + // the quadrilaterals on pyramids and wedges don't necessarily have edge + // lengths of 1 so we cannot skip the denominator + const double c0 = e0 * (shifted_p - vertex) / e0.norm_square(); + const double c1 = e1 * (shifted_p - vertex) / e1.norm_square(); + const Point<3> projected_shifted_p = vertex + c0 * e0 + c1 * e1; + + bool in_quad = false; + if (face_reference_cell == ReferenceCells::Triangle) + { + Tensor<2, 3> shifted_vertex_matrix; + for (unsigned int i = 0; i < 3; ++i) + shifted_vertex_matrix[i] = shifted_vertices[i]; + const Tensor<1, 3> combination_coordinates = + invert(transpose(shifted_vertex_matrix)) * projected_shifted_p; + bool is_convex_combination = true; + for (unsigned int i = 0; i < 3; ++i) + is_convex_combination = is_convex_combination && + (0.0 <= combination_coordinates[i]) && + (combination_coordinates[i] <= 1.0); + in_quad = is_convex_combination; + } + else + in_quad = (0.0 <= c0 && c0 <= 1.0 && 0.0 <= c1 && c1 <= 1.0); + + if (in_quad) + return std::make_pair(projected_shifted_p - shift, + shifted_p.distance_square(projected_shifted_p)); + else + return std::make_pair(Point<3>(), std::numeric_limits::max()); + } +} // namespace + + + +template +Point +ReferenceCell::closest_point(const Point &p) const +{ + AssertDimension(dim, get_dimension()); + + // Handle simple cases first: + if (dim == 0) + return p; + if (contains_point(p, 0.0)) + return p; + if (dim == 1) + return project_to_line(vertex(0), vertex(1), p).first; + + // Find the closest vertex so that we only need to check adjacent faces and + // lines. + Point result; + unsigned int closest_vertex_no = 0; + double closest_vertex_distance_square = vertex(0).distance_square(p); + for (unsigned int i = 1; i < n_vertices(); ++i) + { + const double new_vertex_distance_square = + vertex(i).distance_square(p); + if (new_vertex_distance_square < closest_vertex_distance_square) + { + closest_vertex_no = i; + closest_vertex_distance_square = new_vertex_distance_square; + } + } + + double min_distance_square = std::numeric_limits::max(); + if (dim == 2) + { + for (const unsigned int face_no : + faces_for_given_vertex(closest_vertex_no)) + { + const Point v0 = vertex(line_to_cell_vertices(face_no, 0)); + const Point v1 = vertex(line_to_cell_vertices(face_no, 1)); + + auto pair = project_to_line(v0, v1, p); + if (pair.second < min_distance_square) + { + result = pair.first; + min_distance_square = pair.second; + } + } + } + else + { + // Check faces and then lines. + // + // For reference cells with sloped faces (i.e., all 3D shapes except + // Hexahedra), we might be able to do a valid normal projection to a face + // with a different slope which is on the 'other side' of the reference + // cell. To catch that case we have to unconditionally check lines after + // checking faces. + // + // For pyramids the closest vertex might not be on the closest face: for + // example, the origin is closest to vertex 4 which is not on the bottom + // plane. Get around that by just checking all faces for pyramids. + const std::array all_pyramid_faces{{0, 1, 2, 3, 4}}; + const auto &faces = *this == ReferenceCells::Pyramid ? + ArrayView(all_pyramid_faces) : + faces_for_given_vertex(closest_vertex_no); + for (const unsigned int face_no : faces) + { + auto face_cell = face_reference_cell(face_no); + // We only need the first three points since for quads the last point + // is redundant + std::array, 3> vertices; + for (unsigned int vertex_no = 0; vertex_no < 3; ++vertex_no) + vertices[vertex_no] = vertex(face_to_cell_vertices( + face_no, vertex_no, default_combined_face_orientation())); + + auto pair = project_to_quad(vertices, p, face_cell); + if (pair.second < min_distance_square) + { + result = pair.first; + min_distance_square = pair.second; + } + } + + for (const unsigned int face_no : + faces_for_given_vertex(closest_vertex_no)) + { + auto face_cell = face_reference_cell(face_no); + for (const unsigned int face_line_no : face_cell.line_indices()) + { + const auto cell_line_no = + face_to_cell_lines(face_no, + face_line_no, + default_combined_face_orientation()); + const auto v0 = + vertex(line_to_cell_vertices(cell_line_no, 0)); + const auto v1 = + vertex(line_to_cell_vertices(cell_line_no, 1)); + auto pair = project_to_line(v0, v1, p); + if (pair.second < min_distance_square) + { + result = pair.first; + min_distance_square = pair.second; + } + } + } + } + + Assert(min_distance_square < std::numeric_limits::max(), + ExcInternalError()); + + // If necessary, slightly adjust the computed point so that it is closer to + // being on the surface of the reference cell. Due to roundoff it is difficult + // to place points on sloped surfaces (e.g., for Pyramids) so this check isn't + // perfect but does improve the accuracy of the projected point. + if (!contains_point(result, 0.0)) + { + constexpr unsigned int x_index = 0; + constexpr unsigned int y_index = (dim >= 2 ? 1 : 0); + constexpr unsigned int z_index = (dim >= 3 ? 2 : 0); + switch (this->kind) + { + case ReferenceCells::Vertex: + Assert(false, ExcInternalError()); + break; + // the bounds for each dimension of a hypercube are mutually + // independent: + case ReferenceCells::Line: + case ReferenceCells::Quadrilateral: + case ReferenceCells::Hexahedron: + for (unsigned int d = 0; d < dim; ++d) + result[d] = std_cxx17::clamp(result[d], 0.0, 1.0); + // simplices can use the standard definition of a simplex: + break; + case ReferenceCells::Triangle: + result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0); + result[y_index] = + std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]); + break; + case ReferenceCells::Tetrahedron: + result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0); + result[y_index] = + std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]); + result[z_index] = + std_cxx17::clamp(result[z_index], + 0.0, + 1.0 - result[x_index] - result[y_index]); + break; + // wedges and pyramids are more ad-hoc: + case ReferenceCells::Wedge: + result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0); + result[y_index] = + std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]); + result[z_index] = std_cxx17::clamp(result[z_index], 0.0, 1.0); + break; + case ReferenceCells::Pyramid: + { + result[x_index] = std_cxx17::clamp(result[x_index], -1.0, 1.0); + result[y_index] = std_cxx17::clamp(result[y_index], -1.0, 1.0); + // It suffices to transform everything to the first quadrant to + // adjust z: + const auto x_abs = std::abs(result[x_index]); + const auto y_abs = std::abs(result[y_index]); + + if (y_abs <= x_abs) + result[z_index] = + std_cxx17::clamp(result[z_index], 0.0, 1.0 - x_abs); + else + result[z_index] = + std_cxx17::clamp(result[z_index], 0.0, 1.0 - y_abs); + } + break; + default: + Assert(false, ExcNotImplemented()); + } + } + // We should be within 4 * eps of the cell by this point. The roundoff error + // comes from, e.g., computing (1 - x) + x when moving points onto the top of + // a Pyramid. + Assert(contains_point(result, 4.0 * std::numeric_limits::epsilon()), + ExcInternalError()); + + return result; +} + + + std::ostream & operator<<(std::ostream &out, const ReferenceCell &reference_cell) { diff --git a/source/grid/reference_cell.inst.in b/source/grid/reference_cell.inst.in index 9493dc32e6..6a23a378d1 100644 --- a/source/grid/reference_cell.inst.in +++ b/source/grid/reference_cell.inst.in @@ -29,6 +29,10 @@ for (deal_II_dimension : DIMENSIONS) { template Quadrature ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1D) const; + template const Quadrature &ReferenceCell::get_nodal_type_quadrature() const; + + template Point ReferenceCell::closest_point( + const Point &) const; } diff --git a/tests/grid/closest_point.cc b/tests/grid/closest_point.cc new file mode 100644 index 0000000000..51cd742287 --- /dev/null +++ b/tests/grid/closest_point.cc @@ -0,0 +1,122 @@ +#include + +#include "../tests.h" + +template +void +test(const ReferenceCell &reference_cell, + const Point & p, + const bool print = true) +{ + const Point projected = reference_cell.closest_point(p); + const double distance_square = projected.distance_square(p); + // its difficult to place points exactly on sloped surfaces (e.g., for + // pyramids), but we should always be close + AssertThrow(reference_cell.contains_point( + projected, 2.0 * std::numeric_limits::epsilon()), + ExcInternalError()); + if (print) + deallog << "point = " << p << std::endl + << " d^2 = " << distance_square << std::endl; + + // The distance we get should be at least as good as the distance to the + // nearest vertex + for (unsigned int vertex_no : reference_cell.vertex_indices()) + AssertThrow(distance_square <= + p.distance_square( + reference_cell.template vertex(vertex_no)), + ExcInternalError()); + + // Also generate a lot of points in the unit cell - we should be closer than + // any generated point + const unsigned int n_points_1d = 100; + if (dim == 2) + { + for (unsigned int i = 0; i < n_points_1d; ++i) + for (unsigned int j = 0; j < n_points_1d; ++j) + { + const Point grid_point(i / double(n_points_1d), + j / double(n_points_1d)); + if (reference_cell.contains_point(grid_point)) + AssertThrow(distance_square <= grid_point.distance_square(p), + ExcInternalError()); + } + } + else if (dim == 3) + { + for (unsigned int i = 0; i < n_points_1d; ++i) + for (unsigned int j = 0; j < n_points_1d; ++j) + for (unsigned int k = 0; j < n_points_1d; ++j) + { + // scale things so that we cover the bounding box of both pyramids + // and wedges + const Point grid_point(-1.0 + 2.0 * i / double(n_points_1d), + -1.0 + 2.0 * j / double(n_points_1d), + k / double(n_points_1d)); + if (reference_cell.contains_point(grid_point)) + AssertThrow(distance_square <= grid_point.distance_square(p), + ExcInternalError()); + } + } +} + +int +main() +{ + initlog(); + + for (const auto &reference_cell : {ReferenceCells::Line}) + { + deallog << "reference cell = " << reference_cell.to_string() << std::endl; + test(reference_cell, Point<1>(-1)); + test(reference_cell, Point<1>(-0.5)); + test(reference_cell, Point<1>(0.5)); + test(reference_cell, Point<1>(1.0)); + test(reference_cell, Point<1>(1.25)); + } + + for (const auto &reference_cell : + {ReferenceCells::Triangle, ReferenceCells::Quadrilateral}) + { + deallog << "reference cell = " << reference_cell.to_string() << std::endl; + test(reference_cell, Point<2>(-0.5, 0.5)); + test(reference_cell, Point<2>(-1, -1)); + test(reference_cell, Point<2>(0.5, -0.5)); + test(reference_cell, Point<2>(1.0, 0.0)); + test(reference_cell, Point<2>(1.0, 1.0)); + test(reference_cell, Point<2>(0.25, 0.25)); + + for (unsigned int i = 0; i < 1e3; ++i) + test(reference_cell, + Point<2>(random_value(-4.0, 4.0), random_value(-4.0, 4.0)), + false); + } + + for (const auto &reference_cell : {ReferenceCells::Tetrahedron, + ReferenceCells::Pyramid, + ReferenceCells::Wedge, + ReferenceCells::Hexahedron}) + { + deallog << "reference cell = " << reference_cell.to_string() << std::endl; + test(reference_cell, Point<3>(-0.5, 0.5, 0.5)); + test(reference_cell, Point<3>(-1, -1, -1)); + test(reference_cell, Point<3>(0.5, -0.5, -0.5)); + test(reference_cell, Point<3>(1.0, 0.0, 0.0)); + test(reference_cell, Point<3>(1.0, 1.0, 1.0)); + test(reference_cell, Point<3>(-0.75, 0.5, 1.5)); + test(reference_cell, Point<3>(2.0, 2.0, 1.0)); + test(reference_cell, Point<3>(0.25, 0.25, 0.25)); + test(reference_cell, Point<3>(1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0)); + test(reference_cell, Point<3>(0.7, 0.9, 0.1)); + test(reference_cell, Point<3>(0.08, 0.18, -0.07)); + test(reference_cell, Point<3>(0.0, 0.0, 0.0)); + + for (unsigned int i = 0; i < 1e3; ++i) + test(reference_cell, + Point<3>(random_value(-4.0, 4.0), + random_value(-4.0, 4.0), + random_value(-4.0, 4.0)), + false); + } + deallog << "OK!" << std::endl; +} diff --git a/tests/grid/closest_point.output b/tests/grid/closest_point.output new file mode 100644 index 0000000000..0f0aeec7e6 --- /dev/null +++ b/tests/grid/closest_point.output @@ -0,0 +1,139 @@ + +DEAL::reference cell = Line +DEAL::point = -1.00000 +DEAL:: d^2 = 1.00000 +DEAL::point = -0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = 0.500000 +DEAL:: d^2 = 0.00000 +DEAL::point = 1.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = 1.25000 +DEAL:: d^2 = 0.0625000 +DEAL::reference cell = Tri +DEAL::point = -0.500000 0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = -1.00000 -1.00000 +DEAL:: d^2 = 2.00000 +DEAL::point = 0.500000 -0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = 1.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = 1.00000 1.00000 +DEAL:: d^2 = 0.500000 +DEAL::point = 0.250000 0.250000 +DEAL:: d^2 = 0.00000 +DEAL::reference cell = Quad +DEAL::point = -0.500000 0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = -1.00000 -1.00000 +DEAL:: d^2 = 2.00000 +DEAL::point = 0.500000 -0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = 1.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = 1.00000 1.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.250000 0.250000 +DEAL:: d^2 = 0.00000 +DEAL::reference cell = Tet +DEAL::point = -0.500000 0.500000 0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = -1.00000 -1.00000 -1.00000 +DEAL:: d^2 = 3.00000 +DEAL::point = 0.500000 -0.500000 -0.500000 +DEAL:: d^2 = 0.500000 +DEAL::point = 1.00000 0.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = 1.00000 1.00000 1.00000 +DEAL:: d^2 = 1.33333 +DEAL::point = -0.750000 0.500000 1.50000 +DEAL:: d^2 = 1.06250 +DEAL::point = 2.00000 2.00000 1.00000 +DEAL:: d^2 = 5.50000 +DEAL::point = 0.250000 0.250000 0.250000 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.333333 0.333333 0.333333 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.700000 0.900000 0.100000 +DEAL:: d^2 = 0.190000 +DEAL::point = 0.0800000 0.180000 -0.0700000 +DEAL:: d^2 = 0.00490000 +DEAL::point = 0.00000 0.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::reference cell = Pyramid +DEAL::point = -0.500000 0.500000 0.500000 +DEAL:: d^2 = 0.00000 +DEAL::point = -1.00000 -1.00000 -1.00000 +DEAL:: d^2 = 1.00000 +DEAL::point = 0.500000 -0.500000 -0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = 1.00000 0.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = 1.00000 1.00000 1.00000 +DEAL:: d^2 = 0.666667 +DEAL::point = -0.750000 0.500000 1.50000 +DEAL:: d^2 = 0.875000 +DEAL::point = 2.00000 2.00000 1.00000 +DEAL:: d^2 = 3.00000 +DEAL::point = 0.250000 0.250000 0.250000 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.333333 0.333333 0.333333 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.700000 0.900000 0.100000 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.0800000 0.180000 -0.0700000 +DEAL:: d^2 = 0.00490000 +DEAL::point = 0.00000 0.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::reference cell = Wedge +DEAL::point = -0.500000 0.500000 0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = -1.00000 -1.00000 -1.00000 +DEAL:: d^2 = 3.00000 +DEAL::point = 0.500000 -0.500000 -0.500000 +DEAL:: d^2 = 0.500000 +DEAL::point = 1.00000 0.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = 1.00000 1.00000 1.00000 +DEAL:: d^2 = 0.500000 +DEAL::point = -0.750000 0.500000 1.50000 +DEAL:: d^2 = 0.812500 +DEAL::point = 2.00000 2.00000 1.00000 +DEAL:: d^2 = 4.50000 +DEAL::point = 0.250000 0.250000 0.250000 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.333333 0.333333 0.333333 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.700000 0.900000 0.100000 +DEAL:: d^2 = 0.180000 +DEAL::point = 0.0800000 0.180000 -0.0700000 +DEAL:: d^2 = 0.00490000 +DEAL::point = 0.00000 0.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::reference cell = Hex +DEAL::point = -0.500000 0.500000 0.500000 +DEAL:: d^2 = 0.250000 +DEAL::point = -1.00000 -1.00000 -1.00000 +DEAL:: d^2 = 3.00000 +DEAL::point = 0.500000 -0.500000 -0.500000 +DEAL:: d^2 = 0.500000 +DEAL::point = 1.00000 0.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = 1.00000 1.00000 1.00000 +DEAL:: d^2 = 0.00000 +DEAL::point = -0.750000 0.500000 1.50000 +DEAL:: d^2 = 0.812500 +DEAL::point = 2.00000 2.00000 1.00000 +DEAL:: d^2 = 2.00000 +DEAL::point = 0.250000 0.250000 0.250000 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.333333 0.333333 0.333333 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.700000 0.900000 0.100000 +DEAL:: d^2 = 0.00000 +DEAL::point = 0.0800000 0.180000 -0.0700000 +DEAL:: d^2 = 0.00490000 +DEAL::point = 0.00000 0.00000 0.00000 +DEAL:: d^2 = 0.00000 +DEAL::OK!