From 89ac95566944b6833c7b5659b5516acea0d617e9 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Thu, 14 Sep 2023 12:57:01 +0200 Subject: [PATCH 01/23] Optimize t8_forest_element_point_inside for axis-aligned geometries --- src/t8_forest/t8_forest_cxx.cxx | 23 ++++++++++++++++++++--- src/t8_forest/t8_forest_general.h | 9 +++++++-- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/t8_forest/t8_forest_cxx.cxx b/src/t8_forest/t8_forest_cxx.cxx index 2e2faf96af..c851c7a569 100644 --- a/src/t8_forest/t8_forest_cxx.cxx +++ b/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,22 @@ 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, 0, 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*/ + } + } switch (element_shape) { case T8_ECLASS_VERTEX: { /* A point is 'inside' a vertex if they have the same coordinates */ @@ -1380,10 +1396,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; } diff --git a/src/t8_forest/t8_forest_general.h b/src/t8_forest/t8_forest_general.h index cab8656932..0fc70611db 100644 --- a/src/t8_forest/t8_forest_general.h +++ b/src/t8_forest/t8_forest_general.h @@ -752,6 +752,8 @@ 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 @@ -759,7 +761,7 @@ t8_forest_iterate (t8_forest_t forest); */ 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,6 +776,8 @@ 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 @@ -781,7 +785,8 @@ t8_forest_element_point_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t */ 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. */ From d7e6a883e37c2c79c955fb9fd6c1796c2e962cb1 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Thu, 14 Sep 2023 13:24:56 +0200 Subject: [PATCH 02/23] Added test to check axis_aligned feature of t8_forest_point_inside --- test/t8_geometry/t8_gtest_point_inside.cxx | 32 ++++++++++++++++++---- tutorials/general/t8_tutorial_search.cxx | 2 +- 2 files changed, 27 insertions(+), 7 deletions(-) diff --git a/test/t8_geometry/t8_gtest_point_inside.cxx b/test/t8_geometry/t8_gtest_point_inside.cxx index dfc85970c1..fa69906737 100644 --- a/test/t8_geometry/t8_gtest_point_inside.cxx +++ b/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,7 +112,7 @@ 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."; @@ -120,16 +120,34 @@ TEST (t8_point_inside, test_point_inside_specific_quad) } /* *INDENT-OFF* */ -class geometry_point_inside: public testing::TestWithParam> { +class geometry_point_inside: public testing::TestWithParam> { 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 @@ -137,6 +155,7 @@ class geometry_point_inside: public testing::TestWithParamcoordinates, 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. From fb4d634226522e254fc22b05f97576515287c953 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Thu, 14 Sep 2023 15:16:55 +0200 Subject: [PATCH 03/23] use second vertex for v_max --- src/t8_forest/t8_forest_cxx.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/t8_forest/t8_forest_cxx.cxx b/src/t8_forest/t8_forest_cxx.cxx index c851c7a569..45c4f38166 100644 --- a/src/t8_forest/t8_forest_cxx.cxx +++ b/src/t8_forest/t8_forest_cxx.cxx @@ -1249,7 +1249,7 @@ t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, c /*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, 0, v_max); + 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 */ @@ -1258,6 +1258,7 @@ t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, c 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: { From 533554af890236d68089f07522d1351d626455c1 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 15 Sep 2023 13:26:39 +0200 Subject: [PATCH 04/23] updated test for axis-algined geometries. TODO: Update cmesh_is_commited to be aware of geometry --- src/t8_cmesh/t8_cmesh_examples.c | 73 ++++++++++------------ test/t8_geometry/t8_gtest_point_inside.cxx | 42 ++++++++----- 2 files changed, 58 insertions(+), 57 deletions(-) diff --git a/src/t8_cmesh/t8_cmesh_examples.c b/src/t8_cmesh/t8_cmesh_examples.c index 6c91a04ec7..7c8968c8f8 100644 --- a/src/t8_cmesh/t8_cmesh_examples.c +++ b/src/t8_cmesh/t8_cmesh_examples.c @@ -941,29 +941,24 @@ 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 */ - if (use_axis_aligned_geom && eclass == T8_ECLASS_QUAD) { - t8_vec_axy (box, vertices, 1.0); /* Vertex 0*/ - t8_vec_axpyz (box, box_dir, vertices + 9, 1.0); /* Vertex 3 */ - } - else { - t8_vec_axy (box, vertices, 1.0); /* Vertex 0 */ - t8_vec_axpyz (box, box_dir, vertices + 6, 1.0); /* Vertex 2 */ + /* Reduce box along x axis */ + t8_resize_box (2, box, box_dir, 0, 1, box_quads); + t8_update_box_face_edges (2, box, box_dir, 0, box_quads); - /* Reduce box along x axis */ - t8_resize_box (2, box, box_dir, 0, 1, box_quads); - t8_update_box_face_edges (2, box, box_dir, 0, box_quads); - - t8_vec_axy (box, vertices + 3, 1.0); /* Vertex 1 */ - t8_vec_axpyz (box, box_dir, vertices + 9, 1.0); /* Vertex 3 */ + t8_vec_axy (box, vertices + 3, 1.0); /* Vertex 1 */ + t8_vec_axpyz (box, box_dir, vertices + 9, 1.0); /* Vertex 3 */ + if (use_axis_aligned_geom && eclass == T8_ECLASS_QUAD) { + memcpy (vertices + 3, vertices + 9, 3 * sizeof (double)); } /* Map vertices of current quad on to respective trees inside. */ if (eclass == T8_ECLASS_QUAD) { /* No mapping is required. */ const t8_locidx_t tree_id = quad_y_id * quads_x + quad_x_id; - t8_cmesh_set_tree_vertices (cmesh, tree_id, vertices, - (use_axis_aligned_geom && eclass == T8_ECLASS_QUAD) ? 2 : 4); + t8_cmesh_set_tree_vertices (cmesh, tree_id, vertices, (use_axis_aligned_geom) ? 2 : 4); } else { T8_ASSERT (eclass == T8_ECLASS_TRIANGLE); @@ -1088,42 +1083,38 @@ t8_cmesh_set_vertices_3D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub for (t8_locidx_t hex_z_id = 0; hex_z_id < hexs_z; hex_z_id++) { for (t8_locidx_t hex_y_id = 0; hex_y_id < hexs_y; hex_y_id++) { for (t8_locidx_t hex_x_id = 0; hex_x_id < hexs_x; hex_x_id++) { - if (use_axis_aligned_geom && eclass == T8_ECLASS_HEX) { - t8_vec_axy (box, vertices, 1.0); /* Vertex 0 */ - t8_vec_axpyz (box, box_dir + 12, vertices + 21, 1.0); /* Vertex 7 */ - } - else { - t8_vec_axy (box, vertices, 1.0); /* Vertex 0 */ - t8_vec_axpyz (box, box_dir + 12, vertices + 6, 1.0); /* Vertex 2 */ + memcpy (vertices, box, 3 * sizeof (double)); /* Vertex 0 */ + t8_vec_axpyz (box, box_dir + 12, vertices + 6, 1.0); /* Vertex 2 */ - /* Reduce box along z axis and face 4. */ - t8_resize_box (3, box, box_dir, 4, 1, box_hexs); - t8_update_box_face_edges (3, box, box_dir, 4, box_hexs); + /* Reduce box along z axis and face 4. */ + t8_resize_box (3, box, box_dir, 4, 1, box_hexs); + t8_update_box_face_edges (3, box, box_dir, 4, box_hexs); - t8_vec_axy (box, vertices + 12, 1.0); /* Vertex 4 */ - t8_vec_axpyz (box, box_dir + 12, vertices + 18, 1.0); /* Vertex 6 */ + t8_vec_axy (box, vertices + 12, 1.0); /* Vertex 4 */ + t8_vec_axpyz (box, box_dir + 12, vertices + 18, 1.0); /* Vertex 6 */ - /* Reduce box along x axis and face 0. */ - t8_resize_box (3, box, box_dir, 0, 1, box_hexs); - t8_update_box_face_edges (3, box, box_dir, 0, box_hexs); + /* Reduce box along x axis and face 0. */ + t8_resize_box (3, box, box_dir, 0, 1, box_hexs); + t8_update_box_face_edges (3, box, box_dir, 0, box_hexs); - t8_vec_axy (box, vertices + 15, 1.0); /* Vertex 5 */ - t8_vec_axpyz (box, box_dir + 12, vertices + 21, 1.0); /* Vertex 7 */ + t8_vec_axy (box, vertices + 15, 1.0); /* Vertex 5 */ + t8_vec_axpyz (box, box_dir + 12, vertices + 21, 1.0); /* Vertex 7 */ - /* Increase box along z axis and and face 4 */ - t8_resize_box (3, box, box_dir, 4, -1, box_hexs); - t8_update_box_face_edges (3, box, box_dir, 4, box_hexs); + /* Increase box along z axis and and face 4 */ + t8_resize_box (3, box, box_dir, 4, -1, box_hexs); + t8_update_box_face_edges (3, box, box_dir, 4, box_hexs); - t8_vec_axy (box, vertices + 3, 1.0); /* Vertex 1 */ - t8_vec_axpyz (box, box_dir + 12, vertices + 9, 1.0); /* Vertex 3 */ + t8_vec_axy (box, vertices + 3, 1.0); /* Vertex 1 */ + t8_vec_axpyz (box, box_dir + 12, vertices + 9, 1.0); /* Vertex 3 */ + if (use_axis_aligned_geom && eclass == T8_ECLASS_HEX) { + memcpy (vertices + 3, vertices + 21, 3 * sizeof (double)); } /* Map vertices of current hex on to respective trees inside. */ const t8_locidx_t hex_id = hex_z_id * hexs_y * hexs_x + hex_y_id * hexs_x + hex_x_id; if (eclass == T8_ECLASS_HEX) { /* No mapping is required. */ - t8_cmesh_set_tree_vertices (cmesh, hex_id, vertices, - (use_axis_aligned_geom && eclass == T8_ECLASS_HEX) ? 2 : 8); + t8_cmesh_set_tree_vertices (cmesh, hex_id, vertices, use_axis_aligned_geom ? 2 : 8); } else if (eclass == T8_ECLASS_TET) { const t8_locidx_t tree_id_0 = hex_id * 6; @@ -1274,9 +1265,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_locidx_t tree_x = 0; tree_x < polygons_x; tree_x++) { t8_cmesh_set_tree_vertices (cmesh, tree_x, vertices, 2); diff --git a/test/t8_geometry/t8_gtest_point_inside.cxx b/test/t8_geometry/t8_gtest_point_inside.cxx index 58fe21ff53..44d175d6e4 100644 --- a/test/t8_geometry/t8_gtest_point_inside.cxx +++ b/test/t8_geometry/t8_gtest_point_inside.cxx @@ -33,6 +33,7 @@ #include #include #include +#include /* In this test we define a triangle in the x-y plane * and a point that lies in a triangle that is parallel @@ -180,23 +181,26 @@ TEST_P (geometry_point_inside, test_point_inside) * finding function that may be biased towards certrain coordinate axis. */ double *const tree_vertices = t8_cmesh_get_tree_vertices (cmesh, 0); - const int num_vertices = t8_eclass_num_vertices[eclass]; + const int num_vertices = (use_axis_aligned_geom && (eclass == T8_ECLASS_HEX || eclass == T8_ECLASS_QUAD)) + ? 2 + : t8_eclass_num_vertices[eclass]; /* Translate all points by the same vector to move the element a bit. */ - double translate_all_points[3] = { -0.1, 0.3, 0.15 }; + double translate_all_points[3] = { 1, 0, 0 }; t8_cmesh_translate_coordinates (tree_vertices, tree_vertices, num_vertices, translate_all_points); /* Translate points 0 and 1 (if it exists) extra in order to move the 2D elements * and 3D faces outside of axis perpendicular planes. */ - double translate_points_0_1[3] = { 0.1, -0.1, 0.3 }; - t8_cmesh_translate_coordinates (tree_vertices, tree_vertices, 1, translate_points_0_1); - if (num_vertices > 2) { - t8_cmesh_translate_coordinates (tree_vertices + 3, tree_vertices + 3, 1, translate_points_0_1); + if (!use_axis_aligned_geom) { + double translate_points_0_1[3] = { 0.1, -0.1, 0.3 }; + t8_cmesh_translate_coordinates (tree_vertices, tree_vertices, 1, translate_points_0_1); + if (num_vertices > 2) { + t8_cmesh_translate_coordinates (tree_vertices + 3, tree_vertices + 3, 1, translate_points_0_1); + } } /* Build a uniform forest */ t8_forest_t forest = t8_forest_new_uniform (cmesh, default_scheme, level, 1, sc_MPI_COMM_WORLD); const t8_locidx_t num_trees = t8_forest_get_num_local_trees (forest); - for (t8_locidx_t itree = 0; itree < num_trees; ++itree) { t8_log_indent_push (); const t8_locidx_t num_elements = t8_forest_get_tree_num_elements (forest, itree); @@ -211,7 +215,7 @@ TEST_P (geometry_point_inside, test_point_inside) const int num_corners = eclass_scheme->t8_element_num_corners (element); /* For each corner get its coordinates */ for (int icorner = 0; icorner < num_corners; ++icorner) { - t8_forest_element_coordinate (forest, 0, element, icorner, element_vertices[icorner]); + t8_forest_element_coordinate (forest, itree, element, icorner, element_vertices[icorner]); } /* Allocate the barycentric coordinates */ @@ -249,7 +253,7 @@ TEST_P (geometry_point_inside, test_point_inside) } /* Corrected number of points due to possible rounding errors in pow */ const int num_points = sc_intpow (num_steps, num_corners - 1); - double *test_point = T8_ALLOC (double, num_points * 3); + double *test_point = T8_ALLOC_ZERO (double, num_points * 3); int *point_is_inside = T8_ALLOC (int, num_points); int *point_is_recognized_as_inside = T8_ALLOC (int, num_points); double step = (barycentric_range_upper_bound - barycentric_range_lower_bound) / (num_steps - 1); @@ -270,7 +274,6 @@ TEST_P (geometry_point_inside, test_point_inside) } for (int icorner = 0; icorner < num_corners - 1; ++icorner) { int this_step = (ipoint / sc_intpow (num_steps, icorner)) % num_steps; - /* Set barycentric coordinates */ barycentric_coordinates[icorner] = barycentric_range_lower_bound + this_step * step; dampening = (1 - Sum) * (1 - Sum); @@ -285,7 +288,6 @@ TEST_P (geometry_point_inside, test_point_inside) /* The point is inside if and only if all barycentric coordinates are >= 0. */ point_is_inside[ipoint] = point_is_inside[ipoint] && barycentric_coordinates[icorner] >= 0; } - /* Ensure that sum over all bar. coordinates is 1 */ barycentric_coordinates[num_corners - 1] = 1 - Sum; @@ -294,6 +296,7 @@ TEST_P (geometry_point_inside, test_point_inside) test_point[ipoint * 3 + icoord] += barycentric_coordinates[num_corners - 1] * element_vertices[num_corners - 1][icoord]; } + /* The point is inside if and only if all barycentric coordinates are >= 0. */ point_is_inside[ipoint] = point_is_inside[ipoint] && barycentric_coordinates[num_corners - 1] >= 0; @@ -301,12 +304,18 @@ 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, - use_axis_aligned_geom, tolerance); + if (eclass == T8_ECLASS_LINE || eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_HEX) { + t8_forest_element_point_batch_inside (forest, 0, element, test_point, num_points, point_is_recognized_as_inside, + use_axis_aligned_geom, tolerance); + } + else { + t8_forest_element_point_batch_inside (forest, 0, element, test_point, num_points, point_is_recognized_as_inside, + 0, 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] - << ") should " << (point_is_inside[ipoint] ? "" : "not") << "be inside the " << t8_eclass_to_string[eclass] + << ") should " << (point_is_inside[ipoint] ? "" : "not ") << "be inside the " << t8_eclass_to_string[eclass] << " element, but is not detected as such."; } /* End loop over points. */ t8_debugf ("%i (%.2f%%) test points are inside the element\n", num_in, (100.0 * num_in) / num_points); @@ -323,9 +332,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::Range (0,1))); + testing::Range (0, 2))); + #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::Range (0, 1))); + testing::Range (0, 2))); #endif From 84a2a8c2ca91b1d5bdf700559366dda06935caba Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 15 Sep 2023 15:04:13 +0200 Subject: [PATCH 05/23] Adapt commit-check and forest_point_inside-check to axis-aligned-geom --- src/t8_cmesh/t8_cmesh.c | 10 +++++++++- src/t8_forest/t8_forest_cxx.cxx | 3 ++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/t8_cmesh/t8_cmesh.c b/src/t8_cmesh/t8_cmesh.c index fd60612ed9..1c21568885 100644 --- a/src/t8_cmesh/t8_cmesh.c +++ b/src/t8_cmesh/t8_cmesh.c @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -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_geom_is_linear_axis_aligned (t8_cmesh_get_tree_geometry (cmesh, gtree_id))) { + /* 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); } diff --git a/src/t8_forest/t8_forest_cxx.cxx b/src/t8_forest/t8_forest_cxx.cxx index 45c4f38166..8907fc0bb5 100644 --- a/src/t8_forest/t8_forest_cxx.cxx +++ b/src/t8_forest/t8_forest_cxx.cxx @@ -38,6 +38,7 @@ #include #if T8_ENABLE_DEBUG #include +#include #endif /* We want to export the whole implementation to be callable from "C" */ @@ -1241,7 +1242,7 @@ t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, c 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); const t8_geometry_c *geometry = t8_cmesh_get_tree_geometry (cmesh, cgtreeid); - T8_ASSERT (t8_geom_is_linear (geometry)); + T8_ASSERT (t8_geom_is_linear (geometry) || t8_geom_is_linear_axis_aligned (geometry)); #endif if (geom_is_axis_aligned) { double v_min[3]; From 1f0799a9280bc8d7f9a88782d3db2497224e267a Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 15 Sep 2023 15:04:40 +0200 Subject: [PATCH 06/23] Use tolerance in check --- src/t8_forest/t8_forest_cxx.cxx | 10 ++++++---- test/t8_geometry/t8_gtest_point_inside.cxx | 2 +- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/t8_forest/t8_forest_cxx.cxx b/src/t8_forest/t8_forest_cxx.cxx index 8907fc0bb5..b6d8a95961 100644 --- a/src/t8_forest/t8_forest_cxx.cxx +++ b/src/t8_forest/t8_forest_cxx.cxx @@ -1254,10 +1254,12 @@ t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, c 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*/ + is_inside[ipoint] = v_min[0] - tolerance <= points[ipoint * 3] && points[ipoint * 3] <= v_max[0] + tolerance + && /* check x-coordinate*/ + v_min[1] - tolerance <= points[ipoint * 3 + 1] + && points[ipoint * 3 + 1] <= v_max[1] + tolerance && /* check y-coordinate*/ + v_min[2] - tolerance <= points[ipoint * 3 + 2] + && points[ipoint * 3 + 2] <= v_max[2] + tolerance; /* check z-coordinate*/ } return; } diff --git a/test/t8_geometry/t8_gtest_point_inside.cxx b/test/t8_geometry/t8_gtest_point_inside.cxx index 44d175d6e4..96e8218d65 100644 --- a/test/t8_geometry/t8_gtest_point_inside.cxx +++ b/test/t8_geometry/t8_gtest_point_inside.cxx @@ -185,7 +185,7 @@ TEST_P (geometry_point_inside, test_point_inside) ? 2 : t8_eclass_num_vertices[eclass]; /* Translate all points by the same vector to move the element a bit. */ - double translate_all_points[3] = { 1, 0, 0 }; + double translate_all_points[3] = { -0.1, 0.3, 0.15 }; t8_cmesh_translate_coordinates (tree_vertices, tree_vertices, num_vertices, translate_all_points); /* Translate points 0 and 1 (if it exists) extra in order to move the 2D elements * and 3D faces outside of axis perpendicular planes. */ From 672217000e6e473df4051b1a2657d929ac6b398b Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 15 Sep 2023 15:41:28 +0200 Subject: [PATCH 07/23] Better formatting --- src/t8_forest/t8_forest_cxx.cxx | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/t8_forest/t8_forest_cxx.cxx b/src/t8_forest/t8_forest_cxx.cxx index b6d8a95961..a489ff352b 100644 --- a/src/t8_forest/t8_forest_cxx.cxx +++ b/src/t8_forest/t8_forest_cxx.cxx @@ -1254,12 +1254,13 @@ t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, c 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] - tolerance <= points[ipoint * 3] && points[ipoint * 3] <= v_max[0] + tolerance - && /* check x-coordinate*/ - v_min[1] - tolerance <= points[ipoint * 3 + 1] - && points[ipoint * 3 + 1] <= v_max[1] + tolerance && /* check y-coordinate*/ - v_min[2] - tolerance <= points[ipoint * 3 + 2] - && points[ipoint * 3 + 2] <= v_max[2] + tolerance; /* check z-coordinate*/ + /* check x-coordinate */ + /* check y-coordinate */ + /* check z-coordinate */ + is_inside[ipoint] + = v_min[0] - tolerance <= points[ipoint * 3] && points[ipoint * 3] <= v_max[0] + tolerance + && v_min[1] - tolerance <= points[ipoint * 3 + 1] && points[ipoint * 3 + 1] <= v_max[1] + tolerance + && v_min[2] - tolerance <= points[ipoint * 3 + 2] && points[ipoint * 3 + 2] <= v_max[2] + tolerance; } return; } From 4a711ed78989654b6527a11e53c54e77e4956d20 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 22 Sep 2023 15:15:18 +0200 Subject: [PATCH 08/23] Add point_inside_check into geometry class --- example/geometry/t8_example_geometries.cxx | 32 ++ src/t8_geometry/t8_geometry_base.hxx | 45 +++ .../t8_geometry_analytic.hxx | 17 + .../t8_geometry_examples.hxx | 8 + .../t8_geometry_linear.cxx | 351 ++++++++++++++++++ .../t8_geometry_linear.hxx | 14 + .../t8_geometry_linear_axis_aligned.cxx | 26 ++ .../t8_geometry_linear_axis_aligned.hxx | 14 + .../t8_geometry_occ.hxx | 17 + .../t8_geometry_zero.hxx | 17 + src/t8_geometry/t8_geometry_with_vertices.hxx | 18 + 11 files changed, 559 insertions(+) diff --git a/example/geometry/t8_example_geometries.cxx b/example/geometry/t8_example_geometries.cxx index c661d0c685..64eb84c91e 100644 --- a/example/geometry/t8_example_geometries.cxx +++ b/example/geometry/t8_example_geometries.cxx @@ -118,6 +118,14 @@ class t8_geometry_sincos: public t8_geometry { SC_ABORT_NOT_REACHED (); } + 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /* Load tree data is empty since we have no tree data. * We need to provide an implementation anyways. */ void @@ -214,6 +222,14 @@ class t8_geometry_cylinder: public t8_geometry { out_coords[2] = sin (ref_coords[0] * 2 * M_PI); } + 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /* Jacobian, not implemented. */ void t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, @@ -341,6 +357,14 @@ class t8_geometry_moving: public t8_geometry { out_coords[2] = 0; } + 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /* Jacobian, not implemented. */ void t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, @@ -390,6 +414,14 @@ class t8_geometry_cube_zdistorted: public t8_geometry { out_coords[2] = ref_coords[2] * (0.8 + 0.2 * sin (ref_coords[0] * 2 * M_PI) * cos (ref_coords[1] * 2 * M_PI)); } + 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /* Jacobian, not implemented. */ void t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, diff --git a/src/t8_geometry/t8_geometry_base.hxx b/src/t8_geometry/t8_geometry_base.hxx index 2d26e87df8..0b4b8945cd 100644 --- a/src/t8_geometry/t8_geometry_base.hxx +++ b/src/t8_geometry/t8_geometry_base.hxx @@ -30,6 +30,7 @@ #include #include +#include T8_EXTERN_C_BEGIN (); @@ -95,6 +96,50 @@ struct t8_geometry t8_geom_load_tree_data (t8_cmesh_t cmesh, t8_gloidx_t gtreeid) = 0; + /** Query whether a batch of points lies inside an element. + * \param [in] forest The forest. + * \param [in] ltree_id The forest local id of the tree in which the element is. + * \param [in] element The element. + * \param [in] points 3-dimensional coordinates of the points to check + * \param [in] num_points The number of points to check + * \param [in, out] is_inside An array of length \a num_points, filled with 0/1 on output. True (non-zero) if a \a point + * 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] 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. + */ + 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) + = 0; + + /** Query whether a single points lies inside an element. + * \param [in] forest The forest. + * \param [in] ltree_id The forest local id of the tree in which the element is. + * \param [in] element The element. + * \param [in] points 3-dimensional coordinates of the points to check + * \param [in] num_points The number of points to check + * \param [in, out] is_inside An array of length \a num_points, filled with 0/1 on output. True (non-zero) if a \a point + * 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] 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. + */ + inline int + t8_geom_point_inside_element (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, + const double *points, const int num_points, const double tolerance) + { + int is_inside = 0; + t8_geom_point_batch_inside_element (forest, ltreeid, element, points, 1, &is_inside, tolerance); + return is_inside; + } /** * Get the dimension of this geometry. * \return The dimension. diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx index 03a98b183f..623ea7bb4f 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx @@ -119,6 +119,23 @@ struct t8_geometry_analytic: public t8_geometry t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, double *jacobian) const; + /** + * \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 wether 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /** Update a possible internal data buffer for per tree data. * This function is called before the first coordinates in a new tree are * evaluated. You can use it for example to load the vertex coordinates of the diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx index 11ee18423f..f3893d4130 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx @@ -74,6 +74,14 @@ class t8_geometry_squared_disk: public t8_geometry { SC_ABORT_NOT_REACHED (); } + 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /* Load tree data is empty since we have no tree data. * We need to provide an implementation anyways. */ void diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx index a0ba7ea8f4..2625abe028 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx @@ -23,6 +23,8 @@ #include #include #include +#include +#include t8_geometry_linear::t8_geometry_linear (int dim): t8_geometry_with_vertices (dim, "") { @@ -54,6 +56,355 @@ t8_geometry_linear::t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtr SC_ABORT ("Not implemented."); } +/** 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 */ + 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; +} + +#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_vec_axpyz (p_0, p_1, A, -1); + + /* B = p2 - p0 */ + double B[3]; + t8_vec_axpyz (p_0, p_2, B, -1); + + /* C = p3 - p0 */ + double C[3]; + t8_vec_axpyz (p_0, p_3, C, -1); + + /* n1 = A x B */ + double A_cross_B[3]; + t8_vec_cross (A, B, A_cross_B); + + /* n2 = A x C */ + double A_cross_C[3]; + t8_vec_cross (A, C, A_cross_C); + + /* n1 x n2 */ + double n1_cross_n2[3]; + t8_vec_cross (A_cross_B, A_cross_C, n1_cross_n2); + + /* || n1 x n2 || */ + const double norm = t8_vec_norm (n1_cross_n2); + return norm < tolerance; +} +#endif + +void +t8_geometry_linear::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 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); + 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 + * >= 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 (); + } +} + T8_EXTERN_C_BEGIN (); /* Satisfy the C interface from t8_geometry_linear.h. diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx index dacea3abb0..6ba0f51634 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx @@ -80,6 +80,20 @@ struct t8_geometry_linear: public t8_geometry_with_vertices t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, double *jacobian) const; + /** + * \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 wether 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); + /* Load tree data is inherited from t8_geometry_with_vertices. */ }; diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx index f11615e350..7f5c733808 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx @@ -56,6 +56,32 @@ t8_geometry_linear_axis_aligned::t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8 SC_ABORT ("Not implemented."); } +void +t8_geometry_linear_axis_aligned::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) +{ + 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 */ + /* check x-coordinate */ + /* check y-coordinate */ + /* check z-coordinate */ + is_inside[ipoint] + = v_min[0] - tolerance <= points[ipoint * 3] && points[ipoint * 3] <= v_max[0] + tolerance + && v_min[1] - tolerance <= points[ipoint * 3 + 1] && points[ipoint * 3 + 1] <= v_max[1] + tolerance + && v_min[2] - tolerance <= points[ipoint * 3 + 2] && points[ipoint * 3 + 2] <= v_max[2] + tolerance; + } + return; +} + T8_EXTERN_C_BEGIN (); /* Satisfy the C interface from t8_geometry_linear_axis_aligned.h. diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx index b9ff4df069..629e869870 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx @@ -82,6 +82,20 @@ struct t8_geometry_linear_axis_aligned: 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; + + /** + * \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 wether 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); }; #endif /* !T8_GEOMETRY_LINEAR_AXIS_ALIGNED_HXX! */ diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx index 8efe5e15bd..62f0253c16 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx @@ -108,6 +108,23 @@ struct t8_geometry_occ: public t8_geometry_with_vertices t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, double *jacobian) const; + /** + * \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 wether 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /** Update a possible internal data buffer for per tree data. * This function is called before the first coordinates in a new tree are * evaluated. You can use it for example to load the vertex coordinates of the diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx index d09f1082bf..d632968c45 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx @@ -76,6 +76,23 @@ struct t8_geometry_zero: public t8_geometry t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, double *jacobian) const; + /** + * \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 wether 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /** Update a possible internal data buffer for per tree data. * This function is called before the first coordinates in a new tree are * evaluated. diff --git a/src/t8_geometry/t8_geometry_with_vertices.hxx b/src/t8_geometry/t8_geometry_with_vertices.hxx index b24b88a95a..aa13b0c9e4 100644 --- a/src/t8_geometry/t8_geometry_with_vertices.hxx +++ b/src/t8_geometry/t8_geometry_with_vertices.hxx @@ -29,6 +29,7 @@ #define T8_GEOMETRY_WITH_VERTICES_HXX #include +#include #include #include #include @@ -63,6 +64,23 @@ class t8_geometry_with_vertices: public t8_geometry { { } + /** + * \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 wether 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) + { + SC_ABORTF ("Function not yet implemented"); + } + /** Update a possible internal data buffer for per tree data. * This function is called before the first coordinates in a new tree are * evaluated. From e6b7db53579e93dff9ee1a8e657a4603b0480665 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 22 Sep 2023 15:26:45 +0200 Subject: [PATCH 09/23] remove const from t8_geom_handler_getter, as we want to use it --- src/t8_cmesh.h | 2 +- src/t8_cmesh/t8_cmesh_geometry.cxx | 2 +- src/t8_geometry/t8_geometry.cxx | 4 ++-- src/t8_geometry/t8_geometry.h | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/t8_cmesh.h b/src/t8_cmesh.h index b086521423..273fef21b1 100644 --- a/src/t8_cmesh.h +++ b/src/t8_cmesh.h @@ -464,7 +464,7 @@ t8_gloidx_t t8_cmesh_get_first_treeid (t8_cmesh_t cmesh); /* TODO: Comment */ -const t8_geometry_c * +t8_geometry_c * t8_cmesh_get_tree_geometry (t8_cmesh_t cmesh, t8_gloidx_t gtreeid); /** Query whether a given t8_locidx_t belongs to a local tree of a cmesh. diff --git a/src/t8_cmesh/t8_cmesh_geometry.cxx b/src/t8_cmesh/t8_cmesh_geometry.cxx index eef53b022e..34b43d9d81 100644 --- a/src/t8_cmesh/t8_cmesh_geometry.cxx +++ b/src/t8_cmesh/t8_cmesh_geometry.cxx @@ -53,7 +53,7 @@ t8_cmesh_set_tree_geometry (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const char *g t8_cmesh_set_attribute_string (cmesh, gtreeid, t8_get_package_id (), T8_CMESH_GEOMETRY_ATTRIBUTE_KEY, geom_name); } -const t8_geometry_c * +t8_geometry_c * t8_cmesh_get_tree_geometry (t8_cmesh_t cmesh, t8_gloidx_t gtreeid) { T8_ASSERT (t8_cmesh_is_committed (cmesh)); diff --git a/src/t8_geometry/t8_geometry.cxx b/src/t8_geometry/t8_geometry.cxx index 32f12ac933..18a6666314 100644 --- a/src/t8_geometry/t8_geometry.cxx +++ b/src/t8_geometry/t8_geometry.cxx @@ -253,14 +253,14 @@ t8_geom_handler_get_num_geometries (const t8_geometry_handler_t *geom_handler) return geom_handler->registered_geometries.elem_count; } -const t8_geometry_c * +t8_geometry_c * t8_geom_handler_get_unique_geometry (const t8_geometry_handler_t *geom_handler) { T8_ASSERT (t8_geom_handler_is_committed (geom_handler)); T8_ASSERT (t8_geom_handler_get_num_geometries (geom_handler) == 1); sc_array *geometries = (sc_array *) &geom_handler->registered_geometries; - return *(const t8_geometry_c **) sc_array_index_int (geometries, 0); + return *(t8_geometry_c **) sc_array_index_int (geometries, 0); } static inline void diff --git a/src/t8_geometry/t8_geometry.h b/src/t8_geometry/t8_geometry.h index 1bcb782a8f..2d5dc81771 100644 --- a/src/t8_geometry/t8_geometry.h +++ b/src/t8_geometry/t8_geometry.h @@ -125,7 +125,7 @@ t8_geom_handler_get_num_geometries (const t8_geometry_handler_t *geom_handler); * \note Most cmeshes will have only one geometry and this function is an optimization * for that special case. It is used for example in \ref t8_cmesh_get_tree_geometry. */ -const t8_geometry_c * +t8_geometry_c * t8_geom_handler_get_unique_geometry (const t8_geometry_handler_t *geom_handler); /** From fb58d513fbbdd2a690d2fc35bcaefbb80ead4f88 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 22 Sep 2023 15:27:17 +0200 Subject: [PATCH 10/23] typos --- .../t8_geometry_implementations/t8_geometry_analytic.hxx | 2 +- .../t8_geometry_implementations/t8_geometry_linear.hxx | 2 +- .../t8_geometry_linear_axis_aligned.hxx | 2 +- src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx | 2 +- .../t8_geometry_implementations/t8_geometry_zero.hxx | 2 +- src/t8_geometry/t8_geometry_with_vertices.hxx | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx index 623ea7bb4f..d6d96446f1 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_analytic.hxx @@ -125,7 +125,7 @@ struct t8_geometry_analytic: public t8_geometry * \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 wether the point is inside or not + * \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 diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx index 6ba0f51634..e3744c88c8 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx @@ -86,7 +86,7 @@ struct t8_geometry_linear: public t8_geometry_with_vertices * \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 wether the point is inside or not + * \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 diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx index 629e869870..a2b2b8746f 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx @@ -89,7 +89,7 @@ struct t8_geometry_linear_axis_aligned: public t8_geometry_with_vertices * \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 wether the point is inside or not + * \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 diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx index 62f0253c16..3d4d6d2eb7 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx @@ -114,7 +114,7 @@ struct t8_geometry_occ: public t8_geometry_with_vertices * \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 wether the point is inside or not + * \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 diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx index d632968c45..d0107ee6df 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx @@ -82,7 +82,7 @@ struct t8_geometry_zero: public t8_geometry * \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 wether the point is inside or not + * \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 diff --git a/src/t8_geometry/t8_geometry_with_vertices.hxx b/src/t8_geometry/t8_geometry_with_vertices.hxx index aa13b0c9e4..34caf87718 100644 --- a/src/t8_geometry/t8_geometry_with_vertices.hxx +++ b/src/t8_geometry/t8_geometry_with_vertices.hxx @@ -70,7 +70,7 @@ class t8_geometry_with_vertices: public t8_geometry { * \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 wether the point is inside or not + * \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 From b753503d3680b29393c6d876fbc13bf2d6d0cc8b Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 22 Sep 2023 15:28:24 +0200 Subject: [PATCH 11/23] Remove point_inside_check from forest --- src/t8_forest/t8_forest_cxx.cxx | 324 +------------------------------- 1 file changed, 2 insertions(+), 322 deletions(-) diff --git a/src/t8_forest/t8_forest_cxx.cxx b/src/t8_forest/t8_forest_cxx.cxx index a489ff352b..e01b4650a9 100644 --- a/src/t8_forest/t8_forest_cxx.cxx +++ b/src/t8_forest/t8_forest_cxx.cxx @@ -1066,337 +1066,17 @@ 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 */ - 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 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); - 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); - const t8_geometry_c *geometry = t8_cmesh_get_tree_geometry (cmesh, cgtreeid); - T8_ASSERT (t8_geom_is_linear (geometry) || t8_geom_is_linear_axis_aligned (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 */ - /* check x-coordinate */ - /* check y-coordinate */ - /* check z-coordinate */ - is_inside[ipoint] - = v_min[0] - tolerance <= points[ipoint * 3] && points[ipoint * 3] <= v_max[0] + tolerance - && v_min[1] - tolerance <= points[ipoint * 3 + 1] && points[ipoint * 3 + 1] <= v_max[1] + tolerance - && v_min[2] - tolerance <= points[ipoint * 3 + 2] && points[ipoint * 3 + 2] <= v_max[2] + tolerance; - } - return; - } - 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 - * >= 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 (); - } + 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); } int From c4633209cb386428dc11ae16b9f5f07f4cb8da41 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 6 Oct 2023 14:35:37 +0200 Subject: [PATCH 12/23] Use the geometry as an argument to pass to hypercube_pad --- example/remove/t8_example_empty_trees.cxx | 4 ++- src/t8_cmesh/t8_cmesh_examples.c | 25 +++++-------------- src/t8_cmesh/t8_cmesh_examples.h | 6 ++--- .../t8_gtest_cmesh_set_join_by_vertices.cxx | 4 ++- test/t8_geometry/t8_gtest_point_inside.cxx | 6 ++++- 5 files changed, 20 insertions(+), 25 deletions(-) diff --git a/example/remove/t8_example_empty_trees.cxx b/example/remove/t8_example_empty_trees.cxx index b7b7465e6d..ca03fe0b5f 100644 --- a/example/remove/t8_example_empty_trees.cxx +++ b/example/remove/t8_example_empty_trees.cxx @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -49,9 +50,10 @@ t8_strip_of_quads (t8_gloidx_t num_trees, t8_gloidx_t empty_tree, const char **v { const double boundary_coords[12] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0 }; + const t8_geometry_c *geometry = t8_geometry_linear_axis_aligned_new (2); t8_cmesh_t cmesh - = t8_cmesh_new_hypercube_pad (T8_ECLASS_QUAD, sc_MPI_COMM_WORLD, boundary_coords, num_trees, 1, 0, 1); + = t8_cmesh_new_hypercube_pad (T8_ECLASS_QUAD, sc_MPI_COMM_WORLD, boundary_coords, num_trees, 1, 0, geometry); t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default_cxx (), 0, 0, sc_MPI_COMM_WORLD); diff --git a/src/t8_cmesh/t8_cmesh_examples.c b/src/t8_cmesh/t8_cmesh_examples.c index 17278a5c2b..ed4b118115 100644 --- a/src/t8_cmesh/t8_cmesh_examples.c +++ b/src/t8_cmesh/t8_cmesh_examples.c @@ -1219,7 +1219,7 @@ t8_cmesh_set_vertices_3D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub t8_cmesh_t t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const double *boundary, t8_locidx_t polygons_x, - t8_locidx_t polygons_y, t8_locidx_t polygons_z, const int use_axis_aligned_geom) + t8_locidx_t polygons_y, t8_locidx_t polygons_z, const t8_geometry_c *geometry) { SC_CHECK_ABORT (eclass != T8_ECLASS_PYRAMID, "Pyramids are not yet supported."); const int dim = t8_eclass_to_dimension[eclass]; @@ -1240,22 +1240,9 @@ t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const do t8_cmesh_t cmesh; t8_cmesh_init (&cmesh); - if (use_axis_aligned_geom && (eclass == T8_ECLASS_HEX || eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_LINE)) { - /* We use axis aligned geometries */ - const t8_geometry_c *axis_aligned_geom = t8_geometry_linear_axis_aligned_new (dim); - t8_cmesh_register_geometry (cmesh, axis_aligned_geom); - } -#if T8_ENABLE_DEBUG - else if (use_axis_aligned_geom - && !(eclass == T8_ECLASS_HEX || eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_LINE)) { - SC_ABORTF ("Axis aligned geometry is not available for eclass %s!\n", t8_eclass_to_string[eclass]); - } -#endif - else { - /* We use standard linear geometry */ - const t8_geometry_c *linear_geom = t8_geometry_linear_new (dim); - t8_cmesh_register_geometry (cmesh, linear_geom); - } + const int is_axis_aligned = t8_geom_is_linear_axis_aligned (geometry); + + t8_cmesh_register_geometry (cmesh, geometry); /* Number of trees inside each polygon of given eclass. */ const t8_locidx_t num_trees_for_single_hypercube[T8_ECLASS_COUNT] = { 1, 1, 1, 2, 1, 6, 2, -1 }; @@ -1269,11 +1256,11 @@ t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const do /* Set the vertices of all trees. */ if (dim == 3) { T8_ASSERT (eclass == T8_ECLASS_HEX || eclass == T8_ECLASS_TET || eclass == T8_ECLASS_PRISM); - t8_cmesh_set_vertices_3D (cmesh, eclass, boundary, polygons_x, polygons_y, polygons_z, use_axis_aligned_geom); + t8_cmesh_set_vertices_3D (cmesh, eclass, boundary, polygons_x, polygons_y, polygons_z, is_axis_aligned); } else if (dim == 2) { T8_ASSERT (eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_TRIANGLE); - t8_cmesh_set_vertices_2D (cmesh, eclass, boundary, polygons_x, polygons_y, use_axis_aligned_geom); + t8_cmesh_set_vertices_2D (cmesh, eclass, boundary, polygons_x, polygons_y, is_axis_aligned); } else if (dim == 1) { T8_ASSERT (eclass == T8_ECLASS_LINE); diff --git a/src/t8_cmesh/t8_cmesh_examples.h b/src/t8_cmesh/t8_cmesh_examples.h index 0b0b423285..147ad29d9d 100644 --- a/src/t8_cmesh/t8_cmesh_examples.h +++ b/src/t8_cmesh/t8_cmesh_examples.h @@ -30,6 +30,7 @@ #include #include #include +#include T8_EXTERN_C_BEGIN (); @@ -106,8 +107,7 @@ t8_cmesh_new_hypercube (t8_eclass_t eclass, sc_MPI_Comm comm, int do_bcast, int * Only required if \a eclass is 2D or 3D. * \param [in] polygons_z The number of polygons along the z-axis. * Only required if \a eclass is 3D. - * \param [in] use_axis_aligned_geom Flag if cmesh uses the axis aligned_geometry. Only available for - * T8_ECLASS_LINE/QUAD/HEX + * \param [in] geometry The geometry to use. If the geometry is axis_aligned only two points per tree are stored * \return A committed t8_cmesh structure with * \a polygons_x * \a polygons_z * \a polygons_y many * sub-hypercubes of class \a eclass. @@ -134,7 +134,7 @@ t8_cmesh_new_hypercube (t8_eclass_t eclass, sc_MPI_Comm comm, int do_bcast, int */ t8_cmesh_t t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const double *boundary, t8_locidx_t polygons_x, - t8_locidx_t polygons_y, t8_locidx_t polygons_z, const int use_axis_aligned_geom); + t8_locidx_t polygons_y, t8_locidx_t polygons_z, const t8_geometry_c *geometry); /** Hybercube with 6 Tets, 6 Prism, 4 Hex. * \param [in] comm The mpi communicator to be used. diff --git a/test/t8_cmesh/t8_gtest_cmesh_set_join_by_vertices.cxx b/test/t8_cmesh/t8_gtest_cmesh_set_join_by_vertices.cxx index 12d4398f85..b78c5eb046 100644 --- a/test/t8_cmesh/t8_gtest_cmesh_set_join_by_vertices.cxx +++ b/test/t8_cmesh/t8_gtest_cmesh_set_join_by_vertices.cxx @@ -27,6 +27,7 @@ #include #include #include +#include #include @@ -241,7 +242,8 @@ TEST (t8_cmesh_set_join_by_vertices, test_cmesh_set_join_by_vertices) const double boundary_coords[24] = { 1, 0, 0, 4, 0, 0, 0, 6, 0, 5, 5, 0, -1, -2, 8, 9, 0, 10, 0, 8, 9, 10, 10, 10 }; t8_eclass_t eclass = T8_ECLASS_HEX; - t8_cmesh_t cmesh = t8_cmesh_new_hypercube_pad (eclass, comm, boundary_coords, 2, 2, 2, 0); + t8_geometry_c *geometry = new t8_geometry_linear (3); + t8_cmesh_t cmesh = t8_cmesh_new_hypercube_pad (eclass, comm, boundary_coords, 2, 2, 2, geometry); test_with_cmesh (cmesh); t8_cmesh_destroy (&cmesh); } diff --git a/test/t8_geometry/t8_gtest_point_inside.cxx b/test/t8_geometry/t8_gtest_point_inside.cxx index 96e8218d65..2996dd5873 100644 --- a/test/t8_geometry/t8_gtest_point_inside.cxx +++ b/test/t8_geometry/t8_gtest_point_inside.cxx @@ -33,6 +33,7 @@ #include #include #include +#include #include /* In this test we define a triangle in the x-y plane @@ -129,6 +130,7 @@ class geometry_point_inside: public testing::TestWithParam (GetParam ()); level = std::get<1> (GetParam ()); use_axis_aligned_geom = std::get<2> (GetParam ()); + const int dim = t8_eclass_to_dimension[eclass]; /* Construct a cube coarse mesh */ if (use_axis_aligned_geom && (eclass == T8_ECLASS_LINE || eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_HEX)) { @@ -144,7 +146,8 @@ class geometry_point_inside: public testing::TestWithParam Date: Thu, 18 Jan 2024 13:56:10 +0100 Subject: [PATCH 13/23] Update interface --- example/geometry/t8_example_geometries.cxx | 32 ------------------- src/t8_cmesh.h | 2 +- src/t8_cmesh/t8_cmesh_geometry.cxx | 2 +- src/t8_forest/t8_forest_cxx.cxx | 10 +++--- src/t8_forest/t8_forest_general.h | 9 ++---- src/t8_geometry/t8_geometry.cxx | 2 +- src/t8_geometry/t8_geometry.h | 2 +- src/t8_geometry/t8_geometry_base.hxx | 8 +++-- .../t8_geometry_examples.hxx | 8 ----- .../t8_geometry_linear.cxx | 3 +- .../t8_geometry_linear.hxx | 2 +- .../t8_geometry_linear_axis_aligned.cxx | 2 +- .../t8_geometry_linear_axis_aligned.hxx | 2 +- .../t8_geometry_occ.hxx | 17 ---------- src/t8_geometry/t8_geometry_with_vertices.hxx | 17 ---------- tutorials/general/t8_tutorial_search.cxx | 2 +- 16 files changed, 21 insertions(+), 99 deletions(-) diff --git a/example/geometry/t8_example_geometries.cxx b/example/geometry/t8_example_geometries.cxx index bd4ea32b84..1922b350ea 100644 --- a/example/geometry/t8_example_geometries.cxx +++ b/example/geometry/t8_example_geometries.cxx @@ -119,14 +119,6 @@ struct t8_geometry_sincos: public t8_geometry SC_ABORT_NOT_REACHED (); } - 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) - { - SC_ABORTF ("Function not yet implemented"); - } - /* Load tree data is empty since we have no tree data. * We need to provide an implementation anyways. */ void @@ -245,14 +237,6 @@ struct t8_geometry_cylinder: public t8_geometry out_coords[2] = sin (ref_coords[0] * 2 * M_PI); } - 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) - { - SC_ABORTF ("Function not yet implemented"); - } - /* Jacobian, not implemented. */ void t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, @@ -402,14 +386,6 @@ struct t8_geometry_moving: public t8_geometry out_coords[2] = 0; } - 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) - { - SC_ABORTF ("Function not yet implemented"); - } - /* Jacobian, not implemented. */ void t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, @@ -470,14 +446,6 @@ struct t8_geometry_cube_zdistorted: public t8_geometry out_coords[2] = ref_coords[2] * (0.8 + 0.2 * sin (ref_coords[0] * 2 * M_PI) * cos (ref_coords[1] * 2 * M_PI)); } - 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) - { - SC_ABORTF ("Function not yet implemented"); - } - /* Jacobian, not implemented. */ void t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, diff --git a/src/t8_cmesh.h b/src/t8_cmesh.h index 7ac50358ab..206bb7ff31 100644 --- a/src/t8_cmesh.h +++ b/src/t8_cmesh.h @@ -460,7 +460,7 @@ t8_gloidx_t t8_cmesh_get_first_treeid (t8_cmesh_t cmesh); /* TODO: Comment */ -t8_geometry_c * +const t8_geometry_c * t8_cmesh_get_tree_geometry (t8_cmesh_t cmesh, t8_gloidx_t gtreeid); /** Query whether a given t8_locidx_t belongs to a local tree of a cmesh. diff --git a/src/t8_cmesh/t8_cmesh_geometry.cxx b/src/t8_cmesh/t8_cmesh_geometry.cxx index 34b43d9d81..eef53b022e 100644 --- a/src/t8_cmesh/t8_cmesh_geometry.cxx +++ b/src/t8_cmesh/t8_cmesh_geometry.cxx @@ -53,7 +53,7 @@ t8_cmesh_set_tree_geometry (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const char *g t8_cmesh_set_attribute_string (cmesh, gtreeid, t8_get_package_id (), T8_CMESH_GEOMETRY_ATTRIBUTE_KEY, geom_name); } -t8_geometry_c * +const t8_geometry_c * t8_cmesh_get_tree_geometry (t8_cmesh_t cmesh, t8_gloidx_t gtreeid) { T8_ASSERT (t8_cmesh_is_committed (cmesh)); diff --git a/src/t8_forest/t8_forest_cxx.cxx b/src/t8_forest/t8_forest_cxx.cxx index f6d341410e..4c74d2f3c6 100644 --- a/src/t8_forest/t8_forest_cxx.cxx +++ b/src/t8_forest/t8_forest_cxx.cxx @@ -1091,24 +1091,22 @@ t8_forest_element_face_normal (t8_forest_t forest, t8_locidx_t ltreeid, const t8 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 int geom_is_axis_aligned, const double tolerance) + const double *points, int num_points, int *is_inside, const double tolerance) { /* 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_geometry_c *geometry = t8_cmesh_get_tree_geometry (cmesh, cgtreeid); + 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); } int t8_forest_element_point_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, - const double point[3], const int geom_is_axis_aligned, const double tolerance) + const double point[3], const double tolerance) { int is_inside = 0; - t8_forest_element_point_batch_inside (forest, ltreeid, element, point, 1, &is_inside, geom_is_axis_aligned, - tolerance); + t8_forest_element_point_batch_inside (forest, ltreeid, element, point, 1, &is_inside, tolerance); return is_inside; } diff --git a/src/t8_forest/t8_forest_general.h b/src/t8_forest/t8_forest_general.h index 04a3ab87fa..8b06def117 100644 --- a/src/t8_forest/t8_forest_general.h +++ b/src/t8_forest/t8_forest_general.h @@ -752,8 +752,6 @@ 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 @@ -761,7 +759,7 @@ t8_forest_iterate (t8_forest_t forest); */ int t8_forest_element_point_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, - const double point[3], const int geom_is_axis_aligned, const double tolerance); + const double point[3], 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 @@ -776,8 +774,6 @@ 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 @@ -785,8 +781,7 @@ t8_forest_element_point_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t */ 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 int geom_is_axis_aligned, const double tolerance); + const double *points, int num_points, int *is_inside, const double tolerance); /* TODO: if set level and partition/adapt/balance all give NULL, then * refine uniformly and partition/adapt/balance the unfiform forest. */ diff --git a/src/t8_geometry/t8_geometry.cxx b/src/t8_geometry/t8_geometry.cxx index 4fa8b761e4..8ba392d219 100644 --- a/src/t8_geometry/t8_geometry.cxx +++ b/src/t8_geometry/t8_geometry.cxx @@ -253,7 +253,7 @@ t8_geom_handler_get_num_geometries (const t8_geometry_handler_t *geom_handler) return geom_handler->registered_geometries.elem_count; } -t8_geometry_c * +const t8_geometry_c * t8_geom_handler_get_unique_geometry (const t8_geometry_handler_t *geom_handler) { T8_ASSERT (t8_geom_handler_is_committed (geom_handler)); diff --git a/src/t8_geometry/t8_geometry.h b/src/t8_geometry/t8_geometry.h index cd75c87299..9394fedf03 100644 --- a/src/t8_geometry/t8_geometry.h +++ b/src/t8_geometry/t8_geometry.h @@ -143,7 +143,7 @@ t8_geom_handler_get_num_geometries (const t8_geometry_handler_t *geom_handler); * \note Most cmeshes will have only one geometry and this function is an optimization * for that special case. It is used for example in \ref t8_cmesh_get_tree_geometry. */ -t8_geometry_c * +const t8_geometry_c * t8_geom_handler_get_unique_geometry (const t8_geometry_handler_t *geom_handler); /** diff --git a/src/t8_geometry/t8_geometry_base.hxx b/src/t8_geometry/t8_geometry_base.hxx index 24d0b927d3..e95564af11 100644 --- a/src/t8_geometry/t8_geometry_base.hxx +++ b/src/t8_geometry/t8_geometry_base.hxx @@ -115,8 +115,10 @@ struct t8_geometry 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) - = 0; + const double tolerance) const + { + SC_ABORTF ("Function not yet implemented"); + }; /** Query whether a single points lies inside an element. * \param [in] forest The forest. @@ -135,7 +137,7 @@ struct t8_geometry */ inline int t8_geom_point_inside_element (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, - const double *points, const int num_points, const double tolerance) + const double *points, const int num_points, const double tolerance) const { int is_inside = 0; t8_geom_point_batch_inside_element (forest, ltreeid, element, points, 1, &is_inside, tolerance); diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx index 3492a76748..702c9acefc 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx @@ -75,14 +75,6 @@ struct t8_geometry_squared_disk: public t8_geometry_with_vertices SC_ABORT_NOT_REACHED (); } - 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) - { - SC_ABORTF ("Function not yet implemented"); - } - /* Load tree data is empty since we have no tree data. * We need to provide an implementation anyways. */ void diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx index 9223853226..1bb303238e 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx @@ -265,7 +265,8 @@ t8_four_points_coplanar (const double p_0[3], const double p_1[3], const double void t8_geometry_linear::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 int num_points, int *is_inside, + const double tolerance) const { 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); diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx index b3ca4bd449..6840b9dc64 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.hxx @@ -102,7 +102,7 @@ struct t8_geometry_linear: public t8_geometry_with_vertices 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 double tolerance) const; /* Load tree data is inherited from t8_geometry_with_vertices. */ }; diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx index e0c1bcdbe0..359aa16123 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx @@ -60,7 +60,7 @@ void t8_geometry_linear_axis_aligned::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 double tolerance) const { double v_min[3]; double v_max[3]; diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx index dc3bddc21d..339e8b446a 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.hxx @@ -105,7 +105,7 @@ struct t8_geometry_linear_axis_aligned: public t8_geometry_with_vertices 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 double tolerance) const; }; #endif /* !T8_GEOMETRY_LINEAR_AXIS_ALIGNED_HXX! */ diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx index 92233a65ea..077d43c632 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_occ.hxx @@ -118,23 +118,6 @@ struct t8_geometry_occ: public t8_geometry_with_vertices t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtreeid, const double *ref_coords, const size_t num_coords, double *jacobian) const; - /** - * \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) - { - SC_ABORTF ("Function not yet implemented"); - } - /** Update a possible internal data buffer for per tree data. * This function is called before the first coordinates in a new tree are * evaluated. You can use it for example to load the vertex coordinates of the diff --git a/src/t8_geometry/t8_geometry_with_vertices.hxx b/src/t8_geometry/t8_geometry_with_vertices.hxx index f99808a8ff..8dd2d8b1a4 100644 --- a/src/t8_geometry/t8_geometry_with_vertices.hxx +++ b/src/t8_geometry/t8_geometry_with_vertices.hxx @@ -65,23 +65,6 @@ struct t8_geometry_with_vertices: public t8_geometry { } - /** - * \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) - { - SC_ABORTF ("Function not yet implemented"); - } - /** Update a possible internal data buffer for per tree data. * This function is called before the first coordinates in a new tree are * evaluated. diff --git a/tutorials/general/t8_tutorial_search.cxx b/tutorials/general/t8_tutorial_search.cxx index 3dee79a2e6..068cb0038e 100644 --- a/tutorials/general/t8_tutorial_search.cxx +++ b/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, 1, tolerance); + = t8_forest_element_point_inside (forest, ltreeid, element, particle->coordinates, tolerance); if (particle_is_inside_element) { if (is_leaf) { /* The particle is inside and this element is a leaf element. From 13ee9175082590a46e5f6af43c94bee34ed0922d Mon Sep 17 00:00:00 2001 From: David Knapp Date: Thu, 18 Jan 2024 14:02:58 +0100 Subject: [PATCH 14/23] remove t8_forest_point_inside renamed batched function to t8_forest_points_inside this function can be used for a single point too --- src/t8_forest/t8_forest_cxx.cxx | 13 ++---------- src/t8_forest/t8_forest_general.h | 25 ++---------------------- tutorials/general/t8_tutorial_search.cxx | 6 +++--- 3 files changed, 7 insertions(+), 37 deletions(-) diff --git a/src/t8_forest/t8_forest_cxx.cxx b/src/t8_forest/t8_forest_cxx.cxx index 4c74d2f3c6..419799c5d3 100644 --- a/src/t8_forest/t8_forest_cxx.cxx +++ b/src/t8_forest/t8_forest_cxx.cxx @@ -1090,8 +1090,8 @@ t8_forest_element_face_normal (t8_forest_t forest, t8_locidx_t ltreeid, const t8 } 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) { /* Check whether the provided geometry is linear */ const t8_cmesh_t cmesh = t8_forest_get_cmesh (forest); @@ -1101,15 +1101,6 @@ t8_forest_element_point_batch_inside (t8_forest_t forest, t8_locidx_t ltreeid, c geometry->t8_geom_point_batch_inside_element (forest, ltreeid, element, points, num_points, is_inside, tolerance); } -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; -} - /* For each tree in a forest compute its first and last descendant */ void t8_forest_compute_desc (t8_forest_t forest) diff --git a/src/t8_forest/t8_forest_general.h b/src/t8_forest/t8_forest_general.h index 8b06def117..de2aada535 100644 --- a/src/t8_forest/t8_forest_general.h +++ b/src/t8_forest/t8_forest_general.h @@ -740,27 +740,6 @@ t8_forest_element_face_neighbor (t8_forest_t forest, t8_locidx_t ltreeid, const void t8_forest_iterate (t8_forest_t forest); -/** Query whether a given point lies inside an element or not. For bilinearly interpolated elements. - * \note For 2D quadrilateral elements this function is only an approximation. It is correct - * if the four vertices lie in the same plane, but it may produce only approximate results if - * the vertices do not lie in the same plane. - * \param [in] forest The forest. - * \param [in] ltree_id The forest local id of the tree in which the element is. - * \param [in] element The element. - * \param [in] point 3-dimensional coordinates of the point to check - * \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. - * \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); - /** 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 * if the four vertices lie in the same plane, but it may produce only approximate results if @@ -780,8 +759,8 @@ t8_forest_element_point_inside (t8_forest_t forest, t8_locidx_t ltreeid, const t * 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); +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); /* TODO: if set level and partition/adapt/balance all give NULL, then * refine uniformly and partition/adapt/balance the unfiform forest. */ diff --git a/tutorials/general/t8_tutorial_search.cxx b/tutorials/general/t8_tutorial_search.cxx index 068cb0038e..8f1095263c 100644 --- a/tutorials/general/t8_tutorial_search.cxx +++ b/tutorials/general/t8_tutorial_search.cxx @@ -189,10 +189,10 @@ t8_tutorial_search_query_callback (t8_forest_t forest, t8_locidx_t ltreeid, cons /* Ensure that the data is actually set. */ T8_ASSERT (particles_per_element != NULL); T8_ASSERT (query != NULL); - + const int num_particles = 1; /* 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_points_inside (forest, ltreeid, element, particle->coordinates, num_particles, + &particle_is_inside_element, tolerance); if (particle_is_inside_element) { if (is_leaf) { /* The particle is inside and this element is a leaf element. From 868544906842ee609ff532cace6e9bc86fbec07e Mon Sep 17 00:00:00 2001 From: David Knapp Date: Thu, 18 Jan 2024 16:29:47 +0100 Subject: [PATCH 15/23] Update test --- test/t8_geometry/t8_gtest_point_inside.cxx | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/test/t8_geometry/t8_gtest_point_inside.cxx b/test/t8_geometry/t8_gtest_point_inside.cxx index 2996dd5873..720f43954d 100644 --- a/test/t8_geometry/t8_gtest_point_inside.cxx +++ b/test/t8_geometry/t8_gtest_point_inside.cxx @@ -72,7 +72,8 @@ 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, 0, tolerance); + int point_is_inside; + t8_forest_element_points_inside (forest, 0, element, test_point, 1, &point_is_inside, tolerance); ASSERT_FALSE (point_is_inside) << "The point is wrongly detected as inside the triangle."; t8_forest_unref (&forest); } @@ -114,7 +115,8 @@ 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, 0, tolerance); + int point_is_inside; + t8_forest_element_points_inside (forest, 0, element, test_point, 1, &point_is_inside, tolerance); ASSERT_FALSE (point_is_inside) << "The point is wrongly detected as inside the quad."; @@ -308,14 +310,8 @@ 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. */ - if (eclass == T8_ECLASS_LINE || eclass == T8_ECLASS_QUAD || eclass == T8_ECLASS_HEX) { - t8_forest_element_point_batch_inside (forest, 0, element, test_point, num_points, point_is_recognized_as_inside, - use_axis_aligned_geom, tolerance); - } - else { - t8_forest_element_point_batch_inside (forest, 0, element, test_point, num_points, point_is_recognized_as_inside, - 0, tolerance); - } + t8_forest_element_points_inside (forest, 0, element, test_point, num_points, point_is_recognized_as_inside, + 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] From c8ec1270ed1bd2e381e9d8e1021da2852ac0f2a8 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Mon, 29 Jan 2024 12:03:18 +0100 Subject: [PATCH 16/23] Fix cp error --- src/t8_cmesh/t8_cmesh_examples.c | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/t8_cmesh/t8_cmesh_examples.c b/src/t8_cmesh/t8_cmesh_examples.c index 64ea4da374..0c7919143f 100644 --- a/src/t8_cmesh/t8_cmesh_examples.c +++ b/src/t8_cmesh/t8_cmesh_examples.c @@ -987,9 +987,6 @@ t8_cmesh_set_vertices_2D (t8_cmesh_t cmesh, const t8_eclass_t eclass, const doub memcpy (vertices, box, 3 * sizeof (double)); /* Vertex 0 */ t8_vec_axpyz (box, box_dir, vertices + 6, 1.0); /* Vertex 2 */ - memcpy (vertices, box, 3 * sizeof (double)); /* Vertex 0 */ - t8_vec_axpyz (box, box_dir, vertices + 6, 1.0); /* Vertex 2 */ - /* Reduce box along x axis */ t8_resize_box (2, box, box_dir, 0, 1, box_quads); t8_update_box_face_edges (2, box, box_dir, 0, box_quads); From 67e126d5d694fa43a63044deb5b2a48fc8154f37 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Mon, 29 Jan 2024 12:54:09 +0100 Subject: [PATCH 17/23] Update geometry check --- src/t8_cmesh/t8_cmesh.c | 2 +- src/t8_cmesh/t8_cmesh_examples.c | 2 +- .../t8_geometry_linear_axis_aligned.cxx | 13 ------------- .../t8_geometry_linear_axis_aligned.h | 8 -------- 4 files changed, 2 insertions(+), 23 deletions(-) diff --git a/src/t8_cmesh/t8_cmesh.c b/src/t8_cmesh/t8_cmesh.c index 00dabe0837..1693262826 100644 --- a/src/t8_cmesh/t8_cmesh.c +++ b/src/t8_cmesh/t8_cmesh.c @@ -565,7 +565,7 @@ t8_cmesh_no_negative_volume (t8_cmesh_t cmesh) /* Vertices are set */ eclass = t8_cmesh_get_tree_class (cmesh, itree); const t8_gloidx_t gtree_id = t8_cmesh_get_global_id (cmesh, itree); - if (t8_geom_is_linear_axis_aligned (t8_cmesh_get_tree_geometry (cmesh, gtree_id))) { + 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]; } diff --git a/src/t8_cmesh/t8_cmesh_examples.c b/src/t8_cmesh/t8_cmesh_examples.c index 457a1dc801..3187a5a142 100644 --- a/src/t8_cmesh/t8_cmesh_examples.c +++ b/src/t8_cmesh/t8_cmesh_examples.c @@ -1279,7 +1279,7 @@ t8_cmesh_new_hypercube_pad (const t8_eclass_t eclass, sc_MPI_Comm comm, const do t8_cmesh_t cmesh; t8_cmesh_init (&cmesh); - const int is_axis_aligned = t8_geom_is_linear_axis_aligned (geometry); + const int is_axis_aligned = t8_geom_get_type (geometry) == T8_GEOMETRY_TYPE_LINEAR_AXIS_ALIGNED; t8_cmesh_register_geometry (cmesh, geometry); diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx index 359aa16123..fc23664338 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.cxx @@ -103,17 +103,4 @@ t8_geometry_linear_axis_aligned_destroy (t8_geometry_c **geom) *geom = NULL; } -int -t8_geom_is_linear_axis_aligned (const t8_geometry_c *geometry) -{ - /* Try to dynamic cast the geometry into linear, axis-aligned geometry. - * This is only successful if geometry pointed to a - * t8_geometry_linear_axis_aligned. - * If successful, then is_linear_geom will be true. - */ - const int is_linear_axis_aligned_geom = (dynamic_cast (geometry) != NULL); - - return is_linear_axis_aligned_geom; -} - T8_EXTERN_C_END (); diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h index 0dd0329390..b741d11ce2 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h @@ -54,14 +54,6 @@ t8_geometry_linear_axis_aligned_new (int dim); void t8_geometry_linear_axis_aligned_destroy (t8_geometry_c **geom); -/** Query whether a given geometry is \ref t8_geometry_linear_axis_aligned. - * \param [in] geometry A geometry. - * \return True (non-zero) if and only if the geometry is of - * type \ref t8_geometry_linear_axis_aligned. - */ -int -t8_geom_is_linear_axis_aligned (const t8_geometry_c *geometry); - T8_EXTERN_C_END (); #endif /* !T8_GEOMETRY_LINEAR_AXIS_ALIGNED_H! */ From 8ddd639f6afd889cdfd8a532acc7238d391cde89 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Fri, 9 Feb 2024 10:21:47 +0100 Subject: [PATCH 18/23] Clean-up --- .../t8_geometry_examples.hxx | 7 ------- test/Makefile.am | 10 ---------- 2 files changed, 17 deletions(-) diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx index d2bb8d1948..0f979e3e09 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_examples.hxx @@ -75,13 +75,6 @@ struct t8_geometry_squared_disk: public t8_geometry_with_vertices SC_ABORT_NOT_REACHED (); } - /* Load tree data is empty since we have no tree data. - * We need to provide an implementation anyways. */ - void - t8_geom_load_tree_data (t8_cmesh_t cmesh, t8_gloidx_t gtreeid) - { - SC_ABORT_NOT_REACHED (); - } /** * Get the type of this geometry. * \return The type. diff --git a/test/Makefile.am b/test/Makefile.am index 0bfd8c213f..b378b009fe 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -28,7 +28,6 @@ t8code_googletest_programs = \ test/t8_gtest_basics \ test/t8_schemes/t8_gtest_ancestor \ test/t8_cmesh/t8_gtest_hypercube \ - test/t8_cmesh/t8_gtest_cmesh_copy \ test/t8_schemes/t8_gtest_element_count_leaves \ test/t8_schemes/t8_gtest_element_ref_coords \ test/t8_geometry/t8_gtest_geometry \ @@ -134,10 +133,6 @@ test_t8_cmesh_t8_gtest_hypercube_SOURCES = \ test/t8_gtest_main.cxx \ test/t8_cmesh/t8_gtest_hypercube.cxx -test_t8_cmesh_t8_gtest_cmesh_copy_SOURCES = \ - test/t8_gtest_main.cxx \ - test/t8_cmesh/t8_gtest_cmesh_copy.cxx - test_t8_cmesh_t8_gtest_cmesh_set_join_by_vertices_SOURCES = \ test/t8_gtest_main.cxx \ test/t8_cmesh/t8_gtest_cmesh_set_join_by_vertices.cxx @@ -362,10 +357,6 @@ test_t8_cmesh_t8_gtest_hypercube_LDADD = $(t8_gtest_target_ld_add) test_t8_cmesh_t8_gtest_hypercube_LDFLAGS = $(t8_gtest_target_ld_flags) test_t8_cmesh_t8_gtest_hypercube_CPPFLAGS = $(t8_gtest_target_cpp_flags) -test_t8_cmesh_t8_gtest_cmesh_copy_LDADD = $(t8_gtest_target_ld_add) -test_t8_cmesh_t8_gtest_cmesh_copy_LDFLAGS = $(t8_gtest_target_ld_flags) -test_t8_cmesh_t8_gtest_cmesh_copy_CPPFLAGS = $(t8_gtest_target_cpp_flags) - test_t8_cmesh_t8_gtest_cmesh_set_join_by_vertices_LDADD = $(t8_gtest_target_ld_add) test_t8_cmesh_t8_gtest_cmesh_set_join_by_vertices_LDFLAGS = $(t8_gtest_target_ld_flags) test_t8_cmesh_t8_gtest_cmesh_set_join_by_vertices_CPPFLAGS = $(t8_gtest_target_cpp_flags) @@ -544,7 +535,6 @@ test_t8_schemes_t8_gtest_init_linear_id_CPPFLAGS += $(t8_gtest_target_mpi_cpp_fl test_t8_gtest_basics_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) test_t8_schemes_t8_gtest_ancestor_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) test_t8_cmesh_t8_gtest_hypercube_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) -test_t8_cmesh_t8_gtest_cmesh_copy_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) test_t8_cmesh_t8_gtest_cmesh_set_join_by_vertices_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) test_t8_schemes_t8_gtest_element_count_leaves_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) test_t8_schemes_t8_gtest_element_ref_coords_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) From dd0c61a109a94a21ba88a1f4973bd88e5f0b0f9d Mon Sep 17 00:00:00 2001 From: David Knapp Date: Tue, 13 Feb 2024 14:28:38 +0100 Subject: [PATCH 19/23] Reduce particle number to 2000 --- tutorials/general/t8_tutorial_search.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tutorials/general/t8_tutorial_search.cxx b/tutorials/general/t8_tutorial_search.cxx index 5e8837793f..07064e08f7 100644 --- a/tutorials/general/t8_tutorial_search.cxx +++ b/tutorials/general/t8_tutorial_search.cxx @@ -373,7 +373,7 @@ main (int argc, char **argv) /* The uniform refinement level of the forest. */ const int level = 5; /* The number of particles to generate. */ - const size_t num_particles = 200000; + const size_t num_particles = 2000; /* The seed for the random number generator. */ const unsigned seed = 0; From a78ae7f721f3588b2ebbaa67ffe62d4455c712ab Mon Sep 17 00:00:00 2001 From: David Knapp Date: Tue, 13 Feb 2024 14:38:46 +0100 Subject: [PATCH 20/23] Moved inside-check-functions into geometry_helper --- src/t8_geometry/t8_geometry_helpers.c | 133 +++++++++++++++ src/t8_geometry/t8_geometry_helpers.h | 47 +++++ .../t8_geometry_linear.cxx | 160 ------------------ 3 files changed, 180 insertions(+), 160 deletions(-) diff --git a/src/t8_geometry/t8_geometry_helpers.c b/src/t8_geometry/t8_geometry_helpers.c index 2a8738b5c3..d216d08bf4 100644 --- a/src/t8_geometry/t8_geometry_helpers.c +++ b/src/t8_geometry/t8_geometry_helpers.c @@ -373,3 +373,136 @@ t8_geom_get_triangle_scaling_factor (int edge_index, const double *tree_vertices double scaling_factor = dist_ref / dist_intersection; return scaling_factor; } + +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; +} + +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; +} + +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. + */ +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 */ + 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; +} \ No newline at end of file diff --git a/src/t8_geometry/t8_geometry_helpers.h b/src/t8_geometry/t8_geometry_helpers.h index 5b3eac2d27..498ed9b8b3 100644 --- a/src/t8_geometry/t8_geometry_helpers.h +++ b/src/t8_geometry/t8_geometry_helpers.h @@ -132,6 +132,53 @@ double t8_geom_get_triangle_scaling_factor (int edge_index, const double *tree_vertices, const double *glob_intersection, const double *glob_ref_point); +/** 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. + */ +int +t8_vertex_point_inside (const double vertex_coords[3], const double point[3], const double tolerance); + +/** + * 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. + */ +int +t8_line_point_inside (const double *p_0, const double *vec, const double *point, const double tolerance); + +/** + * 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. + */ +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); + +/** 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. + */ +int +t8_plane_point_inside (const double point_on_face[3], const double face_normal[3], const double point[3]); + T8_EXTERN_C_END (); #endif /* !T8_GEOMETRY_HELPERS_H! */ diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx index 1bb303238e..242d3716d0 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_linear.cxx @@ -56,166 +56,6 @@ t8_geometry_linear::t8_geom_evaluate_jacobian (t8_cmesh_t cmesh, t8_gloidx_t gtr SC_ABORT ("Not implemented."); } -/** 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 */ - 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; -} - #if T8_ENABLE_DEBUG /* Test whether four given points in 3D are coplanar up to a given tolerance. */ From c8459065ea175817932cee5262cefb79b281989d Mon Sep 17 00:00:00 2001 From: David Knapp Date: Tue, 13 Feb 2024 15:15:41 +0100 Subject: [PATCH 21/23] Update Makefile --- test/Makefile.am | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/test/Makefile.am b/test/Makefile.am index 114d8ffbd0..607f07e6a4 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -290,6 +290,10 @@ test_t8_cmesh_generator_t8_gtest_cmesh_generator_test_SOURCES = \ test/t8_gtest_main.cxx \ test/t8_cmesh_generator/t8_gtest_cmesh_generator_test.cxx +test_t8_cmesh_t8_gtest_cmesh_copy_SOURCES = \ + test/t8_gtest_main.cxx \ + test/t8_cmesh/t8_gtest_cmesh_copy.cxx + #define ld and cpp flags for all targets t8_gtest_target_ld_add = $(LDADD) test/libgtest.la t8_gtest_target_ld_flags = $(AM_LDFLAGS) -pthread @@ -518,6 +522,10 @@ test_t8_cmesh_generator_t8_gtest_cmesh_generator_test_LDADD = $(t8_gtest_target_ test_t8_cmesh_generator_t8_gtest_cmesh_generator_test_LDFLAGS = $(t8_gtest_target_ld_flags) test_t8_cmesh_generator_t8_gtest_cmesh_generator_test_CPPFLAGS = $(t8_gtest_target_cpp_flags) +test_t8_cmesh_t8_gtest_cmesh_copy_LDADD = $(t8_gtest_target_ld_add) +test_t8_cmesh_t8_gtest_cmesh_copy_LDFLAGS = $(t8_gtest_target_ld_flags) +test_t8_cmesh_t8_gtest_cmesh_copy_CPPFLAGS = $(t8_gtest_target_cpp_flags) + # If we did not configure t8code with MPI we need to build Googletest # without MPI support. if !T8_ENABLE_MPI @@ -576,6 +584,7 @@ test_t8_cmesh_t8_gtest_cmesh_tree_vertices_negative_volume_CPPFLAGS += $(t8_gtes test_t8_schemes_t8_gtest_default_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) test_t8_schemes_t8_gtest_child_parent_face_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) test_t8_cmesh_generator_t8_gtest_cmesh_generator_test_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) +test_t8_cmesh_t8_gtest_cmesh_copy_CPPFLAGS += $(t8_gtest_target_mpi_cpp_flags) endif From b6310f64e3e067f55da3b49233dc3a69a16e3d81 Mon Sep 17 00:00:00 2001 From: David Knapp Date: Wed, 14 Feb 2024 13:25:47 +0100 Subject: [PATCH 22/23] Update src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx Co-authored-by: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> --- .../t8_geometry_implementations/t8_geometry_zero.hxx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx index 1219ec0cad..e8cbce80bd 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx @@ -100,7 +100,12 @@ struct t8_geometry_zero: public t8_geometry const double *points, const int num_points, int *is_inside, const double tolerance) { - SC_ABORTF ("Function not yet implemented"); +const int zeros[3] = { 0 }; +for (int i_point = 0; i_point < num_points; ++i_point) +{ + const int offset = i_point * T8_ECLASS_MAX_DIM; + is_inside[i_point] = t8_vertex_point_inside(zeros, points + offset, tolerance) +} } /** Update a possible internal data buffer for per tree data. From 7d7fe2f90431515bc3bc9cd86025e08b243c32fb Mon Sep 17 00:00:00 2001 From: David Knapp Date: Wed, 14 Feb 2024 13:41:56 +0100 Subject: [PATCH 23/23] fix compile errors --- src/t8_geometry/t8_geometry_helpers.c | 7 ------- .../t8_geometry_implementations/t8_geometry_zero.hxx | 12 ++++++------ 2 files changed, 6 insertions(+), 13 deletions(-) diff --git a/src/t8_geometry/t8_geometry_helpers.c b/src/t8_geometry/t8_geometry_helpers.c index d216d08bf4..0abef3326e 100644 --- a/src/t8_geometry/t8_geometry_helpers.c +++ b/src/t8_geometry/t8_geometry_helpers.c @@ -485,13 +485,6 @@ t8_triangle_point_inside (const double p_0[3], const double v[3], const double w 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. - */ int t8_plane_point_inside (const double point_on_face[3], const double face_normal[3], const double point[3]) { diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx index e8cbce80bd..cfcd99cbba 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_zero.hxx @@ -31,6 +31,7 @@ #include #include +#include struct t8_geometry_zero: public t8_geometry { @@ -100,12 +101,11 @@ struct t8_geometry_zero: public t8_geometry const double *points, const int num_points, int *is_inside, const double tolerance) { -const int zeros[3] = { 0 }; -for (int i_point = 0; i_point < num_points; ++i_point) -{ - const int offset = i_point * T8_ECLASS_MAX_DIM; - is_inside[i_point] = t8_vertex_point_inside(zeros, points + offset, tolerance) -} + const double zeros[3] = { 0 }; + for (int i_point = 0; i_point < num_points; ++i_point) { + const int offset = i_point * T8_ECLASS_MAX_DIM; + is_inside[i_point] = t8_vertex_point_inside (zeros, points + offset, tolerance); + } } /** Update a possible internal data buffer for per tree data.