Bugfix: Add tolerance to plane-check

src/t8_geometry/t8_geometry_helpers.c

@@ -658,16 +658,15 @@ t8_triangle_point_inside (const double p_0[3], const double v[3], const double w
 }
 
 int
-t8_plane_point_inside (const double point_on_face[3], const double face_normal[3], const double point[3])
+t8_plane_point_inside (const double point_on_face[3], const double face_normal[3], const double point[3],
+                       const double tolerance)
 {
+  T8_ASSERT (tolerance > 0); /* negative values and zero are not allowed */
   /* Set x = x - p */
   double pof[3] = { point_on_face[0], point_on_face[1], point_on_face[2] };
   t8_axpy (point, pof, -1);
   /* Compute <x-p,n> */
   const double dot_product = t8_dot (pof, face_normal);
-  if (dot_product < 0) {
-    /* The point is on the wrong side of the plane */
-    return 0;
-  }
-  return 1;
+  /* If the dot product is negative, the point is on the wrong side of the plane */
+  return (fabs (dot_product) < tolerance);
 }

src/t8_geometry/t8_geometry_helpers.h

@@ -218,10 +218,12 @@ t8_triangle_point_inside (const double p_0[3], const double v[3], const double w
  * \param[in] point_on_face   A point on the plane
  * \param[in] face_normal     The normal of the face
  * \param[in] point           The point to check
+ * \param[in] tolerance       A double > 0 defining the tolerance
  * \return                    0 if the point is outside, 1 otherwise.
  */
 int
-t8_plane_point_inside (const double point_on_face[3], const double face_normal[3], const double point[3]);
+t8_plane_point_inside (const double point_on_face[3], const double face_normal[3], const double point[3],
+                       const double tolerance);
 
 T8_EXTERN_C_END ();
 

src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx

@@ -232,7 +232,7 @@ t8_geometry_linear::t8_geom_point_batch_inside_element (t8_forest_t forest, t8_l
       double point_on_face[3];
       t8_forest_element_coordinate (forest, ltreeid, element, afacecorner, point_on_face);
       for (int ipoint = 0; ipoint < num_points; ipoint++) {
-        const int is_inside_iface = t8_plane_point_inside (point_on_face, face_normal, &points[ipoint * 3]);
+        const int is_inside_iface = t8_plane_point_inside (point_on_face, face_normal, &points[ipoint * 3], tolerance);
         if (is_inside_iface == 0) {
           /* Point is on the outside of face iface. Update is_inside */
           is_inside[ipoint] = 0;