Add tolerance to avoid empty CGAL intersections due to roundoff

include/deal.II/base/utilities.h

@@ -242,6 +242,23 @@ namespace Utilities
   Number
   truncate_to_n_digits(const Number number, const unsigned int n_digits);
 
+  /**
+   * This function allows to round a floating point number @p number to n_digits
+   * @note This does not imply an exact rounding but it ensures two numbers that
+   * should be the same (if rounded analytically) have the same floating point
+   * representation.
+   */
+  template <typename Number>
+  Number
+  round(const Number number, const unsigned int n_digits);
+
+  /**
+   * Same function as above working directly on given @p Point objects.
+   */
+  template <int dim, typename Number>
+  void
+  round(Point<dim, Number> &point, const unsigned int n_digits);
+
   /**
    * Given a string, convert it to an integer. Throw an assertion if that is
    * not possible.

include/deal.II/cgal/intersections.h

@@ -46,8 +46,19 @@ namespace CGALWrappers
    * @param cell1 Iterator to the second cell.
    * @param mapping0 Mapping for the first cell.
    * @param mapping1 Mapping for the second cell.
-   * @param tol Treshold to decide whether or not a simplex is included.
+   * @param discard_tolerance Treshold to decide whether or not a found simplex is included.
+   * @param vertex_tolerance Tolerance to find an intersection. 
+   * In some cases intersections should only be found if cells are intersecting up to machine 
+   * precision. This is the default case.
+   * In some cases however, we want to find an intersection with a certain tolerance. This is 
+   * for example the case if we constructed a triangulation from two non-fitting 
+   * triangulations to realize required jumps in element sizes. To find intersections up to a 
+   * certain tolerance we slightly modify the vertex positions before asking CGAL about 
+   * intersections. This is done by enforcing the same floating point representation for 
+   * vertices that are equivalent up to a certain tolerance using @c Utilities::round().
    * @return Vector of arrays, where each array identify a simplex by its vertices.
+   * @note Rounding with vertex_tolerance does not have any effect on the triangulation, it is 
+   * only applied to a copy of vertices handed to CGAL.
    */
   template <int dim0, int dim1, int spacedim>
   std::vector<std::array<Point<spacedim>, dim1 + 1>>
@@ -56,7 +67,8 @@ namespace CGALWrappers
     const typename Triangulation<dim1, spacedim>::cell_iterator &cell1,
     const Mapping<dim0, spacedim> &                              mapping0,
     const Mapping<dim1, spacedim> &                              mapping1,
-    const double                                                 tol = 1e-9);
+    const double discard_tolerance = 1e-9,
+    const double vertex_tolerance  = std::numeric_limits<double>::min());
 
 
   /**
@@ -67,11 +79,11 @@ namespace CGALWrappers
   compute_intersection_of_cells(
     const std::array<Point<spacedim>, n_vertices0> &vertices0,
     const std::array<Point<spacedim>, n_vertices1> &vertices1,
-    const double                                    tol = 1e-9)
+    const double                                    discard_tolerance = 1e-9)
   {
     (void)vertices0;
     (void)vertices1;
-    (void)tol;
+    (void)discard_tolerance;
     Assert(false, ExcMessage("No explicit template instantiation available"));
     return {};
   }

include/deal.II/cgal/utilities.h

@@ -606,7 +606,8 @@ namespace CGALWrappers
   std::array<Point<spacedim>, n_vertices>
   get_vertices_in_cgal_order(
     const typename dealii::Triangulation<dim, spacedim>::cell_iterator &cell,
-    const Mapping<dim, spacedim> &                                      mapping)
+    const Mapping<dim, spacedim> &                                      mapping,
+    const double tolerance = std::numeric_limits<double>::min())
   {
     // Elements have to be rectangular or simplices
     Assert(n_vertices == std::pow(2, dim) || n_vertices == dim + 1,
@@ -623,6 +624,10 @@ namespace CGALWrappers
         ReferenceCells::Quadrilateral)
       std::swap(vertices[2], vertices[3]);
 
+    if (tolerance > std::numeric_limits<double>::min())
+      for (auto &v : vertices)
+        Utilities::round(v, -std::log10(tolerance));
+
     return vertices;
   }
 

source/base/utilities.cc

@@ -604,6 +604,27 @@ namespace Utilities
   }
 
 
+
+  template <typename Number>
+  Number
+  round(const Number number, const unsigned int n_digits)
+  {
+    const Number eps =
+      0.5 * std::pow(static_cast<Number>(10.0), static_cast<int>(-n_digits));
+    return truncate_to_n_digits(number + eps, n_digits);
+  }
+
+
+  template <int dim, typename Number>
+  void
+  round(Point<dim, Number> &point, const unsigned int n_digits)
+  {
+    for (unsigned int d = 0; d < dim; ++d)
+      point[d] = round(point[d], n_digits);
+  }
+
+
+
   int
   string_to_int(const std::string &s_)
   {
@@ -1094,6 +1115,24 @@ namespace Utilities
   template float
   truncate_to_n_digits(const float, const unsigned int);
 
+  template double
+  round(const double, const unsigned int);
+  template float
+  round(const float, const unsigned int);
+
+  template void
+  round(Point<1, double> &, const unsigned int);
+  template void
+  round(Point<2, double> &, const unsigned int);
+  template void
+  round(Point<3, double> &, const unsigned int);
+  template void
+  round(Point<1, float> &, const unsigned int);
+  template void
+  round(Point<2, float> &, const unsigned int);
+  template void
+  round(Point<3, float> &, const unsigned int);
+
   template std::vector<std::array<std::uint64_t, 1>>
   inverse_Hilbert_space_filling_curve<1, double>(
     const std::vector<Point<1, double>> &,

source/cgal/intersections.cc

@@ -353,7 +353,7 @@ namespace CGALWrappers
   compute_intersection_of_cells<2, 2, 2, 4, 4>(
     const std::array<Point<2>, 4> &vertices0,
     const std::array<Point<2>, 4> &vertices1,
-    const double                   tol)
+    const double                   discard_tolerance)
   {
     const auto intersection_test =
       internal::compute_intersection(vertices0, vertices1);
@@ -387,7 +387,7 @@ namespace CGALWrappers
             for (const Face_handle &f : cdt.finite_face_handles())
               {
                 if (f->info().in_domain() &&
-                    CGAL::to_double(cdt.triangle(f).area()) > tol)
+                    CGAL::to_double(cdt.triangle(f).area()) > discard_tolerance)
                   {
                     collection.push_back(
                       {{CGALWrappers::cgal_point_to_dealii_point<2>(
@@ -420,7 +420,7 @@ namespace CGALWrappers
   compute_intersection_of_cells<2, 1, 2, 4, 2>(
     const std::array<Point<2>, 4> &vertices0,
     const std::array<Point<2>, 2> &vertices1,
-    const double                   tol)
+    const double                   discard_tolerance)
   {
     std::array<CGALPoint2, 4> pts;
     std::transform(
@@ -440,7 +440,7 @@ namespace CGALWrappers
     for (Face_handle f : cdt.finite_face_handles())
       {
         if (f->info().in_domain() &&
-            CGAL::to_double(cdt.triangle(f).area()) > tol &&
+            CGAL::to_double(cdt.triangle(f).area()) > discard_tolerance &&
             CGAL::do_intersect(segm, cdt.triangle(f)))
           {
             const auto intersection = CGAL::intersection(segm, cdt.triangle(f));
@@ -462,7 +462,7 @@ namespace CGALWrappers
   compute_intersection_of_cells<3, 1, 3, 8, 2>(
     const std::array<Point<3>, 8> &vertices0,
     const std::array<Point<3>, 2> &vertices1,
-    const double                   tol)
+    const double                   discard_tolerance)
   {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 1, 5)
     std::array<CGALPoint3_exact, 8> pts;
@@ -489,7 +489,7 @@ namespace CGALWrappers
             if (const CGALSegment3_exact *s =
                   boost::get<CGALSegment3_exact>(&*intersection))
               {
-                if (s->squared_length() > tol * tol)
+                if (s->squared_length() > discard_tolerance * discard_tolerance)
                   {
                     vertices.push_back(
                       {{CGALWrappers::cgal_point_to_dealii_point<3>(
@@ -509,7 +509,7 @@ namespace CGALWrappers
         "This function requires a version of CGAL greater or equal than 5.1.5."));
     (void)vertices0;
     (void)vertices1;
-    (void)tol;
+    (void)discard_tolerance;
     return {};
 #  endif
   }
@@ -519,7 +519,7 @@ namespace CGALWrappers
   compute_intersection_of_cells<3, 2, 3, 8, 4>(
     const std::array<Point<3>, 8> &vertices0,
     const std::array<Point<3>, 4> &vertices1,
-    const double                   tol)
+    const double                   discard_tolerance)
   {
 #  if DEAL_II_CGAL_VERSION_GTE(5, 1, 5)
     std::array<CGALPoint3_exact, 8> pts_hex;
@@ -563,7 +563,8 @@ namespace CGALWrappers
                 if (const CGALTriangle3_exact *t =
                       boost::get<CGALTriangle3_exact>(&*intersection))
                   {
-                    if (CGAL::to_double(t->squared_area()) > tol * tol)
+                    if (CGAL::to_double(t->squared_area()) >
+                        discard_tolerance * discard_tolerance)
                       {
                         vertices.push_back(
                           {{cgal_point_to_dealii_point<3>((*t)[0]),
@@ -584,7 +585,7 @@ namespace CGALWrappers
                       {
                         const auto triangle = tria_inter.triangle(*it);
                         if (CGAL::to_double(triangle.squared_area()) >
-                            tol * tol)
+                            discard_tolerance * discard_tolerance)
                           {
                             std::array<Point<3>, 3> verts = {
                               {CGALWrappers::cgal_point_to_dealii_point<3>(
@@ -610,7 +611,7 @@ namespace CGALWrappers
         "This function requires a version of CGAL greater or equal than 5.1.5."));
     (void)vertices0;
     (void)vertices1;
-    (void)tol;
+    (void)discard_tolerance;
     return {};
 #  endif
   }
@@ -622,8 +623,7 @@ namespace CGALWrappers
   compute_intersection_of_cells<3, 3, 3, 8, 8>(
     const std::array<Point<3>, 8> &vertices0,
     const std::array<Point<3>, 8> &vertices1,
-
-    const double tol)
+    const double                   discard_tolerance)
   {
     std::array<CGALPoint3_inexact, 8> pts_hex0;
     std::array<CGALPoint3_inexact, 8> pts_hex1;
@@ -673,7 +673,7 @@ namespace CGALWrappers
             namespace PMP = CGAL::Polygon_mesh_processing;
             const bool test_intersection =
               PMP::corefine_and_compute_intersection(surf0, surf1, sm);
-            if (PMP::volume(sm) > tol && test_intersection)
+            if (PMP::volume(sm) > discard_tolerance && test_intersection)
               {
                 // Collect tetrahedrons
                 Triangulation3_inexact tria;
@@ -709,7 +709,8 @@ namespace CGALWrappers
     const typename Triangulation<dim1, spacedim>::cell_iterator &cell1,
     const Mapping<dim0, spacedim> &                              mapping0,
     const Mapping<dim1, spacedim> &                              mapping1,
-    const double                                                 tol)
+    const double discard_tolerance,
+    const double vertex_tolerance)
   {
     Assert(mapping0.get_vertices(cell0).size() == std::pow(2, dim0),
            ExcNotImplemented());
@@ -718,18 +719,17 @@ namespace CGALWrappers
 
     const auto vertices0 =
       CGALWrappers::get_vertices_in_cgal_order<Utilities::pow(2, dim0)>(
-        cell0, mapping0);
+        cell0, mapping0, vertex_tolerance);
     const auto vertices1 =
       CGALWrappers::get_vertices_in_cgal_order<Utilities::pow(2, dim1)>(
-        cell1, mapping1);
+        cell1, mapping1, vertex_tolerance);
 
     return compute_intersection_of_cells<dim0,
                                          dim1,
                                          spacedim,
                                          Utilities::pow(2, dim0),
-                                         Utilities::pow(2, dim1)>(vertices0,
-                                                                  vertices1,
-                                                                  tol);
+                                         Utilities::pow(2, dim1)>(
+      vertices0, vertices1, discard_tolerance);
   }
 
 #  include "intersections.inst"
@@ -751,11 +751,11 @@ std::vector<std::array<Point<spacedim>, dim1 + 1>>
 compute_intersection_of_cells(
   const std::array<Point<spacedim>, n_components0> &vertices0,
   const std::array<Point<spacedim>, n_components1> &vertices1,
-  const double                                      tol)
+  const double                                      discard_tolerance)
 {
   (void)vertices0;
   (void)vertices1;
-  (void)tol;
+  (void)discard_tolerance;
   AssertThrow(false, ExcNeedsCGAL());
 }
 
@@ -766,13 +766,15 @@ compute_intersection_of_cells(
   const typename Triangulation<dim1, spacedim>::cell_iterator &cell1,
   const Mapping<dim0, spacedim> &                              mapping0,
   const Mapping<dim1, spacedim> &                              mapping1,
-  const double                                                 tol)
+  const double discard_tolerance,
+  const double vertex_tolerance)
 {
   (void)cell0;
   (void)cell1;
   (void)mapping0;
   (void)mapping1;
-  (void)tol;
+  (void)discard_tolerance;
+  (void)vertex_tolerance;
   AssertThrow(false, ExcNeedsCGAL());
 }
 

source/cgal/intersections.inst.in

@@ -25,7 +25,8 @@ for (dim0 : DIMENSIONS; dim1 : DIMENSIONS; spacedim : SPACE_DIMENSIONS)
       const typename Triangulation<dim1, spacedim>::cell_iterator &cell1,
       const Mapping<dim0, spacedim> &                              mapping0,
       const Mapping<dim1, spacedim> &                              mapping1,
-      const double                                                 tol);
+      const double                                                 tol,
+      const double discard_tolerance);
 
 #endif
   }

tests/base/utilities_23.cc

@@ -0,0 +1,59 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// test functions in namespace Utilities
+
+#include <deal.II/base/utilities.h>
+
+#include "../tests.h"
+
+template <typename Number>
+void
+test()
+{
+  Number number = 1.23456789;
+  assert(Utilities::round(number, 1) == Utilities::round(Number{1.2}, 1));
+  deallog << "1 digit OK" << std::endl;
+  assert(Utilities::round(number, 2) == Utilities::round(Number{1.23}, 2));
+  deallog << "2 digit OK" << std::endl;
+  assert(Utilities::round(number, 3) == Utilities::round(Number{1.235}, 3));
+  deallog << "3 digit OK" << std::endl;
+  assert(Utilities::round(number, 4) == Utilities::round(Number{1.2346}, 4));
+  deallog << "4 digit OK" << std::endl;
+  assert(Utilities::round(number, 5) == Utilities::round(Number{1.23457}, 5));
+  deallog << "5 digit OK" << std::endl;
+  assert(Utilities::round(number, 6) == Utilities::round(Number{1.234568}, 6));
+  deallog << "6 digit OK" << std::endl;
+  assert(Utilities::round(number, 7) == Utilities::round(Number{1.2345679}, 7));
+  deallog << "7 digit OK" << std::endl;
+  assert(Utilities::round(number, 8) ==
+         Utilities::round(Number{1.23456789}, 8));
+  deallog << "8 digit OK" << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  deallog << "Double" << std::endl;
+  test<double>();
+  deallog << "Float" << std::endl;
+  test<float>();
+  return 0;
+}