 dy = 0.2;
 dz = 0.2;
 size_b = 0.3;
-// Each cylinder volume gets number 1, ..., n_cylinders. These markers are also
-// inherited by surfaces that bound the cylinder except the surface that is
-// shared by 2 cylinders. The label for the cylinder is (volume label +
-// n_cylinders). The outer volume, i.e. bounding box - cylinder union is tagged
-// as 0
+// The cylinders are numbered (the volume) starting from on first growing in 
+// y-direction then in z-direction. The volume id determines also id of all the
+// surfaces that are not shared. For the shared surfaces the y-shared surface is
+// volume_id + (z_cylinders*y_cylinders). In z direction this is 
+// volume_id + 2*(z_cylinders*y_cylinders). Note that in this policy it is the first
+// cylinder which hits the surface that gets to name it.
+// The outer volume, i.e. bounding box - cylinder union is tagged as 0
 //! Cut here
 
 // The first cylinder has center at x0, y0, z0

Surface Loop(sl) = cylinder[];
v = newv;
Volume(v) = {sl};

Point(12) = {R, 0, 0, size_R};
Point(13) = {q, Sqrt(R^2-q^2), 0, size_R};

Point(14) = {0, R+l/2, 0, size_l};
Point(15) = {0, R+l/2, t, size_l};
Point(16) = {t, R+l/2, 0, size_l};
Point(17) = {0, R+(q-t), t, size_l};

Point(18) = {0, 0, 0, size_R};
