Add a 2d geometry to map the disk from unisquare connectivity.

The geometry is named pillow_disk. There are actually four variants of the
geometry which correspond to the four grid mappings shown in figure 3 from
the following publication:

"Logically rectangular grids and finite volume methods for PDEs in circular and
spherical domains", Calhoun et al, SIAM Review, volume 50, Issue 4, January
2008.

See https://doi.org/10.1137/060664094

example/simple/simple2.c

@@ -30,7 +30,8 @@
  *        o three     Refinement on a forest with three trees.
  *        o evil      Check second round of refinement with np=5 level=7
  *        o evil3     Check second round of refinement on three trees
- *        o pillow    Refinement on a 2-tree pillow-shaped domain.
+ *        o pillow      Refinement on a sphere (2-tree with pillow geometry).
+ *        o pillow_disk Refinement on a disk (1-tree with pillow geometry).
  *        o moebius   Refinement on a 5-tree Moebius band.
  *        o star      Refinement on a 6-tree star shaped domain.
  *        o cubed     Refinement on a 6-tree cubed sphere surface.
@@ -62,6 +63,7 @@ typedef enum
   P4EST_CONFIG_EVIL,
   P4EST_CONFIG_EVIL3,
   P4EST_CONFIG_PILLOW,
+  P4EST_CONFIG_PILLOW_DISK,
   P4EST_CONFIG_MOEBIUS,
   P4EST_CONFIG_STAR,
   P4EST_CONFIG_CUBED,
@@ -165,6 +167,15 @@ refine_normal_fn (p4est_t * p4est, p4est_topidx_t which_tree,
   return 1;
 }
 
+static int refine_uniform_fn(p4est_t *p4est, p4est_topidx_t which_tree,
+                             p4est_quadrant_t *quadrant)
+{
+  if ((int) quadrant->level >= refine_level) {
+    return 0;
+  }
+  return 1;
+}
+
 static int
 refine_evil_fn (p4est_t * p4est, p4est_topidx_t which_tree,
                 p4est_quadrant_t * quadrant)
@@ -293,7 +304,7 @@ main (int argc, char **argv)
   usage =
     "Arguments: <configuration> <level>\n"
     "   Configuration can be any of\n"
-    "      unit|brick|three|evil|evil3|pillow|moebius|\n"
+    "      unit|brick|three|evil|evil3|pillow|pillow_disk|moebius|\n"
     "         star|cubed|disk|xdisk|ydisk|pdisk|periodic|\n"
     "         rotwrap|circle|drop|icosahedron|shell2d|disk2d|bowtie|sphere2d\n"
     "   Level controls the maximum depth of refinement\n";
@@ -321,6 +332,9 @@ main (int argc, char **argv)
     else if (!strcmp (argv[1], "pillow")) {
       config = P4EST_CONFIG_PILLOW;
     }
+    else if (!strcmp (argv[1], "pillow_disk")) {
+      config = P4EST_CONFIG_PILLOW_DISK;
+    }
     else if (!strcmp (argv[1], "moebius")) {
       config = P4EST_CONFIG_MOEBIUS;
     }
@@ -392,6 +406,10 @@ main (int argc, char **argv)
     refine_fn = refine_icosahedron_fn;
     coarsen_fn = NULL;
   }
+  else if (config == P4EST_CONFIG_PILLOW_DISK) {
+    refine_fn = refine_uniform_fn;
+    coarsen_fn = NULL;
+  }
   else {
     refine_fn = refine_normal_fn;
     coarsen_fn = NULL;
@@ -416,6 +434,23 @@ main (int argc, char **argv)
     connectivity = p4est_connectivity_new_pillow ();
     geom = p4est_geometry_new_pillow (connectivity, R);
   }
+  else if (config == P4EST_CONFIG_PILLOW_DISK) {
+    double               R = 1.0;        /* disk radius default value */
+    int                  iconfig;
+    pillow_disk_config_t pconfig = FIG3A;
+
+    if (argc >= 4)
+      R = atof (argv[3]);
+    if (argc >= 5) {
+      iconfig = atoi (argv[4]);
+      if (iconfig >=0 && iconfig <=3) {
+        pconfig = (pillow_disk_config_t) iconfig;
+      }
+    }
+
+    connectivity = p4est_connectivity_new_unitsquare ();
+    geom = p4est_geometry_new_pillow_disk (connectivity, R, pconfig);
+  }
   else if (config == P4EST_CONFIG_MOEBIUS) {
     connectivity = p4est_connectivity_new_moebius ();
   }

src/p4est_geometry.c

@@ -51,6 +51,7 @@ typedef enum
   P4EST_GEOMETRY_BUILTIN_DISK2D,
   P4EST_GEOMETRY_BUILTIN_SPHERE2D,
   P4EST_GEOMETRY_BUILTIN_PILLOW,
+  P4EST_GEOMETRY_BUILTIN_PILLOW_DISK,
   P4EST_GEOMETRY_LAST
 }
 p4est_geometry_builtin_type_t;
@@ -92,6 +93,13 @@ typedef struct p4est_geometry_builtin_pillow
   double              R;        /* sphere radius */
 } p4est_geometry_builtin_pillow_t;
 
+typedef struct p4est_geometry_builtin_pillow_disk
+{
+  p4est_geometry_builtin_type_t type;
+  double              R;        /* disk radius */
+  pillow_disk_config_t config;
+} p4est_geometry_builtin_pillow_disk_t;
+
 typedef struct p4est_geometry_builtin
 {
   /** The geom member needs to come first; we cast to p4est_geometry_t * */
@@ -104,6 +112,7 @@ typedef struct p4est_geometry_builtin
     p4est_geometry_builtin_disk2d_t disk2d;
     p4est_geometry_builtin_sphere2d_t sphere2d;
     p4est_geometry_builtin_pillow_t pillow;
+    p4est_geometry_builtin_pillow_disk_t pillow_disk;
   }
   p;
 }
@@ -704,6 +713,112 @@ p4est_geometry_new_pillow (p4est_connectivity_t * conn, double R)
   return (p4est_geometry_t *) builtin;
 }                               /* p4est_geometry_new_pillow */
 
+/**
+ * geometric coordinate transformation for pillow_disk geometry.
+ *
+ * Define the geometric transformation from tree-local reference coordinates to the
+ * physical space.
+ *
+ * \param[in]  geom       associated geometry
+ * \param[in]  which_tree tree id inside forest
+ * \param[in]  rst        tree-local reference coordinates : [0,1]^2.
+ *                        Note: rst[2] is never accessed
+ * \param[out] xyz        Cartesian coordinates in physical space after geometry
+ *
+ */
+static void
+p4est_geometry_pillow_disk_X (p4est_geometry_t * geom,
+                              p4est_topidx_t which_tree,
+                              const double rst[3], double xyz[3])
+{
+  const struct p4est_geometry_builtin_pillow_disk *pillow_disk
+    = &((p4est_geometry_builtin_t *) geom)->p.pillow_disk;
+  double              Rdisk;
+
+  double absx, absy, d, D, R, xp, yp, center;
+  double r, w;
+
+  Rdisk = pillow_disk->R;
+
+  /* remap to [-1, 1] x [-1, 1] */
+  xyz[0] = 2 * rst[0] - 1;
+  xyz[1] = 2 * rst[1] - 1;
+  xyz[2] = 0;
+
+  absx = fabs(xyz[0]);
+  absy = fabs(xyz[1]);
+  d = fmax(absx, absy);
+
+  if (pillow_disk->config == FIG3D) {
+    r = sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
+    r = fmax(r, 1e-10);
+
+    xp = Rdisk * d * xyz[0] / r;
+    yp = Rdisk * d * xyz[1] / r;
+
+    w = d * d;
+
+    xyz[0] = w * xp + (1 - w) * xyz[0] / sqrt(2.0);
+    xyz[1] = w * yp + (1 - w) * xyz[1] / sqrt(2.0);
+    xyz[2] = 0.0;
+
+    return;
+
+  } else {
+
+  switch (pillow_disk->config) {
+  case FIG3A:
+    D = Rdisk * d / sqrt(2);
+    R = Rdisk * d;
+    break;
+  case FIG3B:
+    D = Rdisk * d / sqrt(2);
+    R = Rdisk;
+    break;
+  case FIG3C:
+    D = Rdisk * d * (2 - d) / sqrt(2);
+    R = Rdisk;
+    break;
+  default:
+    D = Rdisk * d / sqrt(2);
+    R = Rdisk * d;
+  }
+
+  center = D - sqrt(R*R - D*D);
+  xp = d > 0 ? D / d * absx : 0;
+  yp = d > 0 ? D / d * absy : 0;
+
+  yp = absy >= absx ? center + sqrt(R*R - xp*xp) : yp;
+  xp = absx >= absy ? center + sqrt(R*R - yp*yp) : xp;
+
+  xyz[0] = sign_d(xyz[0]) * xp;
+  xyz[1] = sign_d(xyz[1]) * yp;
+  xyz[2] = 0;
+
+  }
+}                               /* p4est_geometry_pillow_disk_X */
+
+p4est_geometry_t   *
+p4est_geometry_new_pillow_disk (p4est_connectivity_t * conn,
+                                double R, pillow_disk_config_t config)
+{
+  p4est_geometry_builtin_t *builtin;
+  struct p4est_geometry_builtin_pillow_disk *pillow_disk;
+
+  builtin = P4EST_ALLOC_ZERO (p4est_geometry_builtin_t, 1);
+
+  pillow_disk = &builtin->p.pillow_disk;
+  pillow_disk->type = P4EST_GEOMETRY_BUILTIN_PILLOW_DISK;
+  pillow_disk->R = R;
+  pillow_disk->config = config;
+
+  builtin->geom.name = "p4est_pillow_disk";
+  builtin->geom.user = conn;
+  builtin->geom.X = p4est_geometry_pillow_disk_X;
+
+  return (p4est_geometry_t *) builtin;
+}                               /* p4est_geometry_new_pillow_disk */
+
 #endif /* !P4_TO_P8 */
 
 static void

src/p4est_geometry.h

@@ -181,6 +181,33 @@ p4est_geometry_t   *p4est_geometry_new_sphere2d (p4est_connectivity_t * conn,
  */
 p4est_geometry_t   *p4est_geometry_new_pillow (p4est_connectivity_t *
                                                conn, double R);
+/** Characterize different mapping variants.
+ *
+ * The different mapping correspond to the ones used to produce figure 3 in the
+ * following publication:
+ *
+ * "Logically rectangular grids and finite volume methods for PDEs in circular
+ * and spherical domains", Calhoun et al, SIAM Review, volume 50, Issue 4, January 2008.
+ * https://doi.org/10.1137/060664094
+ */
+typedef enum
+{
+    FIG3A = 0,
+    FIG3B = 1,
+    FIG3C = 2,
+    FIG3D = 3,
+}
+pillow_disk_config_t;
+
+/** Create a geometry for mapping the disk using 2d connectivity unit.
+ *
+ * \param[in] conn      The result of \ref p4est_connectivity_new_unit.
+ * \param[in] R         The radius of the disk.
+ * \param[in] conf      The config to identify a mapping variant
+ */
+p4est_geometry_t   *p4est_geometry_new_pillow_disk (p4est_connectivity_t *
+                                                    conn, double R,
+                                                    pillow_disk_config_t conf);
 
 
 /** Compute node coordinates for a \ref p4est_lnodes structure.