Getting physical coordinates in shell2d geometry

[stichou]

Hello,

I am learning to use p4est and in my particular case I am interested to work with a 2D disk. I am getting inspired from simple2.c and step3.c examples.

If I well understood we cannot access directly the physical coordinates from a quadrant so we have to call first p4est_qcoord_to_vertex which maps quadrant->x and quadrant->y to the physical coordinate system xyz. The problem I have is that the coordinates I get don't make sense with respect to my problem. Indeed, I am creating a disk with inner and outer radius of 10 and 70, respectively, but when I print xyz[0] and xyz[1] I get values between -1 and 1.75. Am I doing something wrong ?

Code

static void
init_fn (p4est_t * p4est, p4est_topidx_t which_tree,
         p4est_quadrant_t * quadrant)
{
  data_t          *data = (data_t *) quadrant->p.user_data;

  p4est_qcoord_to_vertex (p4est->connectivity, which_tree,
                          quadrant->x, quadrant->y, xyz);

  data->rho = which_tree;

  printf("x=%lf, y=%lf\n", xyz[0], xyz[1]);

}
int
main (int argc, char **argv)
{
  ...
    /* Connectivity and geometry */
  conn = p4est_connectivity_new_shell2d ();
  geom = p4est_geometry_new_shell2d (conn, 70.0, 10.0);

  /* Creates the forest */
  p4est = p4est_new_ext (mpicomm, conn, 0, min_level, fill_uniform, sizeof (data_t), init_fn, geom);

  ...
}

Output

x=-1.000000, y=1.000000
x=-0.750000, y=1.000000
x=-1.000000, y=1.250000
x=-0.750000, y=1.250000
x=-0.500000, y=1.000000
x=-0.250000, y=1.000000
x=-0.500000, y=1.250000
x=-0.250000, y=1.250000
x=-1.000000, y=1.500000

...

[pkestene]

p4est_qcoord_to_vertex allows you to map the logical quadrant coordinates to the vertex space (as defined in the chosen tree connectivity)., this is perfectly valid as long as you don't use a geometry_t object which provides an additional mapping.

In that case, if geom != NULL, the way to get the physical coordinates is slightly different, you need to use function geom->X which maps the logical quadrant coordinates (rescaled to unit cube, using doubles not integers) to the physical space.

Have a look, e.g. at example simple2.c : https://github.com/cburstedde/p4est/blob/master/example/simple/simple2.c#L216
this a refine callback function, that uses the physical coordinates to build the refine condition, and it gives you an example on how to compute physical coordinates from quadrant coordinates. In this example xyz is the vector of logical quadrant coordinates (rescaled to unit cube), and XYZ the vector of physical coordinates you're looking for.

[stichou]

Thank you. It's clearer and now my code works.

I am also interested in getting the physical length of a quadrant in all directions. I should employ same technique isn't it ?

double L=P4EST_QUADRANT_LEN (quadrant->level) / P4EST_ROOT_LEN;
dxdydz[0] = L;
dxdydz[1] = L;

geom->X (geom, which_tree, dxdydz, dXdYdZ);

[pkestene]

No, you should create 2 points that are separated by L in the logical space, then transform these2 points, and finally compute distance in the physical space between the transformed point.

[stichou]

Ok thank you.