// define geometry-specific parameters
mm = 1.e-3;
gap = 1*mm; //Frame-Magnet gap [m]
hasFrame = (gap >= 0.1*mm);

DefineConstant[
  cub = {10*mm, Name "Parameters/2Magnet bottom size [m]"}
  hite = {20*mm, Name "Parameters/2Magnet hieght [m]"}
  lc1 = {TotalMemory <= 2048 ? 8*mm : 4*mm, Name "Parameters/3Mesh size on magnets [m]"}
  lc2 = {TotalMemory <= 2048 ? 40*mm : 20*mm, Name "Parameters/4Mesh size at infinity [m]"}
  inf = {100*mm, Name "Parameters/1Air box distance [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)

p1 = newp; Point(p1) = {-cub, -cub, -hite, lc1};
p2 = newp; Point(p2) = { cub, -cub, -hite, lc1};
p3 = newp; Point(p3) = { cub,  cub, -hite, lc1};
p4 = newp; Point(p4) = {-cub,  cub, -hite, 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};
ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
s1 = news; Plane Surface(s1) = {ll1};
mag[] = Extrude {0, 0, 2*hite} { Surface{s1}; };
Physical Volume("Magnet") = {mag[1]};


If (hasFrame)
	//create steel frame around the magnet
	p1 = newp; Point(p1) = {-2*cub, -cub, -hite-gap, lc1};
	p2 = newp; Point(p2) = { 2*cub, -cub, -hite-gap, lc1};
	p3 = newp; Point(p3) = { 2*cub, -cub,  hite+gap, lc1};
	p4 = newp; Point(p4) = {-2*cub, -cub,  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};
	ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
	hite2 = hite + cub + gap;
	p1 = newp; Point(p1) = {-4*cub, -cub, -hite2, lc1};
	p2 = newp; Point(p2) = { 4*cub, -cub, -hite2, lc1};
	p3 = newp; Point(p3) = { 4*cub, -cub,  hite2, lc1};
	p4 = newp; Point(p4) = {-4*cub, -cub,  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};
	ll2 = newll; Line Loop(ll2) = {l1,l2,l3,l4};
	s1 = news; Plane Surface(s1) = {ll2, ll1};
	frame[] = Extrude {0, 2*cub, 0} { Surface{s1}; };
	Physical Volume("Frame") = {frame[1]};
EndIf

//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 ;

solids = newsl; 
If (hasFrame)
   Surface Loop(solids) = CombinedBoundary{ Volume{mag[1], frame[1]}; };
Else
   Surface Loop(solids) = Boundary{ Volume{mag[1]}; };
EndIf
vair = newv; Volume(vair) = {thesurf, solids};
Physical Volume("Air") = {vair}; // air
Physical Surface("Infty") = sloop[];
