//========================================================
// Benchmark "3D unit cell"  
// File: Geometry
// Contributors:
//    Miguel Rojas S.
//    Pablo Stuardo O.
//========================================================

Include "unitcell.dat";

// Making the inner cylinder
p[] += newp; Point(newp) = {0,0,0,GZ};
p[] += newp; Point(newp) = {0,0,H,GZ};

p[] += newp; Point(newp) = {x1,-y1,0,GZ};
p[] += newp; Point(newp) = {x1,-y1,H,GZ};

p[] += newp; Point(newp) = {x1,y1,0,GZ};
p[] += newp; Point(newp) = {x1,y1,H,GZ};

// Making the outer cylinder
p[] += newp; Point(newp) = {x2,-y2,0,GZ};
p[] += newp; Point(newp) = {x2,-y2,H,GZ};

p[] += newp; Point(newp) = {x2,y2,0,GZ};
p[] += newp; Point(newp) = {x2,y2,H,GZ};

l[] += newl; Line(newl) = {p[0],p[1]};
l[] += newl; Line(newl) = {p[1],p[3]};
l[] += newl; Line(newl) = {p[3],p[2]};
l[] += newl; Line(newl) = {p[2],p[0]};
l[] += newl; Line(newl) = {p[3],p[7]};
l[] += newl; Line(newl) = {p[7],p[6]};
l[] += newl; Line(newl) = {p[6],p[2]};

l[] += newl; Line(newl) = {p[1],p[5]};
l[] += newl; Line(newl) = {p[5],p[4]};
l[] += newl; Line(newl) = {p[4],p[0]};
l[] += newl; Line(newl) = {p[5],p[9]};
l[] += newl; Line(newl) = {p[9],p[8]};
l[] += newl; Line(newl) = {p[8],p[4]};

l[] += newl; Circle(newl) = {p[4],p[0],p[2]};
l[] += newl; Circle(newl) = {p[6],p[0],p[8]};
l[] += newl; Circle(newl) = {p[5],p[1],p[3]};
l[] += newl; Circle(newl) = {p[7],p[1],p[9]};

// Declaring surfaces (inner cylinder)
ll[] += newll; Line loop(newll) = {l[0],l[1],l[2],l[3]};
s[] += news; Surface(news) = {-ll[0]};

ll[] += newll; Line loop(newll) = {l[13],l[3],-l[9]};
s[] += news; Surface(news) = {ll[1]};

ll[] += newll; Line loop(newll) = {l[0],l[7],l[8],l[9]};
s[] += news; Surface(news) = {ll[2]};

ll[] += newll; Line loop(newll) = {l[7],l[15],-l[1]};
s[] += news; Surface(news) = {-ll[3]};

ll[] += newll; Line loop(newll) = {l[15],l[2],-l[13],-l[8]};
s[] += news; Surface(news) = {-ll[4]};							//s27

// Declaring surfaces (outer cylinder) 
ll[] += newll; Line loop(newll) = {l[10],l[11],l[12],-l[8]};    //s29
s[] += news; Surface(news) = {ll[5]};

ll[] += newll; Line loop(newll) = {l[14],l[12],l[13],-l[6]};	//s31
s[] += news; Surface(news) = {-ll[6]};

ll[] += newll; Line loop(newll) = {l[4],l[5],l[6],-l[2]};		//s33
s[] += news; Surface(news) = {-ll[7]};

ll[] += newll; Line loop(newll) = {l[15],l[4],l[16],-l[10]};	//s35
s[] += news; Surface(news) = {ll[8]};

ll[] += newll; Line loop(newll) = {l[5],l[14],-l[11],-l[16]};	//s37
s[] += news; Surface(news) = {ll[9]};


// Declaring volumes (inner cylinder)
sl[] += newsl; Surface loop(newsl) = {s[0],s[1],s[2],s[3],s[4]};
v[] += newv; Volume(newv) = {sl[0]};

// Declaring volumes (outer cylinder)
sl[] += newsl; Surface loop(newsl) = {s[4],s[5],s[6],s[7],s[8],s[9]};
v[] += newv; Volume(newv) = {sl[1]};


// Declaring transfinites surfaces
Transfinite Surface {s[0]} = {p[0],p[1],p[2],p[3]} Right;
Transfinite Surface {s[2]} = {p[0],p[1],p[4],p[5]} Right;

Transfinite Surface {s[1]} = {p[0],p[2],p[4]} Right;
Transfinite Surface {s[3]} = {p[1],p[3],p[5]} Right;

Transfinite Surface {s[5]} = {p[4],p[5],p[8],p[9]} Right;
Transfinite Surface {s[7]} = {p[2],p[3],p[6],p[7]} Right;

Transfinite Surface {s[9]} = {p[6],p[7],p[8],p[9]} Right;
Transfinite Surface {s[4]} = {p[2],p[3],p[4],p[5]} Right;


// Physical declaration
Physical Surface(MasterA) = {s[1],s[6]};  
Physical Surface(SlaveA) = {s[3],s[8]};
Physical Surface(MasterB) = {s[0],s[7]};
Physical Surface(SlaveB) = {s[2],s[5]};
Physical Surface(PEC) = {s[9]};
Physical Volume(Vol) = {v[]};

