Improvement: Cad geometry component test

src/t8_forest/t8_forest.cxx

@@ -49,6 +49,7 @@
 #include <t8_data/t8_element_array_iterator.hxx>
 
 #include <algorithm>
+#include <span>
 
 /* We want to export the whole implementation to be callable from "C" */
 T8_EXTERN_C_BEGIN ();
@@ -511,6 +512,33 @@ t8_forest_element_centroid (t8_forest_t forest, t8_locidx_t ltreeid, const t8_el
                                      coordinates);
 }
 
+/* Compute the center of mass of an element. We can use the element reference
+ * coordinates of the centroid.*/
+void
+t8_forest_element_linear_centroid (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element,
+                                   double *coordinates_c)
+{
+  T8_ASSERT (t8_forest_is_committed (forest));
+  const t8_scheme *scheme = t8_forest_get_scheme (forest);
+  const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid);
+  std::span<double, 3> coordinates = std::span<double, 3> (coordinates_c, 3);
+  std::fill (coordinates.begin (), coordinates.end (), 0);
+
+  /* Get the tree's eclass and scheme. */
+  T8_ASSERT (scheme->element_is_valid (tree_class, element));
+
+  /* Get the element class and calculate the centroid using its corners. The centroid is
+    the sum of all corner coordinates divided by the number of corners. */
+  const t8_element_shape_t element_shape = scheme->element_get_shape (tree_class, element);
+  const int num_corners = t8_eclass_num_vertices[element_shape];
+  std::array<double, 3> corner {};
+  for (int icorner = 0; icorner < num_corners; ++icorner) {
+    t8_forest_element_coordinate (forest, ltreeid, element, icorner, corner.data ());
+    t8_axpy (corner.data (), coordinates.data (), 1);
+  }
+  t8_ax (coordinates.data (), 1.0 / num_corners);
+}
+
 /* Compute the length of the line from one corner to a second corner in an element */
 static double
 t8_forest_element_line_length (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, int corner_a,
@@ -871,52 +899,6 @@ t8_forest_element_face_centroid (t8_forest_t forest, t8_locidx_t ltreeid, const
   }
 }
 
-#if T8_ENABLE_DEBUG
-/* Test whether four given points in 3D are coplanar up to a given tolerance.
- */
-static int
-t8_four_points_coplanar (const double p_0[3], const double p_1[3], const double p_2[3], const double p_3[3],
-                         const double tolerance)
-{
-  /* Let p0, p1, p2, p3 be the four points.
-   * The four points are coplanar if the normal vectors to the triangles
-   * p0, p1, p2 and p0, p2, p3 are pointing in the same direction.
-   *
-   * We build the vectors A = p1 - p0, B = p2 - p0 and C = p3 - p0.
-   * The normal vectors to the triangles are n1 = A x B and n2 = A x C.
-   * These are pointing in the same direction if their cross product is 0.
-   * Hence we check if || n1 x n2 || < tolerance. */
-
-  /* A = p1 - p0 */
-  double A[3];
-  t8_axpyz (p_0, p_1, A, -1);
-
-  /* B = p2 - p0 */
-  double B[3];
-  t8_axpyz (p_0, p_2, B, -1);
-
-  /* C = p3 - p0 */
-  double C[3];
-  t8_axpyz (p_0, p_3, C, -1);
-
-  /* n1 = A x B */
-  double A_cross_B[3];
-  t8_cross_3D (A, B, A_cross_B);
-
-  /* n2 = A x C */
-  double A_cross_C[3];
-  t8_cross_3D (A, C, A_cross_C);
-
-  /* n1 x n2 */
-  double n1_cross_n2[3];
-  t8_cross_3D (A_cross_B, A_cross_C, n1_cross_n2);
-
-  /* || n1 x n2 || */
-  const double norm = t8_norm (n1_cross_n2);
-  return norm < tolerance;
-}
-#endif
-
 void
 t8_forest_element_face_normal (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, int face,
                                double normal[3])
@@ -984,7 +966,7 @@ t8_forest_element_face_normal (t8_forest_t forest, t8_locidx_t ltreeid, const t8
     t8_forest_element_coordinate (forest, ltreeid, element, corner_a, vertex_a);
     t8_forest_element_coordinate (forest, ltreeid, element, corner_b, vertex_b);
     /* Compute the center */
-    t8_forest_element_centroid (forest, ltreeid, element, center);
+    t8_forest_element_linear_centroid (forest, ltreeid, element, center);
 
     /* Compute the difference with V_a.
        * Compute the dot products */
@@ -1066,7 +1048,7 @@ t8_forest_element_face_normal (t8_forest_t forest, t8_locidx_t ltreeid, const t8
     norm = t8_norm (normal);
     T8_ASSERT (norm > 1e-14);
     /* Compute the coordinates of the center of the element */
-    t8_forest_element_centroid (forest, ltreeid, element, center);
+    t8_forest_element_linear_centroid (forest, ltreeid, element, center);
     /* Compute center = center - vertex_0 */
     t8_axpy (corner_vertices[0], center, -1);
     /* Compute the dot-product of normal and center */

src/t8_forest/t8_forest_geometrical.h

@@ -57,7 +57,7 @@ t8_forest_element_coordinate (t8_forest_t forest, t8_locidx_t ltree_id, const t8
  *  The point is given in reference coordinates inside the element and gets
  *  converted to reference coordinates inside the tree. After that, the point
  *  is converted to global coordinates inside the domain. If needed, the element
- *  is stretched by the given stretch factors (the resulting mesh is then 
+ *  is stretched by the given stretch factors (the resulting mesh is then
  *  no longer non-overlapping).
  * \param [in]      forest            The forest.
  * \param [in]      ltreeid           The forest local id of the tree in which the element is.
@@ -104,6 +104,20 @@ t8_forest_element_from_ref_coords (t8_forest_t forest, t8_locidx_t ltreeid, cons
 void
 t8_forest_element_centroid (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, double *coordinates);
 
+/** Compute the coordinates of the centroid of an element by its corner coordinates.
+ * This treats every element as a linear element. The centroid is the mean of all the corner coordinates.
+ * \warning This function omits if an element is curved. So for linear elements it produces the same results as
+ * \ref t8_forest_element_centroid. For curved elements the curvature is ignored.
+ * \param [in]      forest     The forest.
+ * \param [in]      ltreeid    The forest local id of the tree in which the element is.
+ * \param [in]      element    The element.
+ * \param [out]     coordinates On input an allocated array to store 3 doubles, on output
+ *                             the x, y and z coordinates of the centroid.
+ */
+void
+t8_forest_element_linear_centroid (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element,
+                                   double *coordinates);
+
 /** Compute the diameter of an element if a geometry for this tree is registered in the forest's cmesh.
  * This is only an approximation.
  * \param [in]      forest     The forest.

src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.cxx

@@ -21,7 +21,7 @@
 */
 
 /** \file t8_geometry_cad.cxx
- * Implements functions declared in \ref t8_geometry_cad.hxx 
+ * Implements functions declared in \ref t8_geometry_cad.hxx
  *  or the C interface \ref t8_geometry_cad.h.
  */
 
@@ -31,6 +31,9 @@
 #include <t8_eclass.h>
 #include <t8_geometry/t8_geometry_helpers.h>
 #include <t8_schemes/t8_default/t8_default_prism/t8_dprism.h>
+#include <t8_schemes/t8_scheme.hxx>
+#include <t8_types/t8_vec.h>
+#include <t8_types/t8_vec.hxx>
 
 #include <BRep_Builder.hxx>
 #include <BRep_Tool.hxx>
@@ -138,6 +141,150 @@ t8_geometry_cad::t8_geom_load_tree_data (t8_cmesh_t cmesh, t8_gloidx_t gtreeid)
   T8_ASSERT (faces != NULL);
 }
 
+void
+t8_geometry_cad::t8_geom_point_batch_inside_element (t8_forest_t forest, t8_locidx_t ltreeid,
+                                                     const t8_element_t *element, const double *points,
+                                                     const int num_points, int *is_inside, const double tolerance) const
+{
+  const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid);
+  const t8_scheme *scheme = t8_forest_get_scheme (forest);
+  const t8_element_shape_t element_shape = scheme->element_get_shape (tree_class, element);
+  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
+     */
+    t8_3D_vec p_0, v;
+
+    /* Compute the vertex coordinates of the line */
+    t8_forest_element_coordinate (forest, ltreeid, element, 0, p_0.data ());
+    /* v = p_1 */
+    t8_forest_element_coordinate (forest, ltreeid, element, 1, v.data ());
+    /* v = p_1 - p_0 */
+    t8_axpy (p_0, v, -1);
+    for (int ipoint = 0; ipoint < num_points; ipoint++) {
+      is_inside[ipoint] = t8_line_point_inside (p_0.data (), v.data (), &points[ipoint * 3], tolerance);
+    }
+    return;
+  }
+  case T8_ECLASS_QUAD: {
+    /* We divide the quad in two triangles and use the triangle check. */
+    t8_3D_vec p_0, p_1, p_2, p_3;
+    /* Compute the vertex coordinates of the quad */
+    t8_forest_element_coordinate (forest, ltreeid, element, 0, p_0.data ());
+    t8_forest_element_coordinate (forest, ltreeid, element, 1, p_1.data ());
+    t8_forest_element_coordinate (forest, ltreeid, element, 2, p_2.data ());
+    t8_forest_element_coordinate (forest, ltreeid, element, 3, p_3.data ());
+
+#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.data (), p_1.data (), p_2.data (), p_3.data (), tolerance)) {
+      t8_debugf ("WARNING: Testing if point is inside a quad that is not coplanar. This test will be inaccurate.\n");
+    }
+#endif
+    t8_3D_vec v;
+    t8_3D_vec w;
+    /* v = v - p_0 = p_1 - p_0 */
+    t8_axpyz (p_0, p_1, v, -1);
+    /* w = w - p_0 = p_2 - p_0 */
+    t8_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.data (), v.data (), w.data (), &points[ipoint * 3], tolerance);
+    }
+    /* If not, check whether the point is inside the second triangle. */
+    /* v = v - p_0 = p_1 - p_0 */
+    t8_axpyz (p_1, p_2, v, -1);
+    /* w = w - p_0 = p_2 - p_0 */
+    t8_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.data (), v.data (), w.data (), &points[ipoint * 3], tolerance);
+      }
+    }
+    return;
+  }
+  case T8_ECLASS_TRIANGLE: {
+    t8_3D_vec p_0, p_1, p_2;
+
+    /* Compute the vertex coordinates of the triangle */
+    t8_forest_element_coordinate (forest, ltreeid, element, 0, p_0.data ());
+    t8_forest_element_coordinate (forest, ltreeid, element, 1, p_1.data ());
+    t8_forest_element_coordinate (forest, ltreeid, element, 2, p_2.data ());
+    t8_3D_vec v;
+    t8_3D_vec w;
+    /* v = v - p_0 = p_1 - p_0 */
+    t8_axpyz (p_0, p_1, v, -1);
+    /* w = w - p_0 = p_2 - p_0 */
+    t8_axpyz (p_0, p_2, w, -1);
+
+    for (int ipoint = 0; ipoint < num_points; ipoint++) {
+      is_inside[ipoint] = t8_triangle_point_inside (p_0.data (), v.data (), w.data (), &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 = scheme->element_get_num_faces (tree_class, 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 = scheme->element_get_face_corner (tree_class, 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 ();
+  }
+}
+
 void
 t8_geometry_cad::t8_geom_evaluate_cad_tri (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords,
                                            const size_t num_coords, double *out_coords) const

src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx

@@ -91,7 +91,7 @@ struct t8_geometry_cad: public t8_geometry_with_vertices
    * \return The type.
    */
   inline t8_geometry_type_t
-  t8_geom_get_type () const
+  t8_geom_get_type () const override
   {
     return T8_GEOMETRY_TYPE_CAD;
   };
@@ -106,7 +106,7 @@ struct t8_geometry_cad: public t8_geometry_with_vertices
    */
   virtual void
   t8_geom_evaluate (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords,
-                    double *out_coords) const;
+                    double *out_coords) const override;
 
   /**
    * Compute the jacobian of the \a t8_geom_evaluate map at a point in the reference space \f$ [0,1]^\mathrm{dim} \f$.
@@ -119,7 +119,7 @@ struct t8_geometry_cad: public t8_geometry_with_vertices
    */
   virtual void
   t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords,
-                             double *jacobian) const;
+                             double *jacobian) const override;
 
   /** Update a possible internal data buffer for per tree data.
    * This function is called before the first coordinates in a new tree are
@@ -129,22 +129,37 @@ struct t8_geometry_cad: public t8_geometry_with_vertices
    * \param [in]  gtreeid    The global tree.
    */
   virtual void
-  t8_geom_load_tree_data (t8_cmesh_t cmesh, t8_gloidx_t gtreeid);
+  t8_geom_load_tree_data (t8_cmesh_t cmesh, t8_gloidx_t gtreeid) override;
 
   /**
    * Check for compatibility of the currently loaded tree with the geometry.
    * This geometry supports all element types, hence it returns true.
    * \return                True if the geometry is compatible with the tree.
    */
   bool
-  t8_geom_check_tree_compatibility () const
+  t8_geom_check_tree_compatibility () const override
   {
     return true;
   }
 
+  /** Checks if points are inside the element. Input a list of points \a points
+   * and it returns a list of flags \a is_inside if the points are contained.
+   * \param[in] forest            The forest of the element.
+   * \param[in] ltreeid           The local tree id of the element's tree
+   * \param[in] element           The element
+   * \param[in] points            points to check
+   * \param[in] num_points        Number of points to check
+   * \param[in, out] is_inside    Array to fill with flags whether the point is inside or not
+   * \param[in] tolerance         Tolerance of the inside-check
+   */
+  virtual void
+  t8_geom_point_batch_inside_element (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element,
+                                      const double *points, const int num_points, int *is_inside,
+                                      const double tolerance) const override;
+
   /**
    * Getter function for the CAD manager.
-   * 
+   *
    * \return The CAD manager of the geometry.
   */
   std::shared_ptr<t8_cad>