]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a function for projecting points onto reference cells. 15218/head
authorDavid Wells <drwells@email.unc.edu>
Mon, 15 May 2023 16:36:25 +0000 (12:36 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sat, 3 Jun 2023 15:11:14 +0000 (11:11 -0400)
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.

include/deal.II/grid/reference_cell.h
source/grid/reference_cell.cc
source/grid/reference_cell.inst.in
tests/grid/closest_point.cc [new file with mode: 0644]
tests/grid/closest_point.output [new file with mode: 0644]

index db61250178306ec1e534aefd11a63ae82ac1a6a7..dccf46541dcd01ba93492cb4ea1542adfe7e6830 100644 (file)
@@ -588,7 +588,15 @@ public:
   bool
   contains_point(const Point<dim> &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 <int dim>
+  Point<dim>
+  closest_point(const Point<dim> &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.
index 2eb2aa2d5192aa881a50f94eac7dadfb092c5c8c..c9ce0979dff8d9d8d90d868d91b6700bb218f4f9 100644 (file)
@@ -16,6 +16,7 @@
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/polynomials_barycentric.h>
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx17/algorithm.h>
 #include <deal.II/base/tensor_product_polynomials.h>
 
 #include <deal.II/fe/fe_pyramid_p.h>
@@ -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 <int dim>
+  std::pair<Point<dim>, double>
+  project_to_line(const Point<dim> &x0,
+                  const Point<dim> &x1,
+                  const Point<dim> &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 <int dim>
+  std::pair<Point<dim>, double>
+  project_to_quad(const std::array<Point<dim>, 3> & /*vertices*/,
+                  const Point<dim> & /*p*/,
+                  const ReferenceCell /*reference_cell*/)
+  {
+    Assert(false, ExcInternalError());
+    return std::make_pair(Point<dim>(),
+                          std::numeric_limits<double>::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<Point<3>, double>
+  project_to_quad(const std::array<Point<3>, 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<Point<3>, 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<double>::max());
+  }
+} // namespace
+
+
+
+template <int dim>
+Point<dim>
+ReferenceCell::closest_point(const Point<dim> &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<dim>(0), vertex<dim>(1), p).first;
+
+  // Find the closest vertex so that we only need to check adjacent faces and
+  // lines.
+  Point<dim>   result;
+  unsigned int closest_vertex_no        = 0;
+  double closest_vertex_distance_square = vertex<dim>(0).distance_square(p);
+  for (unsigned int i = 1; i < n_vertices(); ++i)
+    {
+      const double new_vertex_distance_square =
+        vertex<dim>(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<double>::max();
+  if (dim == 2)
+    {
+      for (const unsigned int face_no :
+           faces_for_given_vertex(closest_vertex_no))
+        {
+          const Point<dim> v0 = vertex<dim>(line_to_cell_vertices(face_no, 0));
+          const Point<dim> v1 = vertex<dim>(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<unsigned int, 5> all_pyramid_faces{{0, 1, 2, 3, 4}};
+      const auto &faces = *this == ReferenceCells::Pyramid ?
+                            ArrayView<const unsigned int>(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<Point<dim>, 3> vertices;
+          for (unsigned int vertex_no = 0; vertex_no < 3; ++vertex_no)
+            vertices[vertex_no] = vertex<dim>(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<dim>(line_to_cell_vertices(cell_line_no, 0));
+              const auto v1 =
+                vertex<dim>(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<double>::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<double>::epsilon()),
+         ExcInternalError());
+
+  return result;
+}
+
+
+
 std::ostream &
 operator<<(std::ostream &out, const ReferenceCell &reference_cell)
 {
index 9493dc32e6d4cd30784b52e563264345d9b5fa38..6a23a378d1089b62b4b56cf6f8f727d9fc4d0f8d 100644 (file)
@@ -29,6 +29,10 @@ for (deal_II_dimension : DIMENSIONS)
   {
     template Quadrature<deal_II_dimension>
     ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1D) const;
+
     template const Quadrature<deal_II_dimension>
       &ReferenceCell::get_nodal_type_quadrature() const;
+
+    template Point<deal_II_dimension> ReferenceCell::closest_point(
+      const Point<deal_II_dimension> &) const;
   }
diff --git a/tests/grid/closest_point.cc b/tests/grid/closest_point.cc
new file mode 100644 (file)
index 0000000..51cd742
--- /dev/null
@@ -0,0 +1,122 @@
+#include <deal.II/grid/reference_cell.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test(const ReferenceCell &reference_cell,
+     const Point<dim> &   p,
+     const bool           print = true)
+{
+  const Point<dim> 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<double>::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<dim>(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<dim> 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<dim> 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 (file)
index 0000000..0f0aeec
--- /dev/null
@@ -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!

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.