From e95ba834e57dfc82d5ec09443415c7237cfc389e Mon Sep 17 00:00:00 2001 From: David Knapp Date: Mon, 23 Feb 2026 15:31:33 +0100 Subject: [PATCH] Add tolerance to plane-check --- src/t8_geometry/t8_geometry_helpers.c | 11 +++++------ src/t8_geometry/t8_geometry_helpers.h | 4 +++- .../t8_geometry_linear.cxx | 2 +- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/t8_geometry/t8_geometry_helpers.c b/src/t8_geometry/t8_geometry_helpers.c index 4fc1e59132..c28c11711e 100644 --- a/src/t8_geometry/t8_geometry_helpers.c +++ b/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 */ 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); } diff --git a/src/t8_geometry/t8_geometry_helpers.h b/src/t8_geometry/t8_geometry_helpers.h index caf57a020a..62575135bb 100644 --- a/src/t8_geometry/t8_geometry_helpers.h +++ b/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 (); diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx index 4b2f2dd161..db64d74ea8 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx +++ b/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;