// define geometry-specific parameters
mm = 1.e-3;
gap = 1*mm; //Magnet-Magnet gap [m], gap > 0

DefineConstant[
  cub = {50*mm, Name "Parameters/1Magnet bottom size [m]"}
  hite = {10*mm, Name "Parameters/2Magnet hieght [m]"}
  fdep = {60*mm, Name "Parameters/3Frame depth [m]"}
  inf = {100*mm, Name "Parameters/4Air box distance [m]"}
  lc1 = {TotalMemory <= 2048 ? 30*mm : 10*mm, Name "Parameters/5Mesh size on magnets [m]"}
  lc2 = {TotalMemory <= 2048 ? 90*mm : 20*mm, Name "Parameters/6Mesh size at infinity [m]"}
];

// change global Gmsh options
Mesh.Optimize = 1; // optimize quality of tetrahedra
Mesh.VolumeEdges = 0; // hide volume edges
Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)

//parameters used by this macro
// cx, cy, cz: center of sphere
// rad: radius of sphere
// slc: mesh size

Macro MakeCircleZ //facing z direction
  cp = newp; Point(cp) = {cx, cy, cz, slc}; //center
  dd[] = {rad, 0, -rad, 0};
  For i In {0:3}
    spoint[i] = newp;
	Point(spoint[i]) = {cx+dd[i], cy+dd[3-i], cz, slc};
	If (i>0)
	   quarc[i-1] = newreg; //quarter arcs
	   Circle(quarc[i-1]) = {spoint[i-1], cp, spoint[i]};
	EndIf
  EndFor
  //closing arc
  quarc[3] = newreg; Circle(quarc[3]) = {spoint[3], cp, spoint[0]};
  zloop = newreg; Line Loop(zloop) = {quarc[]};
  zcircle = newreg; Plane Surface(zcircle) = {zloop};
Return

cx = 0; cy=0; cz = -gap-hite;
rad = cub; slc = lc1;
Call MakeCircleZ;
mag1[] = Extrude {0, 0, hite} { Surface{zcircle}; };
jfl1 = zloop; //joint face loop 1

cx = 0; cy=0; cz = gap+hite;
rad = cub; slc = lc1;
Call MakeCircleZ;
mag2[] = Extrude {0, 0, -hite} { Surface{zcircle}; };
jfl2 = zloop; //joint face 2

//Delete {Volume{mag1[1]};}
//Delete {Volume{mag2[1]};}
Physical Volume("Magnets") = {mag1[1], mag2[1]};

//create steel frame around the magnet
p1 = newp; Point(p1) = {-2*cub, -fdep, -hite-gap, lc1};
p2 = newp; Point(p2) = { 2*cub, -fdep, -hite-gap, lc1};
p3 = newp; Point(p3) = { 2*cub, -fdep,  hite+gap, lc1};
p4 = newp; Point(p4) = {-2*cub, -fdep,  hite+gap, lc1};
l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p3};
l3 = newl; Line(l3) = {p3,p4}; l4 = newl; Line(l4) = {p4,p1};

q1 = newp; Point(q1) = {-2*cub, fdep, -hite-gap, lc1};
q2 = newp; Point(q2) = { 2*cub, fdep, -hite-gap, lc1};
q3 = newp; Point(q3) = { 2*cub, fdep,  hite+gap, lc1};
q4 = newp; Point(q4) = {-2*cub, fdep,  hite+gap, lc1};
e1 = newl; Line(e1) = {q1,q2}; e2 = newl; Line(e2) = {q2,q3};
e3 = newl; Line(e3) = {q3,q4}; e4 = newl; Line(e4) = {q4,q1};

d1 = newl; Line(d1) = {p1,q1}; d2 = newl; Line(d2) = {p2,q2};
d3 = newl; Line(d3) = {p3,q3}; d4 = newl; Line(d4) = {p4,q4};

ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};

lex1 = newll; Line Loop(lex1) = {l1, d2, -e1, -d1};
fex1 = news; Plane Surface(fex1) = {lex1}; 
hof1 = news; Plane Surface(hof1) = {lex1, jfl1}; //holed face 1

lex2 = newll; Line Loop(lex2) = {l2, d3, -e2, -d2};
fex2 = news; Plane Surface(fex2) = {lex2}; 

lex3 = newll; Line Loop(lex3) = {l3, d4, -e3, -d3};
fex3 = news; Plane Surface(fex3) = {lex3}; 
hof2 = news; Plane Surface(hof2) = {lex3, jfl2}; //holed face 2

lex4 = newll; Line Loop(lex4) = {l4, d1, -e4, -d4};
fex4 = news; Plane Surface(fex4) = {lex4}; 

ll1d = newll; Line Loop(ll1d) = {e1,e2,e3,e4};

hite2 = hite + gap + 1.5*cub;
p1 = newp; Point(p1) = {-3.5*cub, -fdep, -hite2, lc1};
p2 = newp; Point(p2) = { 3.5*cub, -fdep, -hite2, lc1};
p3 = newp; Point(p3) = { 3.5*cub, -fdep,  hite2, lc1};
p4 = newp; Point(p4) = {-3.5*cub, -fdep,  hite2, lc1};
l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p3};
l3 = newl; Line(l3) = {p3,p4}; l4 = newl; Line(l4) = {p4,p1};

q1 = newp; Point(q1) = {-3.5*cub, fdep, -hite2, lc1};
q2 = newp; Point(q2) = { 3.5*cub, fdep, -hite2, lc1};
q3 = newp; Point(q3) = { 3.5*cub, fdep,  hite2, lc1};
q4 = newp; Point(q4) = {-3.5*cub, fdep,  hite2, lc1};
e1 = newl; Line(e1) = {q1,q2}; e2 = newl; Line(e2) = {q2,q3};
e3 = newl; Line(e3) = {q3,q4}; e4 = newl; Line(e4) = {q4,q1};

d1 = newl; Line(d1) = {p1,q1}; d2 = newl; Line(d2) = {p2,q2};
d3 = newl; Line(d3) = {p3,q3}; d4 = newl; Line(d4) = {p4,q4};

ll2 = newll; Line Loop(ll2) = {l1,l2,l3,l4};

let1 = newll; Line Loop(let1) = {l1, d2, -e1, -d1};
fet1 = news; Plane Surface(fet1) = {let1}; 

let2 = newll; Line Loop(let2) = {l2, d3, -e2, -d2};
fet2 = news; Plane Surface(fet2) = {let2}; 

let3 = newll; Line Loop(let3) = {l3, d4, -e3, -d3};
fet3 = news; Plane Surface(fet3) = {let3}; 

let4 = newll; Line Loop(let4) = {l4, d1, -e4, -d4};
fet4 = news; Plane Surface(fet4) = {let4}; 

ll2d = newll; Line Loop(ll2d) = {e1,e2,e3,e4};

s1 = news; Plane Surface(s1) = {ll2, ll1};
s1d = news; Plane Surface(s1d) = {ll2d, ll1d};

//the frame, comment out the next five lines to delete it
slframe = newsl; 
Surface Loop(slframe)={ s1, fex1, fex2, fex3, fex4,
                        s1d, fet1, fet2, fet3, fet4};
frame = newv; Volume(frame) = {slframe};
Physical Volume("Frame") = {frame};

//now create the surface enclosing both mags and the frame
solids = newsl; 
Surface Loop(solids) = { s1d, fet1, fet2, fet3, fet4, s1,
    hof1, mag1[2], mag1[3], mag1[4], mag1[5], mag1[0], fex2,
    hof2, mag2[2], mag2[3], mag2[4], mag2[5], mag2[0], fex4};

//parameters used by this macro
// cx, cy, cz: center of sphere
// rad: radius of sphere
// slc: mesh size of sphere

Macro MakeSphereSurface
  cp = newp; Point(cp) = {cx, cy, cz, slc}; //center
  np = newp; Point(np) = {cx, cy, cz+rad, slc}; //north
  sp = newp; Point(sp) = {cx, cy, cz-rad, slc}; //south
  dd[] = {rad, 0, -rad, 0};
  For i In {0:3}
    ap = newp;  pe[i] = ap; //points on the equator
	Point(ap) = {cx+dd[i], cy+dd[3-i], cz, slc};
	na[i] = newreg; Circle(na[i]) = {ap, cp, np}; //north arc
	sa[i] = newreg; Circle(sa[i]) = {ap, cp, sp}; //south arc
	If (i>0)
	   eq[i-1] = newreg; //equator
	   Circle(eq[i-1]) = {pe[i-1], cp, ap};
	EndIf
  EndFor
  //closing arc of the equator
  eq[3] = newreg; Circle(eq[3]) = {ap, cp, pe[0]};
  //now the non-plane surfaces: sloop[]
  For i In {0:3}
    lloop = newreg; Line Loop(lloop)={eq[i], -na[i], na[(i+1)%4]};
	sloop[2*i] = newreg; Ruled Surface(sloop[2*i]) = {lloop};
    lloop = newreg; Line Loop(lloop)={eq[i], -sa[i], sa[(i+1)%4]};
    sloop[2*i+1] = newreg; Ruled Surface(sloop[2*i+1]) = {lloop};
  EndFor
  thesurf = newreg;
  Surface Loop(thesurf) = {sloop[]};
Return

// create air box around magnets
BoundingBox; // recompute model bounding box
cx = (General.MinX + General.MaxX) / 2;
cy = (General.MinY + General.MaxY) / 2;
cz = (General.MinZ + General.MaxZ) / 2;
lx = 2*inf + General.MaxX - General.MinX;
ly = 2*inf + General.MaxY - General.MinZ;
lz = 2*inf + General.MaxZ - General.MinZ;
rad = (lx + ly + lz)/3;
slc = lc2;

Call MakeSphereSurface ;

vair = newv; Volume(vair) = {thesurf, solids};
Physical Volume("Air") = {vair}; // air
Physical Surface("Infty") = sloop[];

// adjust value here for correct merge result
//Geometry.Tolerance = 1e-6; 
//Coherence Mesh;