p4est_bits: compute quadrant boundary coordinates

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,

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.

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

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));
+}

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.

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;