lc = 0.1; // was 0.25
x_center = 0;
y_center = 0;
z_center = 0;

// Define the inner sphere
R = 1;

Point(1) = {x_center, y_center, z_center, lc};
Point(2) = {x_center - R, y_center, z_center, lc};
Point(4) = {x_center, y_center - R, z_center, lc};
Point(5) = {x_center + R, y_center, z_center, lc};
Point(8) = {x_center, y_center, z_center - R, lc};
Point(11) = {x_center, y_center + R, z_center, lc};
Point(14) = {x_center, y_center, z_center + R, lc};
Circle (1) = {2, 1, 4} Plane{0, 0, 1};
Circle (2) = {4, 1, 5} Plane{0, 0, 1};
Circle (3) = {2, 1, 8} Plane{0, 0, 1};
Circle (4) = {4, 1, 8} Plane{0, 0, 1};
Circle (6) = {2, 1, 11} Plane{0, 0, 1};
Circle (7) = {8, 1, 11} Plane{0, 0, 1};
Circle (9) = {2, 1, 14} Plane{0, 0, 1};
Circle (10) = {11, 1, 14} Plane{0, 0, 1};
Circle (13) = {14, 1, 4} Plane{0, 0, 1};
Circle (15) = {8, 1, 5} Plane{0, 0, 1};
Circle (18) = {11, 1, 5} Plane{0, 0, 1};
Circle (21) = {14, 1, 5} Plane{0, 0, 1};
Line Loop (1000005) = {1, 4, -3};
Ruled Surface (5) = {1000005};
Line Loop (1000008) = {3, 7, -6};
Ruled Surface (8) = {1000008};
Line Loop (1000011) = {6, 10, -9};
Ruled Surface (11) = {1000011};
Line Loop (1000014) = {9, 13, -1};
Ruled Surface (14) = {1000014};
Line Loop (1000017) = {-15, -4, 2};
Ruled Surface (17) = {1000017};
Line Loop (1000020) = {-18, -7, 15};
Ruled Surface (20) = {1000020};
Line Loop (1000023) = {-21, -10, 18};
Ruled Surface (23) = {1000023};
Line Loop (1000026) = {-2, -13, 21};
Ruled Surface (26) = {1000026};

Surface Loop(1) = { 5, 8, 11, 14, 17, 20, 23, 26 };

// The cylinder
zcyl=-15;
xcyl=0;
ycyl=0;
radius=1.0;
DZ=0.1;


// The circular cross-section
For j In {0:4}
  // Calculating the points on the circle
  omega=j*2*Pi/4;
  xm[j]=(xcyl+radius*Cos(omega));
  ym[j]=(ycyl+radius*Sin(omega));
  Printf("%g %g",xm[j],ym[j]);
  Point(2000 + j)={xm[j],ym[j],zcyl-DZ/2,lc};
EndFor
Point(2005)={xcyl,ycyl,zcyl-DZ/2,lc};

// Creating a  "Circle"
Circle(2006)={2000,2005,2001};
Circle(2007)={2001,2005,2002};
Circle(2008)={2002,2005,2003};
Circle(2009)={2003,2005,2000};

// Defining a line loop, necessary for making a surface
CircleLoop_zm=newl; Line Loop(CircleLoop_zm) = {2006,2007,2008,2009};

// Defining the surface
CircleSurf_zm=news; Plane Surface(CircleSurf_zm)={CircleLoop_zm};

// When making an extruded line or surface it is automatically meshed identically to the parent surface.
out[]=Extrude{0,0,DZ}{ Line {2006,2007,2008,2009}; };

// out contains the extruded elements. index 0,4,8,12 are the extruded lines which we now use to create the
// other end of the cylinder.
CircleLoop_zp=newl; Line Loop(CircleLoop_zp) ={out[0],out[4],out[8],out[12]};
CircleSurf_zp=news; Plane Surface(CircleSurf_zp)={CircleLoop_zp};

// The two surfaces now need to be forced to get identical nodes. This is done with the Periodic Surface command
Periodic Surface CircleSurf_zm {2006,2007,2008,2009} = CircleSurf_zp {out[0],out[4],out[8],out[12]};

//out also contain the extruded sides of the surface at index 1,5,9,13. We use them to make the entire cylinder
Surface Loop(2) = {out[1],out[5],out[9],out[13],CircleSurf_zm,CircleSurf_zp} ;



// Define the outer sphere
lc = 3.0;
R = 50;

Point(1001) = {x_center, y_center, z_center, lc};
Point(1002) = {x_center - R, y_center, z_center, lc};
Point(1004) = {x_center, y_center - R, z_center, lc};
Point(1005) = {x_center + R, y_center, z_center, lc};
Point(1008) = {x_center, y_center, z_center - R, lc};
Point(1011) = {x_center, y_center + R, z_center, lc};
Point(1014) = {x_center, y_center, z_center + R, lc};
Circle (1001) = {1002, 1001, 1004} Plane{0, 0, 1};
Circle (1002) = {1004, 1001, 1005} Plane{0, 0, 1};
Circle (1003) = {1002, 1001, 1008} Plane{0, 0, 1};
Circle (1004) = {1004, 1001, 1008} Plane{0, 0, 1};
Circle (1006) = {1002, 1001, 1011} Plane{0, 0, 1};
Circle (1007) = {1008, 1001, 1011} Plane{0, 0, 1};
Circle (1009) = {1002, 1001, 1014} Plane{0, 0, 1};
Circle (1010) = {1011, 1001, 1014} Plane{0, 0, 1};
Circle (1013) = {1014, 1001, 1004} Plane{0, 0, 1};
Circle (1015) = {1008, 1001, 1005} Plane{0, 0, 1};
Circle (1018) = {1011, 1001, 1005} Plane{0, 0, 1};
Circle (1021) = {1014, 1001, 1005} Plane{0, 0, 1};
Line Loop (2000005) = {1001, 1004, -1003};
Ruled Surface (1005) = {2000005};
Line Loop (2000008) = {1003, 1007, -1006};
Ruled Surface (1008) = {2000008};
Line Loop (2000011) = {1006, 1010, -1009};
Ruled Surface (1011) = {2000011};
Line Loop (2000014) = {1009, 1013, -1001};
Ruled Surface (1014) = {2000014};
Line Loop (2000017) = {-1015, -1004, 1002};
Ruled Surface (1017) = {2000017};
Line Loop (2000020) = {-1018, -1007, 1015};
Ruled Surface (1020) = {2000020};
Line Loop (2000023) = {-1021, -1010, 1018};
Ruled Surface (1023) = {2000023};
Line Loop (2000026) = {-1002, -1013, 1021};
Ruled Surface (1026) = {2000026};

Surface Loop(3) = { 1005, 1008, 1011, 1014, 1017, 1020, 1023, 1026 };

// Define the volume between the two spheres
Volume(1) = { 3, 2, 1 };


