Skip to content

Commit

Permalink
Merge pull request #751 from DLR-AMR/feature-point_inside_check_for_a…
Browse files Browse the repository at this point in the history 
…xis_aligned_geom

Feature point inside check for axis aligned geom
  • Loading branch information
@Davknapp
Davknapp committed Feb 14, 2024
2 parents b3ee013 + 7d7fe2f commit 800648f
Show file tree
Hide file tree
Showing 18 changed files with 588 additions and 371 deletions.
10 changes: 9 additions & 1 deletion src/t8_cmesh/t8_cmesh.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down Expand Up @@ -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);
}
Expand Down
5 changes: 2 additions & 3 deletions src/t8_cmesh/t8_cmesh_examples.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */

Expand Down Expand Up @@ -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);
Expand Down
318 changes: 5 additions & 313 deletions src/t8_forest/t8_forest_cxx.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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" */
Expand Down Expand Up @@ -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 */
Expand Down
Loading

0 comments on commit 800648f

Please sign in to comment.