From b87d0e79e487994d5a5160b62e0fc2747e88ea8c Mon Sep 17 00:00:00 2001 From: Carsten Burstedde Date: Tue, 21 Jan 2025 21:12:48 +0100 Subject: [PATCH] p4est_bits: compute quadrant boundary coordinates --- src/p4est_bits.c | 42 ++++++++++++++++++++++++++++++++++++++++++ src/p4est_bits.h | 19 ++++++++++++++++++- src/p4est_to_p8est.h | 3 +++ src/p8est_bits.c | 28 ++++++++++++++++++++++++++++ src/p8est_bits.h | 28 +++++++++++++++++++++++++++- test/test_quadrants2.c | 35 +++++++++++++++++++++++++++++++++++ 6 files changed, 153 insertions(+), 2 deletions(-) diff --git a/src/p4est_bits.c b/src/p4est_bits.c index 61ac7fdde..376733bec 100644 --- a/src/p4est_bits.c +++ b/src/p4est_bits.c @@ -1453,6 +1453,31 @@ p4est_quadrant_all_face_neighbors (const p4est_quadrant_t * q, } } +void +p4est_quadrant_face_coordinates (const p4est_quadrant_t * q, int face, + p4est_qcoord_t coords[]) +{ + /* compute half length of qudrant: legal even when at P4EST_QMAXLEVEL */ + const p4est_qcoord_t qh = P4EST_QUADRANT_LEN (q->level); + const p4est_qcoord_t qhh = qh >> 1; + + /* store the coordinate axis normal to the face */ + const int faceh = face >> 1; + + P4EST_ASSERT (0 <= face && face < P4EST_FACES); + P4EST_ASSERT (0 <= faceh && faceh < P4EST_DIM); + P4EST_ASSERT (p4est_quadrant_is_valid (q)); + + /* for the coordinate axis normal to the face, choose shift of 0 or qh */ + /* for a coordinate axis parallel to the face, shift by half */ + coords[0] = q->x + ((faceh == 0) ? ((face & 1) ? qh : 0) : qhh); + coords[1] = q->y + ((faceh == 1) ? ((face & 1) ? qh : 0) : qhh); +#ifdef P4_TO_P8 + coords[2] = q->z + ((faceh == 2) ? ((face & 1) ? qh : 0) : qhh); +#endif + P4EST_ASSERT (p4est_coordinates_is_valid (coords, q->level + 1)); +} + void p4est_quadrant_corner_neighbor (const p4est_quadrant_t * q, int corner, p4est_quadrant_t * r) @@ -1670,6 +1695,23 @@ p4est_quadrant_corner_node (const p4est_quadrant_t * q, P4EST_ASSERT (p4est_quadrant_is_node (r, 0)); } +void +p4est_quadrant_corner_coordinates (const p4est_quadrant_t * q, int corner, + p4est_qcoord_t coords[]) +{ + const p4est_qcoord_t qh = P4EST_QUADRANT_LEN (q->level); + + P4EST_ASSERT (0 <= corner && corner < P4EST_CHILDREN); + P4EST_ASSERT (p4est_quadrant_is_valid (q)); + + coords[0] = q->x + ((corner & 1) ? qh : 0); + coords[1] = q->y + ((corner & 2) ? qh : 0); +#ifdef P4_TO_P8 + coords[2] = q->z + ((corner & 4) ? qh : 0); +#endif + P4EST_ASSERT (p4est_coordinates_is_valid (coords, q->level)); +} + void p4est_quadrant_children (const p4est_quadrant_t * q, p4est_quadrant_t * c0, p4est_quadrant_t * c1, diff --git a/src/p4est_bits.h b/src/p4est_bits.h index ed1228001..472ffdae2 100644 --- a/src/p4est_bits.h +++ b/src/p4est_bits.h @@ -493,6 +493,15 @@ void p4est_quadrant_all_face_neighbors (const p4est_quadrant_t * q, int face, p4est_quadrant_t n[]); +/** Compute the coordinates of a specific quadrant face's midpoint. + * \param [in] q Input quadrant, must be valid. + * \param [in] face The face of which the midpoint coordinates + * are computed. + * \param [out] coord 2D mid-face coordinates are in/on the unit tree. + */ +void p4est_quadrant_face_coordinates + (const p4est_quadrant_t * q, int face, p4est_qcoord_t coords[]); + /** Compute the corner neighbor of a quadrant. * \param [in] q Input quadrant, must be valid. * \param [in] corner The corner across which to generate the neighbor. @@ -541,7 +550,7 @@ void p4est_quadrant_half_corner_neighbor (const p4est_quadrant_t * r); -/** Compute the corner node of a quadrant. +/** Compute a corner node of a quadrant. * \param [in] q Input quadrant, must be valid. * \param [in] corner The corner across which to generate the neighbor. * \param [in,out] r Node that will not be clamped inside. @@ -551,6 +560,14 @@ void p4est_quadrant_corner_node (const p4est_quadrant_t * q, int corner, p4est_quadrant_t * r); +/** Compute the coordinates of a specific quadrant corner. + * \param [in] q Input quadrant, must be valid. + * \param [in] corner The corner for which the coordinates are computed. + * \param [out] coord 2D corner coordinates are in/on the unit tree. + */ +void p4est_quadrant_corner_coordinates + (const p4est_quadrant_t * q, int corner, p4est_qcoord_t coords[]); + /** Compute the 4 children of a quadrant. * \param [in] q Input quadrant. * \param [in,out] c0 First computed child. diff --git a/src/p4est_to_p8est.h b/src/p4est_to_p8est.h index 92391976f..1a0012555 100644 --- a/src/p4est_to_p8est.h +++ b/src/p4est_to_p8est.h @@ -376,12 +376,15 @@ #define p4est_quadrant_face_neighbor_extra p8est_quadrant_face_neighbor_extra #define p4est_quadrant_half_face_neighbors p8est_quadrant_half_face_neighbors #define p4est_quadrant_all_face_neighbors p8est_quadrant_all_face_neighbors +#define p4est_quadrant_face_coordinates p8est_quadrant_face_coordinates #define p4est_quadrant_corner_neighbor p8est_quadrant_corner_neighbor #define p4est_quadrant_corner_neighbor_extra \ p8est_quadrant_corner_neighbor_extra #define p4est_quadrant_half_corner_neighbor \ p8est_quadrant_half_corner_neighbor #define p4est_quadrant_corner_node p8est_quadrant_corner_node +#define p4est_quadrant_corner_coordinates \ + p8est_quadrant_corner_coordinates #define p4est_quadrant_children p8est_quadrant_children #define p4est_quadrant_childrenv p8est_quadrant_childrenv #define p4est_quadrant_childrenpv p8est_quadrant_childrenpv diff --git a/src/p8est_bits.c b/src/p8est_bits.c index eb2352955..78b0c9994 100644 --- a/src/p8est_bits.c +++ b/src/p8est_bits.c @@ -687,3 +687,31 @@ p8est_quadrant_shift_edge (const p4est_quadrant_t * q, P4EST_ASSERT (rdown == NULL || p8est_quadrant_touches_edge (rdown, edge, 1)); } + +void +p8est_quadrant_edge_coordinates (const p8est_quadrant_t * q, int edge, + p4est_qcoord_t coords[]) +{ + /* compute half length of qudrant: legal even when at P4EST_QMAXLEVEL */ + const p4est_qcoord_t qh = P4EST_QUADRANT_LEN (q->level); + const p4est_qcoord_t qhh = qh >> 1; + + /* store the coordinate axis parallel to the edge */ + const int edgeh = edge >> 2; + + /* the y coordinate is the first for an x-edge and the second for a z-edge */ + const int s = (edgeh == 0) ? 1 : 2; + + P4EST_ASSERT (0 <= edge && edge < P8EST_EDGES); + P4EST_ASSERT (0 <= edgeh && edgeh < P4EST_DIM); + P4EST_ASSERT (p4est_quadrant_is_valid (q)); + + /* for the coordinate axis parallel to the edge, shift by half */ + /* for a coordinate axis normal to the edge, choose shift of 0 or qh */ + coords[0] = q->x + ((edgeh == 0) ? qhh : ((edge & 1) ? qh : 0)); + coords[1] = q->y + ((edgeh == 1) ? qhh : ((edge & s) ? qh : 0)); +#ifdef P4_TO_P8 + coords[2] = q->z + ((edgeh == 2) ? qhh : ((edge & 2) ? qh : 0)); +#endif + P4EST_ASSERT (p4est_coordinates_is_valid (coords, q->level + 1)); +} diff --git a/src/p8est_bits.h b/src/p8est_bits.h index 6e7df21c5..f1986c3cc 100644 --- a/src/p8est_bits.h +++ b/src/p8est_bits.h @@ -514,6 +514,15 @@ void p8est_quadrant_all_face_neighbors (const p8est_quadrant_t * q, int face, p8est_quadrant_t n[]); +/** Compute the coordinates of a specific quadrant face's midpoint. + * \param [in] q Input quadrant, must be valid. + * \param [in] face The face of which the midpoint coordinates + * are computed. + * \param [out] coord 3D mid-face coordinates are in/on the unit tree. + */ +void p8est_quadrant_face_coordinates + (const p8est_quadrant_t * q, int face, p4est_qcoord_t coords[]); + /** Compute the edge neighbor of a quadrant. * \param [in] q Input quadrant, must be valid. * \param [in] edge The edge across which to generate the neighbor. @@ -549,6 +558,15 @@ void p8est_quadrant_edge_neighbor_extra (const p8est_quadrant_t p8est_connectivity_t * conn); +/** Compute the coordinates of a specific quadrant edge's midpoint. + * \param [in] q Input quadrant, must be valid. + * \param [in] edge The edge of which the midpoint coordinates + * are computed. + * \param [out] coord 3D mid-edge coordinates are in/on the unit tree. + */ +void p8est_quadrant_edge_coordinates + (const p8est_quadrant_t * q, int edge, p4est_qcoord_t coords[]); + /** Compute the corner neighbor of a quadrant. * \param [in] q Input quadrant, must be valid. * \param [in] corner The corner across which to generate the neighbor. @@ -597,7 +615,7 @@ void p8est_quadrant_half_corner_neighbor (const p8est_quadrant_t * r); -/** Compute the corner node of a quadrant. +/** Compute a corner node of a quadrant. * \param [in] q Input quadrant, must be valid. * \param [in] corner The corner across which to generate the neighbor. * \param [in,out] r Node that will not be clamped inside. @@ -607,6 +625,14 @@ void p8est_quadrant_corner_node (const p8est_quadrant_t * q, int corner, p8est_quadrant_t * r); +/** Compute the coordinates of a specific quadrant corner. + * \param [in] q Input quadrant, must be valid. + * \param [in] corner The corner for which the coordinates are computed. + * \param [out] coord 3D corner coordinates are in/on the unit tree. + */ +void p8est_quadrant_corner_coordinates + (const p8est_quadrant_t * q, int corner, p4est_qcoord_t coords[]); + /** Compute the 8 children of a quadrant. * \param [in] q Input quadrant. * \param [in,out] c0 First computed child. diff --git a/test/test_quadrants2.c b/test/test_quadrants2.c index 1b750fccd..754f3294f 100644 --- a/test/test_quadrants2.c +++ b/test/test_quadrants2.c @@ -148,6 +148,29 @@ check_predecessor_successor (const p4est_quadrant_t * q) SC_CHECK_ABORT (p4est_quadrant_is_equal (&temp2, q), "successor"); } +static void +check_coordinates (const p4est_quadrant_t * q) +{ + int face; +#ifdef P4_TO_P8 + int edge; +#endif + int corner; + p4est_qcoord_t coord[P4EST_DIM]; + + for (face = 0; face < P4EST_FACES; ++face) { + p4est_quadrant_face_coordinates (q, face, coord); + } +#ifdef P4_TO_P8 + for (edge = 0; edge < P8EST_EDGES; ++edge) { + p8est_quadrant_edge_coordinates (q, edge, coord); + } +#endif + for (corner = 0; corner < P4EST_CHILDREN; ++corner) { + p4est_quadrant_corner_coordinates (q, corner, coord); + } +} + #define NEG_ONE_MAXL (~((((p4est_qcoord_t) 1) << P4EST_MAXLEVEL) - 1)) #define NEG_ONE_MAXLM1 (~((((p4est_qcoord_t) 1) << (P4EST_MAXLEVEL - 1)) - 1)) #define NEG_ONE_MAXLP1 \ @@ -200,6 +223,9 @@ main (int argc, char **argv) for (iz = 0; iz < t1->quadrants.elem_count; ++iz) { q1 = p4est_quadrant_array_index (&t1->quadrants, iz); + /* test coordinates of all boundary objects */ + check_coordinates (q1); + /* test the index conversion */ index1 = p4est_quadrant_linear_id (q1, (int) q1->level); p4est_quadrant_set_morton (&r, (int) q1->level, index1); @@ -458,6 +484,15 @@ main (int argc, char **argv) P4EST_QUADRANT_INIT (&P); P4EST_QUADRANT_INIT (&Q); + /* test quadrant coordinates for first and last tree cornre */ + p4est_quadrant_set_morton (&A, 0, 0); + check_coordinates (&A); + p4est_quadrant_first_descendant (&A, &B, P4EST_QMAXLEVEL); + check_coordinates (&B); + p4est_quadrant_last_descendant (&A, &B, P4EST_QMAXLEVEL); + check_coordinates (&B); + + /* go through several hand-crafted quadrants */ A.x = NEG_ONE_MAXL; A.y = NEG_ONE_MAXL; A.level = 0;