// Trying to combine two volumes

// ---------               VOXEL             -------------------------------------

// Voxel: Parameters
L = 600; /* voxel length (x and y). z: L/2 */  
dx1 = 60;

// Create the four points at the base of the voxel
// Point(id. number) = {x,y,z,mesh size approx = dx}
// length/(dx)  = aprox #points in line
// elements in the volume =aprox (lentgh/dx)^3
Point(1) = { 0.5*L,  0.5*L, 0, dx1};
Point(2) = {-0.5*L,  0.5*L, 0, dx1};
Point(3) = {-0.5*L, -0.5*L, 0, dx1};
Point(4) = { 0.5*L, -0.5*L, 0, dx1};

// Create the lines between the points
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};

// Create a line loop with id. = newl. newl is the next line id number available
lvox = newl;
Line Loop(lvox) = {1, 2, 3, 4};

// Create the surface between these lines
svox=news; // news: next available surface id number
Plane Surface(svox) = {lvox};

// Extrude the surface
matvox = Extrude {0, 0, 0.5*L} { Surface{svox}; }; 


// Cylinder 1: Parameters     ------------------------------------------------------------
R = 50;  // For 8 elements on the circle: dx2=1.5R, with L = 60R. To see the cylinder empty, dx2 << R and R bigger
dx2 = 10;
dy=0; // Displacement in y axis
dz=2; // Displacement in z axis 

// Create the four points at the base circumference of the cylinder and the center. 
p0 = newp; Point(p0) = { 0,  0+dy, 0+dz, dx2};
p1 = newp; Point(p1) = { R,  0+dy, 0+dz, dx2};
p2 = newp; Point(p2) = { 0,  R+dy, 0+dz, dx2};
p3 = newp; Point(p3) = {-R,  0+dy, 0+dz, dx2};
p4 = newp; Point(p4) = { 0, -R+dy, 0+dz, dx2};

// Create the arcs between every two points. Arc is a region so search
// the next id number available using newreg. {start,center,end} 
c1 = newreg; Circle(c1) = {p1, p0, p2};
c2 = newreg; Circle(c2) = {p2, p0, p3};
c3 = newreg; Circle(c3) = {p3, p0, p4};
c4 = newreg; Circle(c4) = {p4, p0, p1};
// Create a new line loop by adding the four arcs with id lcirc=newl
lcirc = newl;
Line Loop(lcirc) = {c1, c2, c3, c4};

// Create the surface between these lines
scyl=news; // news: next available surface id number
Plane Surface(scyl) = {lcirc};

// Extrude the surface
matcyl = Extrude {0, 0, 0.25*L-(2*dz)} { Surface{scyl}; }; 


// Cylinder 2: Parameters      -----------------------------------------------
R = 40;
dx3 = 10;

// Create the four points at the base circumference of the cylinder and the center. 
p5 = newp; Point(p5) = { 0,  0,  0+(L/4), dx3};
p6 = newp; Point(p6) = { R,  0,  0+(L/4), dx3};
p7 = newp; Point(p7) = { 0,  0,  R+(L/4), dx3};
p8 = newp; Point(p8) = {-R,  0,  0+(L/4), dx3};
p9 = newp; Point(p9) = { 0,  0, -R+(L/4), dx3};

// Create the arcs between every two points. Arc is a region so search
// the next id number available using newreg. {start,center,end} 
c5 = newreg; Circle(c5) = {p6, p5, p7};
c6 = newreg; Circle(c6) = {p7, p5, p8};
c7 = newreg; Circle(c7) = {p8, p5, p9};
c8 = newreg; Circle(c8) = {p9, p5, p6};
// Create a new line loop by adding the four arcs with id lcirc=newl
lcirc1 = newl;
Line Loop(lcirc1) = {c5, c6, c7, c8};

// Create the surface between these lines
scyl1=news; // news: next available surface id number
Plane Surface(scyl1) = {lcirc1};

// Extrude the surface
matcyl1 = Extrude {0, L/2-dz, 0} { Surface{scyl1}; }; 





//    Create the volume voxel-cylinders
// For the surface loops consider all the elements of the extrusion EXCEPT the volume. 
// Also add the surface which originated the extrusion to complete 6 faces of the volume

 vox=newsl; // news: find the next available surface loop id number
Surface Loop(vox)={matvox[0],svox,matvox[2],matvox[3],matvox[4],matvox[5]}; // Surface loop from the extrusion for the voxel

cyl=newsl; 
Surface Loop(cyl)={matcyl[0],scyl,matcyl[2],matcyl[3],matcyl[4],matcyl[5]}; // Surface loop from the extrusion for the cylinder1

cyl1=newsl; 
Surface Loop(cyl1)={matcyl1[0],scyl1,matcyl1[2],matcyl1[3],matcyl1[4],matcyl1[5]}; // Surface loop from the extrusion for the cylinder2

volume=newv; // newv: find the next available volume id number
Volume(volume)={vox,cyl,cyl1}; // Volume with boundary vox and hole cyl

// For many cylinder have into account:  Likewise, a surface loop 
// defining a hole should not have any surfaces in common with another 
// surface loop defining a hole in the same volume (in which case the 
// two surface loops should be combined


// Physical Volume
Physical Volume(0) = {volume}; 
