Add connectivity and geometry for mapping a spherical shell using a t…

…wo-tree pillow.

example/simple/simple3.c

@@ -32,6 +32,7 @@
  *        o twocubes  Two connected cubes.
  *        o twowrap   Two cubes with periodically identified far ends.
  *        o rotcubes  A collection of six connected rotated cubes.
+ *        o pillow3d  A 2-tree  discretization of a hollow sphere.
  *        o shell     A 24-tree discretization of a hollow sphere.
  *        o sphere    A 13-tree discretization of a solid sphere.
  */
@@ -56,6 +57,7 @@ typedef enum
   P8EST_CONFIG_TWOCUBES,
   P8EST_CONFIG_TWOWRAP,
   P8EST_CONFIG_ROTCUBES,
+  P8EST_CONFIG_PILLOW3D,
   P8EST_CONFIG_SHELL,
   P8EST_CONFIG_SPHERE,
   P8EST_CONFIG_TORUS,
@@ -198,7 +200,7 @@ main (int argc, char **argv)
   usage =
     "Arguments: <configuration> <level>\n"
     "   Configuration can be any of\n"
-    "      unit|brick|periodic|rotwrap|drop|twocubes|twowrap|rotcubes|shell|sphere|torus\n"
+    "      unit|brick|periodic|rotwrap|drop|twocubes|twowrap|rotcubes|pillow3d|shell|sphere|torus\n"
     "   Level controls the maximum depth of refinement\n";
   wrongusage = 0;
   config = P8EST_CONFIG_NULL;
@@ -230,6 +232,9 @@ main (int argc, char **argv)
     else if (!strcmp (argv[1], "rotcubes")) {
       config = P8EST_CONFIG_ROTCUBES;
     }
+    else if (!strcmp (argv[1], "pillow3d")) {
+      config = P8EST_CONFIG_PILLOW3D;
+    }
     else if (!strcmp (argv[1], "shell")) {
       config = P8EST_CONFIG_SHELL;
     }
@@ -282,6 +287,10 @@ main (int argc, char **argv)
   else if (config == P8EST_CONFIG_ROTCUBES) {
     connectivity = p8est_connectivity_new_rotcubes ();
   }
+  else if (config == P8EST_CONFIG_PILLOW3D) {
+    connectivity = p8est_connectivity_new_pillow3d ();
+    geom = p8est_geometry_new_pillow3d (connectivity, 1., .55);
+  }
   else if (config == P8EST_CONFIG_SHELL) {
     connectivity = p8est_connectivity_new_shell ();
     geom = p8est_geometry_new_shell (connectivity, 1., .55);

src/p4est_connectivity.c

@@ -2920,6 +2920,9 @@ p4est_connectivity_new_byname (const char *name)
   else if (!strcmp (name, "rotwrap")) {
     return p8est_connectivity_new_rotwrap ();
   }
+  else if (!strcmp (name, "pillow3d")) {
+    return p8est_connectivity_new_pillow3d ();
+  }
   else if (!strcmp (name, "shell")) {
     return p8est_connectivity_new_shell ();
   }

src/p4est_geometry.c

@@ -36,13 +36,13 @@
 #include <p8est_geometry.h>
 #endif
 
-#ifndef P4_TO_P8
-
 static int sign_d (double value)
 {
   return ((double)0 < value) - (value < (double)0);
 }
 
+#ifndef P4_TO_P8
+
 typedef enum
 {
   P4EST_GEOMETRY_BUILTIN_MAGIC = 0x20F2F8DE,
@@ -659,11 +659,6 @@ p4est_geometry_pillow_X (p4est_geometry_t * geom,
 
   Rsphere = pillow->R;
 
-  /* transform from the tree-local reference coordinates into the logical vertex space
-   * using bi/trilinear transformation.
-   */
-  p4est_geometry_connectivity_X (geom, which_tree, rst, xyz);
-
   sgnz = 2.0 * which_tree - 1;
 
   /* remap to [-1, 1] x [-1, 1] */

src/p8est_connectivity.c

@@ -605,6 +605,42 @@ p8est_connectivity_new_rotcubes (void)
                                       corner_to_tree, corner_to_corner);
 }
 
+p4est_connectivity_t *p8est_connectivity_new_pillow3d(void)
+{
+  const p4est_topidx_t num_vertices = 8;
+  const p4est_topidx_t num_trees = 2;
+  const p4est_topidx_t num_ett = 0;
+  const p4est_topidx_t num_ctt = 0;
+  const double        vertices[8 * 3] = {
+    0, 0, 0,
+    0, 1, 0,
+    1, 0, 0,
+    1, 1, 0,
+    0, 0, 1,
+    0, 1, 1,
+    1, 0, 1,
+    1, 1, 1,
+  };
+  const p4est_topidx_t tree_to_vertex[2 * 8] = {
+    0, 1, 2, 3, 4, 5, 6, 7,
+    0, 1, 2, 3, 4, 5, 6, 7,
+  };
+  const p4est_topidx_t tree_to_tree[2 * 6] = {
+    1, 1, 1, 1, 0, 0,
+    0, 0, 0, 0, 1, 1,
+  };
+  const int8_t        tree_to_face[2 * 6] = {
+    0, 1, 2, 3, 4, 5,
+    0, 1, 2, 3, 4, 5,
+  };
+
+  return p4est_connectivity_new_copy (num_vertices, num_trees, 0, 0,
+                                      vertices, tree_to_vertex,
+                                      tree_to_tree, tree_to_face,
+                                      NULL, &num_ett, NULL, NULL,
+                                      NULL, &num_ctt, NULL, NULL);
+}
+
 p4est_connectivity_t *
 p8est_connectivity_new_shell (void)
 {

src/p8est_connectivity.h

@@ -714,6 +714,12 @@ p8est_connectivity_t *p8est_connectivity_new_twowrap (void);
  */
 p8est_connectivity_t *p8est_connectivity_new_rotcubes (void);
 
+/** Create a connectivity structure for two trees on top of each other.
+ * This connectivity is meant to be used with \ref p8est_geometry_new_pillow3d
+ * to map a spherical shell.
+ */
+p8est_connectivity_t *p8est_connectivity_new_pillow3d (void);
+
 /** An m by n by p array with periodicity in x, y, and z if
  * periodic_a, periodic_b, and periodic_c are true, respectively.
  */

src/p8est_geometry.c

@@ -35,13 +35,21 @@
 typedef enum
 {
   P8EST_GEOMETRY_BUILTIN_MAGIC = 0x30F3F8DF,
+  P8EST_GEOMETRY_BUILTIN_PILLOW3D,
   P8EST_GEOMETRY_BUILTIN_SHELL,
   P8EST_GEOMETRY_BUILTIN_SPHERE,
   P8EST_GEOMETRY_BUILTIN_TORUS,
   P8EST_GEOMETRY_LAST
 }
 p8est_geometry_builtin_type_t;
 
+typedef struct p8est_geometry_builtin_pillow3d
+{
+  p8est_geometry_builtin_type_t type;
+  double              R2, R1;   /* outer/inner sphere radius */
+}
+p8est_geometry_builtin_pillow3d_t;
+
 typedef struct p8est_geometry_builtin_shell
 {
   p8est_geometry_builtin_type_t type;
@@ -77,6 +85,7 @@ typedef struct p8est_geometry_builtin
   union
   {
     p8est_geometry_builtin_type_t type;
+    p8est_geometry_builtin_pillow3d_t pillow3d;
     p8est_geometry_builtin_shell_t shell;
     p8est_geometry_builtin_sphere_t sphere;
     p8est_geometry_builtin_torus_t torus;
@@ -85,6 +94,81 @@ typedef struct p8est_geometry_builtin
 }
 p8est_geometry_builtin_t;
 
+static void
+p8est_geometry_pillow3d_X (p8est_geometry_t * geom,
+                           p4est_topidx_t which_tree,
+                           const double rst[3], double xyz[3])
+{
+  const struct p8est_geometry_builtin_pillow3d *pillow3d
+    = &((p8est_geometry_builtin_t *) geom)->p.pillow3d;
+  double              xc, yc, zc, R1, R2, abs_xc, abs_yc;;
+  double              xp, yp, zp;
+  double              d, D, R, Rz, center;
+
+  /* remap into [-1, 1] x [-1, 1] x [0, 1] */
+  xc = 2 * rst[0] - 1;
+  yc = 2 * rst[1] - 1;
+  zc = rst[2];
+
+  abs_xc = fabs(xc);
+  abs_yc = fabs(yc);
+
+  d = fmax(abs_xc, abs_yc);
+  d = fmax(d, 1e-10);
+
+  /* D = d / sqrt(2.0); */
+  D = d * (2 - d) / sqrt(2.0);
+  R = 1;
+
+  /* D = sin(M_PI * d / 2) / sqrt(2.0); */
+  /* R = sin(M_PI * d / 2); */
+
+  center = D - sqrt(R * R - D * D);
+
+  xp = D / d * abs_xc;
+  yp = D / d * abs_yc;
+
+  yp = abs_yc >= abs_xc ? center + sqrt(R * R - xp * xp) : yp;
+  xp = abs_xc >= abs_yc ? center + sqrt(R * R - yp * yp) : xp;
+
+  xp *= sign_d(xc);
+  yp *= sign_d(yc);
+  zp = sqrt(1 - (xp * xp + yp * yp));
+
+  /* tree 0 => upper hemisphere */
+  /* tree 1 => lower hemisphere */
+  zp = which_tree == 1 ? -zp : zp;
+
+  R1 = pillow3d->R1;
+  R2 = pillow3d->R2;
+
+  Rz = R1 + zc * (R2 - R1);
+
+  xyz[0] = Rz * xp;
+  xyz[1] = Rz * yp;
+  xyz[2] = Rz * zp;
+}
+
+p8est_geometry_t   *
+p8est_geometry_new_pillow3d (p8est_connectivity_t * conn, double R2, double R1)
+{
+  p8est_geometry_builtin_t *builtin;
+  struct p8est_geometry_builtin_pillow3d *pillow3d;
+
+  builtin = P4EST_ALLOC_ZERO (p8est_geometry_builtin_t, 1);
+
+  pillow3d = &builtin->p.pillow3d;
+  pillow3d->type = P8EST_GEOMETRY_BUILTIN_PILLOW3D;
+  pillow3d->R2 = R2;
+  pillow3d->R1 = R1;
+
+  builtin->geom.name = "p8est_pillow3d";
+  builtin->geom.user = conn;
+  builtin->geom.X = p8est_geometry_pillow3d_X;
+
+  return (p8est_geometry_t *) builtin;
+}
+
 static void
 p8est_geometry_shell_X (p8est_geometry_t * geom,
                         p4est_topidx_t which_tree,

src/p8est_geometry.h

@@ -108,6 +108,20 @@ void                p8est_geometry_connectivity_X (p8est_geometry_t *geom,
                                                    const double abc[3],
                                                    double xyz[3]);
 
+/** Create a geometry structure for the spherical shell of 2 trees.
+ * \param [in] conn Result of p8est_connectivity_new_pillow3d.
+ *                  We do NOT take ownership and expect it to stay alive.
+ * \param [in] R2   The outer radius of the shell.
+ * \param [in] R1   The inner radius of the shell.
+ * \return          Geometry structure; use with \ref p4est_geometry_destroy.
+ *
+ * \note this coordinate transformation is describe in "Logically rectangular
+ * grids and finite volume methods for PDEs in circular and spherical domains",
+ * Calhoun et al., https://doi.org/10.1137/060664094
+ */
+p8est_geometry_t   *p8est_geometry_new_pillow3d (p8est_connectivity_t * conn,
+                                                 double R2, double R1);
+
 /** Create a geometry structure for the spherical shell of 24 trees.
  * \param [in] conn Result of p8est_connectivity_new_shell or equivalent.
  *                  We do NOT take ownership and expect it to stay alive.