From 4b3a07f396e54bcdc2009b023701f3158d9026c6 Mon Sep 17 00:00:00 2001 From: Pierre Kestener Date: Fri, 18 Oct 2024 22:16:30 +0200 Subject: [PATCH] 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 | 39 ++++++++++++- src/p4est_geometry.c | 115 +++++++++++++++++++++++++++++++++++++++ src/p4est_geometry.h | 27 +++++++++ 3 files changed, 179 insertions(+), 2 deletions(-) diff --git a/example/simple/simple2.c b/example/simple/simple2.c index ed191ec83..a4a2b8a11 100644 --- a/example/simple/simple2.c +++ b/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: \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 (); } diff --git a/src/p4est_geometry.c b/src/p4est_geometry.c index 8e6e97185..64eb74f6b 100644 --- a/src/p4est_geometry.c +++ b/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 diff --git a/src/p4est_geometry.h b/src/p4est_geometry.h index 4f30447d1..2303e38e0 100644 --- a/src/p4est_geometry.h +++ b/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.