Feature point inside check for axis aligned geom

src/t8_forest/t8_forest_cxx.cxx

@@ -1228,7 +1228,8 @@ t8_plane_point_inside (const double point_on_face[3], const double face_normal[3
 
 void
 t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element,
-                                      const double *points, int num_points, int *is_inside, const double tolerance)
+                                      const double *points, int num_points, int *is_inside,
+                                      const int geom_is_axis_aligned, const double tolerance)
 {
   const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid);
   t8_eclass_scheme_c *ts = t8_forest_get_eclass_scheme (forest, tree_class);
@@ -1242,7 +1243,23 @@ t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, c
   const t8_geometry_c *geometry = t8_cmesh_get_tree_geometry (cmesh, cgtreeid);
   T8_ASSERT (t8_geom_is_linear (geometry));
 #endif
+  if (geom_is_axis_aligned) {
+    double v_min[3];
+    double v_max[3];
 
+    /*Geometry is fully described by v_min and v_max*/
+    t8_forest_element_coordinate (forest, ltreeid, element, 0, v_min);
+    t8_forest_element_coordinate (forest, ltreeid, element, 1, v_max);
+
+    for (int ipoint = 0; ipoint < num_points; ipoint++) {
+      /* A point is inside if it is inbetween the x/y/z-coordinates of v_min and v_max */
+      is_inside[ipoint]
+        = v_min[0] <= points[ipoint * 3] && points[ipoint * 3] <= v_max[0] &&         /* check x-coordinate*/
+          v_min[1] <= points[ipoint * 3 + 1] && points[ipoint * 3 + 1] <= v_max[1] && /* check y-coordinate*/
+          v_min[2] <= points[ipoint * 3 + 2] && points[ipoint * 3 + 2] <= v_max[2];   /* check z-coordinate*/
+    }
+    return;
+  }
   switch (element_shape) {
   case T8_ECLASS_VERTEX: {
     /* A point is 'inside' a vertex if they have the same coordinates */
@@ -1380,10 +1397,11 @@ t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, c
 
 int
 t8_forest_element_point_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element,
-                                const double point[3], const double tolerance)
+                                const double point[3], const int geom_is_axis_aligned, const double tolerance)
 {
   int is_inside = 0;
-  t8_forest_element_point_batch_inside (forest, ltreeid, element, point, 1, &is_inside, tolerance);
+  t8_forest_element_point_batch_inside (forest, ltreeid, element, point, 1, &is_inside, geom_is_axis_aligned,
+                                        tolerance);
   return is_inside;
 }
 

src/t8_forest/t8_forest_general.h

@@ -752,14 +752,16 @@ t8_forest_iterate (t8_forest_t forest);
  *                             If this value is larger we detect more points.
  *                             If it is zero we probably do not detect points even if they are inside
  *                             due to rounding errors.
+ * \param [in]      geom_is_axis_aligned Flag to tell if the used geometry is an axis-aligned geometry. Simplifies the
+ *                            check for this type of geometry. Can only be used for lines, quads and hexs. 
  * \return          True (non-zero) if \a point lies within \a element, false otherwise.
  *                  The return value is also true if the point lies on the element boundary.
  *                  Thus, this function may return true for different leaf elements, if they
  *                  are neighbors and the point lies on the common boundary.
  */
 int
 t8_forest_element_point_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element,
-                                const double point[3], const double tolerance);
+                                const double point[3], const int geom_is_axis_aligned, const double tolerance);
 
 /** Query whether a batch of points lies inside an element. For bilinearly interpolated elements.
  * \note For 2D quadrilateral elements this function is only an approximation. It is correct
@@ -774,14 +776,17 @@ t8_forest_element_point_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t
  *                              lies within an \a element, false otherwise. The return value is also true if the point 
  *                              lies on the element boundary. Thus, this function may return true for different leaf 
  *                              elements, if they are neighbors and the point lies on the common boundary.
+ * \param [in]      geom_is_axis_aligned Flag to tell if the used geometry is an axis-aligned geometry. Simplifies the
+ *                             check for this type of geometry. Can only be used for lines, quads and hexs. 
  * \param [in]      tolerance   Tolerance that we allow the point to not exactly match the element.
  *                              If this value is larger we detect more points.
  *                              If it is zero we probably do not detect points even if they are inside
  *                              due to rounding errors.
  */
 void
 t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element,
-                                      const double *points, int num_points, int *is_inside, const double tolerance);
+                                      const double *points, int num_points, int *is_inside,
+                                      const int geom_is_axis_aligned, const double tolerance);
 
 /* TODO: if set level and partition/adapt/balance all give NULL, then
  * refine uniformly and partition/adapt/balance the unfiform forest. */

test/t8_geometry/t8_gtest_point_inside.cxx

@@ -70,7 +70,7 @@ TEST (t8_point_inside, test_point_inside_specific_triangle)
 
   t8_element_t *element = t8_forest_get_element (forest, 0, NULL);
 
-  const int point_is_inside = t8_forest_element_point_inside (forest, 0, element, test_point, tolerance);
+  const int point_is_inside = t8_forest_element_point_inside (forest, 0, element, test_point, 0, tolerance);
   ASSERT_FALSE (point_is_inside) << "The point is wrongly detected as inside the triangle.";
   t8_forest_unref (&forest);
 }
@@ -112,30 +112,50 @@ TEST (t8_point_inside, test_point_inside_specific_quad)
 
   t8_element_t *element = t8_forest_get_element (forest, 0, NULL);
 
-  const int point_is_inside = t8_forest_element_point_inside (forest, 0, element, test_point, tolerance);
+  const int point_is_inside = t8_forest_element_point_inside (forest, 0, element, test_point, 0, tolerance);
 
   ASSERT_FALSE (point_is_inside) << "The point is wrongly detected as inside the quad.";
 
   t8_forest_unref (&forest);
 }
 
-class geometry_point_inside: public testing::TestWithParam<std::tuple<t8_eclass, int>> {
+/* *INDENT-OFF* */
+class geometry_point_inside: public testing::TestWithParam<std::tuple<t8_eclass, int, int>> {
  protected:
   void
   SetUp () override
   {
     eclass = std::get<0> (GetParam ());
     level = std::get<1> (GetParam ());
+    use_axis_aligned_geom = std::get<2> (GetParam ());
 
     /* Construct a cube coarse mesh */
-    cmesh = t8_cmesh_new_from_class (eclass, sc_MPI_COMM_WORLD);
+    if (use_axis_aligned_geom && (eclass == T8_ECLASS_LINE || eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_HEX)) {
+      /* clang-format off */
+      const double boundaries[24] = { 
+        0, 0, 0, 
+        1, 0, 0, 
+        0, 1, 0, 
+        1, 1, 0, 
+        0, 0, 1, 
+        1, 0, 1, 
+        0, 1, 1,
+        1, 1, 1 
+      };
+      /* clang-format on */
+      cmesh = t8_cmesh_new_hypercube_pad (eclass, sc_MPI_COMM_WORLD, boundaries, 1, 1, 1, use_axis_aligned_geom);
+    }
+    else {
+      cmesh = t8_cmesh_new_from_class (eclass, sc_MPI_COMM_WORLD);
+    }
   }
   void
   TearDown () override
   {
   }
   t8_eclass_t eclass;
   int level;
+  int use_axis_aligned_geom;
   t8_cmesh_t cmesh;
 };
 
@@ -282,7 +302,7 @@ TEST_P (geometry_point_inside, test_point_inside)
       /* We now check whether the point inside function correctly sees whether
          * the point is inside the element or not. */
       t8_forest_element_point_batch_inside (forest, 0, element, test_point, num_points, point_is_recognized_as_inside,
-                                            tolerance);
+                                            use_axis_aligned_geom, tolerance);
       for (int ipoint = 0; ipoint < num_points; ipoint++) {
         ASSERT_EQ (!point_is_recognized_as_inside[ipoint], !point_is_inside[ipoint])
           << "Testing point #" << ipoint << "(" << test_point[0] << "," << test_point[1] << "," << test_point[2]
@@ -302,8 +322,10 @@ TEST_P (geometry_point_inside, test_point_inside)
 
 #if T8_ENABLE_LESS_TESTS
 INSTANTIATE_TEST_SUITE_P (t8_gtest_point_inside, geometry_point_inside,
-                          testing::Combine (testing::Range (T8_ECLASS_LINE, T8_ECLASS_COUNT), testing::Range (0, 4)));
+                          testing::Combine (testing::Range (T8_ECLASS_LINE, T8_ECLASS_COUNT), testing::Range (0, 4),
+                                            testing::Range (0,1)));
 #else
 INSTANTIATE_TEST_SUITE_P (t8_gtest_point_inside, geometry_point_inside,
-                          testing::Combine (testing::Range (T8_ECLASS_LINE, T8_ECLASS_COUNT), testing::Range (0, 6)));
+                          testing::Combine (testing::Range (T8_ECLASS_LINE, T8_ECLASS_COUNT), testing::Range (0, 6),
+                                            testing::Range (0, 1)));
 #endif

tutorials/general/t8_tutorial_search.cxx

@@ -192,7 +192,7 @@ t8_tutorial_search_query_callback (t8_forest_t forest, t8_locidx_t ltreeid, cons
 
   /* Test whether this particle is inside this element. */
   particle_is_inside_element
-    = t8_forest_element_point_inside (forest, ltreeid, element, particle->coordinates, tolerance);
+    = t8_forest_element_point_inside (forest, ltreeid, element, particle->coordinates, 1, tolerance);
   if (particle_is_inside_element) {
     if (is_leaf) {
       /* The particle is inside and this element is a leaf element.