Radial/cylindrical coordinates

This implements the radial coordinate mapping illustrated below.

Radial coordinate mapping

Radial coordinate mapping

It works in 1D, 2D and 3D. The 3D version corresponds to cylindrical coordinates since the z-coordinate is unchanged.

The only parameter is dθ, the total angle of the sector.

double dtheta = π/3.;

For convenience we add definitions for the radial and angular coordinates (r,θ).

map {
  double r = x, θ = y*dtheta/L0;

event metric (i = 0) {

We initialise the scale factors, taking care to first allocate the fields if they are still constant.

  if (is_constant(cm)) {
    scalar * l = list_copy (all);
    cm = new scalar;
    free (all);
    all = list_concat ({cm}, l);
    free (l);
  if (is_constant(fm.x)) {
    scalar * l = list_copy (all);
    fm = new face vector;
    free (all);
    all = list_concat ((scalar *){fm}, l);
    free (l);

The area (in 2D) of a mapped element is the area of an annulus defined by the two radii rΔ/2 and r+Δ/2, divided by the total number of sectors N=2πL0/(dθΔ), this gives π[(r+Δ/2)2(rΔ/2)2]2πL0/(dθΔ)=2πrΔ2πL0/(dθΔ)=rdθL0Δ2 By definition, the (area) metric factor cm is the mapped area divided by the unmapped area Δ2.

  scalar cmv = cm;
    cmv[] = r*dtheta/L0;

It is important to set proper boundary conditions, in particular when refining the grid.

  cm[left] = dirichlet (r*dtheta/L0);
  cm[right] = dirichlet (r*dtheta/L0);

The (length) metric factor fm is the ratio of the mapped length of a face to its unmapped length Δ. In the present case, it is unity for all dimensions except for the x coordinates for which it is the ratio of the arclength by the unmapped length Δ. We also set a small minimal value to avoid division by zero, in the case of a vanishing inner radius.

  face vector fmv = fm;
    fmv.x[] = 1.;
    fmv.x[] = max(r*dtheta/L0, 1e-20);
  boundary ({cm, fm});