SetFactory("Built-in");

// Create the initial volume.
azi_frac = 6;
Point(1) = {0, 1, 0, 1.0};
Point(2) = {0, 0, 0, 1.0};
Line(1) = {2, 1};
Extrude {{1, 0, 0}, {0, 0, 0}, 2*Pi/azi_frac} {
  Curve{1}; 
}
Extrude {2, 0, 0} {
  Curve{3}; Curve{1}; Curve{2}; 
}
Curve Loop(2) = {13, -5, -9};
Plane Surface(17) = {2};
Surface Loop(1) = {8, 4, 12, 17, 16};
Volume(1) = {1};

// Print out all the surfaces.
all_surfs[] = Surface{:};
For i In {0 : #all_surfs[]-1}
  Printf("all_surfs[%g] = %g", i, all_surfs[i]);
EndFor

// Grab only the surfaces that are used in the definition of the volume.
// Should be the same as all_surfs.
src_surfs = Boundary{ Volume{1}; };
For i In {0 : #src_surfs[]-1}
  Printf("src_surfs[%g] = %g", i, src_surfs[i]);
EndFor

// We'll rotate all the surfaces of the initial volume, except for the
// one that will eventually be shared by the original volume and the new
// volume.
rot_surfs[] = {src_surfs[0], src_surfs[1], src_surfs[3], src_surfs[4]};

azi_frac_to_mesh = 3;
For i In {1 : azi_frac_to_mesh-1}

  // Rotate src_surfs, saving all the new surfaces in new_surfs.
  rot_angle = 2*Pi/azi_frac;
  new_surfs[] = Rotate {{1, 0, 0}, {0, 0, 0}, rot_angle} {Duplicata{Surface{rot_surfs[]};}};
  For s In {0 : #new_surfs[]-1}
    Printf("azi_frac %g, new_surfs[%g] = %g", i, s, new_surfs[s]);
  EndFor

  // Create the new volume.
  sl = newsl;
  Surface Loop(sl) = {new_surfs[], rot_surfs[3]};
  iso_v[i] = newv;
  Volume(iso_v[i]) = {sl};

  // Reset the rot_surfs for the next loop.
  rot_surfs[] = new_surfs[];

EndFor
