Feature point inside check for axis aligned geom

src/t8_cmesh/t8_cmesh.c

@@ -24,6 +24,7 @@
 #include <t8_cmesh.h>
 #include <t8_cmesh/t8_cmesh_geometry.h>
 #include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.h>
+#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h>
 #include <t8_refcount.h>
 #include <t8_data/t8_shmem.h>
 #include <t8_vec.h>
@@ -533,7 +534,14 @@ t8_cmesh_no_negative_volume (t8_cmesh_t cmesh)
     if (vertices != NULL) {
       /* Vertices are set */
       eclass = t8_cmesh_get_tree_class (cmesh, itree);
-      ret = t8_cmesh_tree_vertices_negative_volume (eclass, vertices, t8_eclass_num_vertices[eclass]);
+      const t8_gloidx_t gtree_id = t8_cmesh_get_global_id (cmesh, itree);
+      if (t8_geometry_get_type (cmesh, gtree_id) == T8_GEOMETRY_TYPE_LINEAR_AXIS_ALIGNED) {
+        /* Tree has negative volume if the diagonal goes from v_max to v_min and not vice versa */
+        ret = vertices[3] < vertices[0] && vertices[4] < vertices[1] && vertices[5] < vertices[2];
+      }
+      else {
+        ret = t8_cmesh_tree_vertices_negative_volume (eclass, vertices, t8_eclass_num_vertices[eclass]);
+      }
       if (ret) {
         t8_debugf ("Detected negative volume in tree %li\n", (long) itree);
       }

src/t8_cmesh/t8_cmesh_examples.c

@@ -987,7 +987,6 @@ t8_cmesh_set_vertices_2D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub
    */
   for (t8_locidx_t quad_y_id = 0; quad_y_id < quads_y; quad_y_id++) {
     for (t8_locidx_t quad_x_id = 0; quad_x_id < quads_x; quad_x_id++) {
-
       memcpy (vertices, box, 3 * sizeof (double));    /* Vertex 0 */
       t8_vec_axpyz (box, box_dir, vertices + 6, 1.0); /* Vertex 2 */
 
@@ -1316,9 +1315,9 @@ t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const do
 
     double vertices[6];
     /* Set first vertex to lower end of line */
-    t8_vec_axy (boundary, vertices, 1.0);
+    memcpy (vertices, boundary, 3 * sizeof (double));
     /* Set second vertex to lower end of line + line_dir */
-    t8_vec_axpyz (vertices + 3, boundary, line_dir, 1.0);
+    t8_vec_axpyz (line_dir, boundary, vertices + 3, 1.0);
 
     for (t8_gloidx_t tree_x = 0; tree_x < polygons_x; tree_x++) {
       t8_cmesh_set_tree_vertices (cmesh, tree_x, vertices, 2);

src/t8_forest/t8_forest_cxx.cxx

@@ -38,6 +38,7 @@
 #include <t8_geometry/t8_geometry_base.hxx>
 #if T8_ENABLE_DEBUG
 #include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.h>
+#include <t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h>
 #endif
 
 /* We want to export the whole implementation to be callable from "C" */
@@ -1088,325 +1089,16 @@ t8_forest_element_face_normal (t8_forest_t forest, t8_locidx_t ltreeid, const t8
   }
 }
 
-/**
- * Check if a point lies inside a vertex
- * 
- * \param[in] vertex_coords The coordinates of the vertex
- * \param[in] point         The coordinates of the point to check
- * \param[in] tolerance     A double > 0 defining the tolerance
- * \return                  0 if the point is outside, 1 otherwise.  
- */
-static int
-t8_vertex_point_inside (const double vertex_coords[3], const double point[3], const double tolerance)
-{
-  T8_ASSERT (tolerance > 0);
-  if (t8_vec_dist (vertex_coords, point) > tolerance) {
-    return 0;
-  }
-  return 1;
-}
-
-/**
- * Check if a point is inside a line that is defined by a starting point \a p_0
- * and a vector \a vec
- * 
- * \param[in] p_0         Starting point of the line
- * \param[in] vec         Direction of the line (not normalized)
- * \param[in] point       The coordinates of the point to check
- * \param[in] tolerance   A double > 0 defining the tolerance
- * \return                0 if the point is outside, 1 otherwise.  
- */
-static int
-t8_line_point_inside (const double *p_0, const double *vec, const double *point, const double tolerance)
-{
-  T8_ASSERT (tolerance > 0);
-  double b[3];
-  /* b = p - p_0 */
-  t8_vec_axpyz (p_0, point, b, -1);
-  double x = 0; /* Initialized to prevent compiler warning. */
-  int i;
-  /* So x is the solution to
-  * vec * x = b.
-  * We can compute it as
-  * x = b[i] / vec[i]
-  * if any vec[i] is not 0.
-  *
-  * Otherwise the line is degenerated (which should not happen).
-  */
-  for (i = 0; i < 3; ++i) {
-    if (vec[i] != 0) {
-      x = b[i] / vec[i];
-      break; /* found a non-zero coordinate. We can stop now. */
-    }
-  }
-
-  /* If i == 3 here, then vec = 0 and hence the line is degenerated. */
-  SC_CHECK_ABORT (i < 3, "Degenerated line element. Both endpoints are the same.");
-
-  if (x < -tolerance || x > 1 + tolerance) {
-    /* x is not an admissible solution. */
-    return 0;
-  }
-
-  /* we can check whether x gives us a solution by
-     * checking whether
-     *  vec * x = b
-     * is actually true.
-     */
-  double vec_check[3] = { vec[0], vec[1], vec[2] };
-  t8_vec_ax (vec_check, x);
-  if (t8_vec_dist (vec_check, b) > tolerance) {
-    /* Point does not lie on the line. */
-    return 0;
-  }
-  /* The point is on the line. */
-  return 1;
-}
-
-/**
- * Check if a point is inside of a triangle described by a point \a p_0 and two vectors \a v and \a w. 
- * 
- * \param[in] p_0         The first vertex of a triangle
- * \param[in] v           The vector from p_0 to p_1 (second vertex in the triangle)
- * \param[in] w           The vector from p_0 to p_2 (third vertex in the triangle)
- * \param[in] point       The coordinates of the point to check
- * \param[in] tolerance   A double > 0 defining the tolerance
- * \return                0 if the point is outside, 1 otherwise.  
- */
-static int
-t8_triangle_point_inside (const double p_0[3], const double v[3], const double w[3], const double point[3],
-                          const double tolerance)
-{
-  /* A point p is inside the triangle that is spanned
-   * by the point p_0 and vectors v and w if and only if the linear system
-   * vx + wy = point - p_0
-   * has a solution with 0 <= x,y and x + y <= 1.
-   *
-   * We check whether such a solution exists by computing
-   * certain determinants of 2x2 submatrizes of the 3x3 matrix
-   *
-   *  | v w e_3 | with v = p_1 - p_0, w = p_2 - p_0, and e_3 = (0 0 1)^t (third unit vector)
-   */
-
-  T8_ASSERT (tolerance > 0); /* negative values and zero are not allowed */
-  double b[3];
-  /* b = point - p_0 */
-  t8_vec_axpyz (p_0, point, b, -1);
-
-  /* Let d = det (v w e_3) */
-  const double det_vwe3 = v[0] * w[1] - v[1] * w[0];
-
-  /* The system has a solution, we need to compute it and
-   * check whether 0 <= x,y and x + y <= 1 */
-  /* x = det (b w e_3) / d
-   * y = det (v b e_3) / d
-   */
-  const double x = (b[0] * w[1] - b[1] * w[0]) / det_vwe3;
-  const double y = (v[0] * b[1] - v[1] * b[0]) / det_vwe3;
-
-  if (x < -tolerance || y < -tolerance || x + y > 1 + tolerance) {
-    /* The solution is not admissible.
-     * x < 0 or y < 0 or x + y > 1 */
-    return 0;
-  }
-  /* The solution may be admissible, but we have to
-   * check whether the result of
-   *  (p_1 - p_0)x + (p_2 - p_0)y ( = vx + wy)
-   * is actually p - p_0.
-   * Since the system of equations is overrepresented (3 equations, 2 variables)
-   * this may actually break.
-   * If it breaks, it will break in the z coordinate of the result.
-   */
-  const double z = v[2] * x + w[2] * y;
-  /* Must match the last coordinate of b = p - p_0 */
-  if (fabs (z - b[2]) > tolerance) {
-    /* Does not match. Point lies outside. */
-    return 0;
-  }
-  /* All checks passed. Point lies inside. */
-  return 1;
-}
-
-/** Check if a point lays on the inner side of a plane of a bilinearly interpolated volume element. 
- * the plane is described by a point and the normal of the face. 
- * \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
- * \return                    0 if the point is outside, 1 otherwise.                   
- */
-static int
-t8_plane_point_inside (const double point_on_face[3], const double face_normal[3], const double point[3])
-{
-  /* Set x = x - p */
-  double pof[3] = { point_on_face[0], point_on_face[1], point_on_face[2] };
-  t8_vec_axpy (point, pof, -1);
-  /* Compute <x-p,n> */
-  const double dot_product = t8_vec_dot (pof, face_normal);
-  if (dot_product < 0) {
-    /* The point is on the wrong side of the plane */
-    return 0;
-  }
-  return 1;
-}
-
 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)
+t8_forest_element_points_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 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);
-  const t8_element_shape_t element_shape = ts->t8_element_shape (element);
-
-#if T8_ENABLE_DEBUG
   /* Check whether the provided geometry is linear */
   const t8_cmesh_t cmesh = t8_forest_get_cmesh (forest);
   const t8_locidx_t cltreeid = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid);
   const t8_gloidx_t cgtreeid = t8_cmesh_get_global_id (cmesh, cltreeid);
-  T8_ASSERT (t8_geometry_get_type (cmesh, cgtreeid) == T8_GEOMETRY_TYPE_LINEAR);
-#endif
-
-  switch (element_shape) {
-  case T8_ECLASS_VERTEX: {
-    /* A point is 'inside' a vertex if they have the same coordinates */
-    double vertex_coords[3];
-    /* Get the vertex coordinates */
-    t8_forest_element_coordinate (forest, ltreeid, element, 0, vertex_coords);
-    /* Check whether the point and the vertex are within tolerance distance
-       * to each other */
-    for (int ipoint = 0; ipoint < num_points; ipoint++) {
-      is_inside[ipoint] = t8_vertex_point_inside (vertex_coords, &points[ipoint * 3], tolerance);
-    }
-    return;
-  }
-  case T8_ECLASS_LINE: {
-    /* A point p is inside a line that is defined by the edge nodes
-     * p_0 and p_1
-     * if and only if the linear system
-     * (p_1 - p_0)x = p - p_0
-     * has a solution x with 0 <= x <= 1
-     */
-    double p_0[3], v[3];
-
-    /* Compute the vertex coordinates of the line */
-    t8_forest_element_coordinate (forest, ltreeid, element, 0, p_0);
-    /* v = p_1 */
-    t8_forest_element_coordinate (forest, ltreeid, element, 1, v);
-    /* v = p_1 - p_0 */
-    t8_vec_axpy (p_0, v, -1);
-    for (int ipoint = 0; ipoint < num_points; ipoint++) {
-      is_inside[ipoint] = t8_line_point_inside (p_0, v, &points[ipoint * 3], tolerance);
-    }
-    return;
-  }
-  case T8_ECLASS_QUAD: {
-    /* We divide the quad in two triangles and use the triangle check. */
-    double p_0[3], p_1[3], p_2[3], p_3[3];
-    /* Compute the vertex coordinates of the quad */
-    t8_forest_element_coordinate (forest, ltreeid, element, 0, p_0);
-    t8_forest_element_coordinate (forest, ltreeid, element, 1, p_1);
-    t8_forest_element_coordinate (forest, ltreeid, element, 2, p_2);
-    t8_forest_element_coordinate (forest, ltreeid, element, 3, p_3);
-
-#if T8_ENABLE_DEBUG
-    /* Issue a warning if the points of the quad do not lie in the same plane */
-    if (!t8_four_points_coplanar (p_0, p_1, p_2, p_3, tolerance)) {
-      t8_debugf ("WARNING: Testing if point is inside a quad that is not coplanar. This test will be inaccurate.\n");
-    }
-#endif
-    double v[3];
-    double w[3];
-    /* v = v - p_0 = p_1 - p_0 */
-    t8_vec_axpyz (p_0, p_1, v, -1);
-    /* w = w - p_0 = p_2 - p_0 */
-    t8_vec_axpyz (p_0, p_2, w, -1);
-    /* Check whether the point is inside the first triangle. */
-    for (int ipoint = 0; ipoint < num_points; ipoint++) {
-      is_inside[ipoint] = t8_triangle_point_inside (p_0, v, w, &points[ipoint * 3], tolerance);
-    }
-    /* If not, check whether the point is inside the second triangle. */
-    /* v = v - p_0 = p_1 - p_0 */
-    t8_vec_axpyz (p_1, p_2, v, -1);
-    /* w = w - p_0 = p_2 - p_0 */
-    t8_vec_axpyz (p_1, p_3, w, -1);
-    for (int ipoint = 0; ipoint < num_points; ipoint++) {
-      if (!is_inside[ipoint]) {
-        /* point_inside is true if the point was inside the first or second triangle. Otherwise it is false. */
-        is_inside[ipoint] = t8_triangle_point_inside (p_1, v, w, &points[ipoint * 3], tolerance);
-      }
-    }
-    return;
-  }
-  case T8_ECLASS_TRIANGLE: {
-    double p_0[3], p_1[3], p_2[3];
-
-    /* Compute the vertex coordinates of the triangle */
-    t8_forest_element_coordinate (forest, ltreeid, element, 0, p_0);
-    t8_forest_element_coordinate (forest, ltreeid, element, 1, p_1);
-    t8_forest_element_coordinate (forest, ltreeid, element, 2, p_2);
-    double v[3];
-    double w[3];
-    /* v = v - p_0 = p_1 - p_0 */
-    t8_vec_axpyz (p_0, p_1, v, -1);
-    /* w = w - p_0 = p_2 - p_0 */
-    t8_vec_axpyz (p_0, p_2, w, -1);
-
-    for (int ipoint = 0; ipoint < num_points; ipoint++) {
-      is_inside[ipoint] = t8_triangle_point_inside (p_0, v, w, &points[ipoint * 3], tolerance);
-    }
-    return;
-  }
-  case T8_ECLASS_TET:
-  case T8_ECLASS_HEX:
-  case T8_ECLASS_PRISM:
-  case T8_ECLASS_PYRAMID: {
-    /* For bilinearly interpolated volume elements, a point is inside an element
-     * if and only if it lies on the inner side of each face.
-     * The inner side is defined as the side where the outside normal vector does not
-     * point to.
-     * The point is on this inner side if and only if the scalar product of
-     * a point on the plane minus the point with the outer normal of the face
-     * is >= 0.
-     *
-     * In other words, let p be the point to check, n the outer normal and x a point
-     * on the plane, then p is on the inner side if and only if
-     *  <x - p, n> >= 0
-     */
-
-    const int num_faces = ts->t8_element_num_faces (element);
-    /* Assume that every point is inside of the element */
-    for (int ipoint = 0; ipoint < num_points; ipoint++) {
-      is_inside[ipoint] = 1;
-    }
-    for (int iface = 0; iface < num_faces; ++iface) {
-      double face_normal[3];
-      /* Compute the outer normal n of the face */
-      t8_forest_element_face_normal (forest, ltreeid, element, iface, face_normal);
-      /* Compute a point x on the face */
-      const int afacecorner = ts->t8_element_get_face_corner (element, iface, 0);
-      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]);
-        if (is_inside_iface == 0) {
-          /* Point is on the outside of face iface. Update is_inside */
-          is_inside[ipoint] = 0;
-        }
-      }
-    }
-    return;
-  }
-  default:
-    SC_ABORT_NOT_REACHED ();
-  }
-}
-
-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)
-{
-  int is_inside = 0;
-  t8_forest_element_point_batch_inside (forest, ltreeid, element, point, 1, &is_inside, tolerance);
-  return is_inside;
+  const t8_geometry_c *geometry = t8_cmesh_get_tree_geometry (cmesh, cgtreeid);
+  geometry->t8_geom_point_batch_inside_element (forest, ltreeid, element, points, num_points, is_inside, tolerance);
 }
 
 /* For each tree in a forest compute its first and last descendant */