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

Basilisk source code

## root / src / axi.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 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106  /** # Axisymmetric coordinates For problems with a symmetry of revolution around the $z$-axis of a [cylindrical coordinate system](http://en.wikipedia.org/wiki/Cylindrical_coordinate_system). The longitudinal coordinate ($z$-axis) is *x* and the radial coordinate ($\rho$- or $r$-axis) is *y*. Note that *y* (and so *Y0*) cannot be negative. We first define a macro which will be used in some geometry-specific code (e.g. [curvature computation](curvature.h)). */ #define AXI 1 /** On trees we need refinement functions. */ #if TREE static void refine_cm_axi (Point point, scalar cm) { fine(cm,0,0) = fine(cm,1,0) = y - Delta/4.; fine(cm,0,1) = fine(cm,1,1) = y + Delta/4.; } static void refine_face_x_axi (Point point, scalar fm) { if (!is_refined(neighbor(-1))) { fine(fm,0,0) = y - Delta/4.; fine(fm,0,1) = y + Delta/4.; } if (!is_refined(neighbor(1)) && neighbor(1).neighbors) { fine(fm,2,0) = y - Delta/4.; fine(fm,2,1) = y + Delta/4.; } fine(fm,1,0) = y - Delta/4.; fine(fm,1,1) = y + Delta/4.; } static void refine_face_y_axi (Point point, scalar fm) { if (!is_refined(neighbor(0,-1))) fine(fm,0,0) = fine(fm,1,0) = max(y - Delta/2., 1e-20); if (!is_refined(neighbor(0,1)) && neighbor(0,1).neighbors) fine(fm,0,2) = fine(fm,1,2) = y + Delta/2.; fine(fm,0,1) = fine(fm,1,1) = y; } #endif event metric (i = 0) { /** By default *cm* is a constant scalar field. To make it variable, we need to allocate a new field. We also move it at the begining of the list of variables: this is important to ensure the metric is defined before other fields. */ if (is_constant(cm)) { scalar * l = list_copy (all); cm = new scalar; free (all); all = list_concat ({cm}, l); free (l); } /** The volume/area of a cell is proportional to $r$ (i.e. $y$). We need to set boundary conditions at the top and bottom so that *cm* is interpolated properly when refining/coarsening the mesh. */ scalar cmv = cm; foreach() cmv[] = y; cm[top] = dirichlet(y); cm[bottom] = dirichlet(y); /** We do the same for the length scale factors. The "length" of faces on the axis of revolution is zero ($y=r=0$ on the axis). To avoid division by zero we set it to epsilon (note that mathematically the limit is well posed). */ if (is_constant(fm.x)) { scalar * l = list_copy (all); fm = new face vector; free (all); all = list_concat ((scalar *){fm}, l); free (l); } face vector fmv = fm; foreach_face() fmv.x[] = max(y, 1e-20); fm.t[top] = dirichlet(y); fm.t[bottom] = dirichlet(y); /** We set our refinement/prolongation functions on trees. */ #if TREE cm.refine = cm.prolongation = refine_cm_axi; fm.x.prolongation = refine_face_x_axi; fm.y.prolongation = refine_face_y_axi; #endif boundary ({cm, fm}); }