

MAGNET = 100;
AIR = 101;
AIRINF = 102;

// define variables for geometry parameters
DefineConstant[ H = {5, Min 0.01, Max 30, Step 1, // height of PM
Name "Parameters/Geometry/Magnet height (mm)"} ] ;

DefineConstant[ Ri = {2, Min 0.01, Max 10, Step 1, // inner radius
Name "Parameters/Geometry/Inner radius (mm)"} ] ;

DefineConstant[ Ro = {5, Min 0.01, Max 10, Step 1, // outer radius
Name "Parameters/Geometry/Outer radius (mm)"} ] ;

DefineConstant[ Air_r = {10, Min 10, Max 100, Step 1,
Name "Parameters/Geometry/Radius of Air (mm)"} ] ;

lm = 2; // resolution of magnet
lr = 5; // resolution of air region

// one segment
po = newp; Point(po) = {Ri, 0, -H/2, lm};
Point(po+1) = {Ro, 0, -H/2, lm};
Point(po+2) = {Ro, 0, H/2, lm};
Point(po+3) = {Ri, 0, H/2, lm};
l = newl; Line(l) = {po, po+1};
Line(l+1) = {po+1,po+2};
Line(1+2) = {po+2,po+3};
Line(l+3) = {po+3,po};
LLSo = newreg; Line Loop(LLSo) = {l,l+1,l+2,l+3};
So = newreg; Plane Surface(So) = {LLSo};
ext_mag = Extrude{{0,0,1},{0,0,0},Pi/2}{Surface{So};};
SL1 = newreg; Surface Loop(SL1) = {ext_mag[0], ext_mag[2], So, ext_mag[3], ext_mag[4], ext_mag[5]};


// Air
p = newp; Point(p) = {  0 ,    0, 0,  lr/2};
Point(p+1) = {  Air_r,    0,  0,  lr};
Point(p+2) = {  0 ,   Air_r,  0,  lr};
c = newc; Circle(c) = {p+1,p,p+2};
Point(p+3) = { -Air_r,    0,  0,  lr};
Point(p+4) = {   0,  -Air_r,  0,  lr};
Circle(c+1) = {p+2,p,p+3};
Circle(c+2) = {p+3,p,p+4};
Circle(c+3) = {p+4,p,p+1};
Point(p+5) = {   0,    0,  Air_r, lr};
Point(p+6) = {   0,    0,  -Air_r, lr};
Circle(c+4) = {p+2,p,p+5};
Circle(c+5) = {p+5,p,p+4};
Circle(c+6) = {p+4,p,p+6};
Circle(c+7) = {p+6,p,p+2};
Circle(c+8) = {p+1,p,p+6};
Circle(c+9) = {p+6,p,p+3};
Circle(c+10) = {p+3,p,p+5};
Circle(c+11) = {p+5,p,p+1};
LLA1 = newll; Line Loop(LLA1) = {c+1,c+7,-(c+9)};
SA1 = news; Ruled Surface(SA1) = {LLA1};
LLA2 = newll; Line Loop(LLA2) = {c+9,c+2,c+6};
SA2 = news; Ruled Surface(SA2) = {LLA2};
LLA3 = newll; Line Loop(LLA3) = {-(c+7),-(c+8),c};
SA3 = news; Ruled Surface(SA3) = {LLA3};
LLA4 = newll; Line Loop(LLA4) = {-(c+10),-(c+1),c+4};
SA4 = news; Ruled Surface(SA4) = {LLA4};
LLA5 = newll; Line Loop(LLA5) = {-(c+4),-(c+11),-c};
SA5 = news; Ruled Surface(SA5) = {LLA5};
LLA6 = newll; Line Loop(LLA6) = {-(c+2),c+10,c+5};
SA6 = news; Ruled Surface(SA6) = {LLA6};
LLA7 = newll; Line Loop(LLA7) = {-(c+6),c+3,c+8};
SA7 = news; Ruled Surface(SA7) = {LLA7};
LLA8 = newll; Line Loop(LLA8) = {-(c+3),c+11,-(c+5)};
SA8 = news; Ruled Surface(SA8) = {LLA8};
SL2 = newsl; Surface Loop(SL2) = {SA1,SA2,SA3,SA4,SA5,SA6,SA7,SA8};
V1 = newreg; Volume(V1) = {SL2, SL1};

// Physical entities
Physical Volume(MAGNET) = {ext_mag[1]};
Physical Volume(AIR) = {V1};
Physical Surface(AIRINF) = {SL2};

Mesh.Algorithm=2;
Mesh.Algorithm3D=1;

Coherence;
