From 8862c0e27d1f91581e8cca803f1d9c08668ffca3 Mon Sep 17 00:00:00 2001 From: Pierre Kestener Date: Sun, 13 Oct 2024 19:33:48 +0200 Subject: [PATCH] Add connectivity and geometry for mapping a spherical shell using a two-tree pillow. --- example/simple/simple3.c | 11 +++++- src/p4est_connectivity.c | 3 ++ src/p4est_geometry.c | 9 +---- src/p8est_connectivity.c | 36 +++++++++++++++++ src/p8est_connectivity.h | 6 +++ src/p8est_geometry.c | 84 ++++++++++++++++++++++++++++++++++++++++ src/p8est_geometry.h | 14 +++++++ 7 files changed, 155 insertions(+), 8 deletions(-) diff --git a/example/simple/simple3.c b/example/simple/simple3.c index f832ee918..986d68612 100644 --- a/example/simple/simple3.c +++ b/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: \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); diff --git a/src/p4est_connectivity.c b/src/p4est_connectivity.c index 9ba5a5b6b..72322718e 100644 --- a/src/p4est_connectivity.c +++ b/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 (); } diff --git a/src/p4est_geometry.c b/src/p4est_geometry.c index 9a080f5d0..8e6e97185 100644 --- a/src/p4est_geometry.c +++ b/src/p4est_geometry.c @@ -36,13 +36,13 @@ #include #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] */ diff --git a/src/p8est_connectivity.c b/src/p8est_connectivity.c index 977d2b45c..4314dbad2 100644 --- a/src/p8est_connectivity.c +++ b/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) { diff --git a/src/p8est_connectivity.h b/src/p8est_connectivity.h index 14e2ec282..517b31112 100644 --- a/src/p8est_connectivity.h +++ b/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. */ diff --git a/src/p8est_geometry.c b/src/p8est_geometry.c index a600c4cdb..dc66464c2 100644 --- a/src/p8est_geometry.c +++ b/src/p8est_geometry.c @@ -35,6 +35,7 @@ typedef enum { P8EST_GEOMETRY_BUILTIN_MAGIC = 0x30F3F8DF, + P8EST_GEOMETRY_BUILTIN_PILLOW3D, P8EST_GEOMETRY_BUILTIN_SHELL, P8EST_GEOMETRY_BUILTIN_SPHERE, P8EST_GEOMETRY_BUILTIN_TORUS, @@ -42,6 +43,13 @@ typedef enum } 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, diff --git a/src/p8est_geometry.h b/src/p8est_geometry.h index 158d74f72..c4122047b 100644 --- a/src/p8est_geometry.h +++ b/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.