# hub.darcs.net :: basilisk -> basilisk -> files

Basilisk source code

## root / src / radial.h.page

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83  /** # Radial/cylindrical coordinates This implements the radial coordinate mapping illustrated below. ![Radial coordinate mapping](figures/radial.svg) 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\theta$, the total angle of the sector. */ double dtheta = pi/3.; /** For convenience we add definitions for the radial and angular coordinates $(r, \theta)$. */ map { double r = x, theta = y*dtheta/L0; NOT_UNUSED(r); NOT_UNUSED(theta); } 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-\Delta/2$ and $r+\Delta/2$, divided by the total number of sectors $N=2\pi L0/(d\theta\Delta)$, this gives $$\frac{\pi [(r + \Delta / 2)^2 - (r - \Delta / 2)^2]}{2 \pi L 0 / (d \theta \Delta)} = \frac{2 \pi r \Delta}{2 \pi L 0 / (d \theta \Delta)} = \frac{rd \theta}{L 0} \Delta^2$$ By definition, the (area) metric factor cm is the mapped area divided by the unmapped area $\Delta^2$. */ scalar cmv = cm; foreach() 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 $\Delta$. 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 $\Delta$. We also set a small minimal value to avoid division by zero, in the case of a vanishing inner radius. */ face vector fmv = fm; foreach_face() fmv.x[] = 1.; foreach_face(x) fmv.x[] = max(r*dtheta/L0, 1e-20); boundary ({cm, fm}); }