Open Huite opened 1 year ago
This might be worthwhile when dealing with a fully triangular grid.
Coming from xugrid:
@nb.njit(inline="always") def point_in_triangle(p: Point, triangle: Triangle) -> bool: """Unrolled half-plane check.""" a = cross_product(to_vector(triangle.a, triangle.b), to_vector(triangle.a, p)) > 0 b = cross_product(to_vector(triangle.b, triangle.c), to_vector(triangle.b, p)) > 0 c = cross_product(to_vector(triangle.c, triangle.a), to_vector(triangle.c, p)) > 0 return (a == b) and (b == c) @nb.njit(inline="always", parallel=True, cache=True) def points_in_triangles( points: FloatArray, face_indices: IntArray, faces: IntArray, vertices: FloatArray, ): n_points = len(points) inside = np.empty(n_points, dtype=np.bool_) for i in nb.prange(n_points): face_index = face_indices[i] face = faces[face_index] triangle = as_triangle(vertices, face) point = as_point(points[i]) inside[i] = point_in_triangle(point, triangle) return inside
This might be worthwhile when dealing with a fully triangular grid.
Coming from xugrid: