Basilisk source code (http://basilisk.fr/src/)

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});
}