// include parameters common to geometry and solver
Include "magnets_data.pro";

// define geometry-specific parameters
DefineConstant[
  lc1 = {TotalMemory <= 20000 ? 1000*mm : 50*mm, Name "Parameters/2Mesh size on magnets [m]"}
  lc2 = {TotalMemory <= 20000 ? 4000*mm : 5000*mm, Name "Parameters/2Mesh size at infinity [m]"}
  inf = {10000*mm, Name "Parameters/1Air box distance [m]"}
  

  //Density kg/m^3 - Same density for all Grades  
  dens_NdFeB = 7500 	//NdFeB
  dens_Fe = 4900		//FERRITE 
  dens_A = 7300		//ALNICO
];

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

// create magnets
For i In {1:NumMagnets}
  If(M~{i} == 0) // cylinder
    p1 = newp; Point(p1) = {X~{i}, Y~{i}-L~{i}/2, Z~{i}, lc1};
    p2 = newp; Point(p2) = {X~{i}+R~{i}, Y~{i}-L~{i}/2, Z~{i}, lc1};
    p3 = newp; Point(p3) = {X~{i}, Y~{i}+L~{i}/2, Z~{i}, lc1};
    p4 = newp; Point(p4) = {X~{i}+R~{i}, Y~{i}+L~{i}/2, Z~{i}, lc1};
    l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p4};
    l3 = newl; Line(l3) = {p4,p3}; l4 = newl; Line(l4) = {p3,p1};
    ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
    s1 = news; Plane Surface(s1) = {ll1};
    e1[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{s1}; };
    e2[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e1[0]}; };
    e3[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e2[0]}; };
    e4[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e3[0]}; };
    Magnet~{i}[] = {e1[1], e2[1], e3[1], e4[1]};

    Vol~{i}=  Pi*R~{i}*R~{i}*L~{i};
    Printf(" Magnet (%g) Cylinder", i);    
  ElseIf(M~{i} == 1) // paralleleped
    p1 = newp; Point(p1) = {X~{i}-Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
    p2 = newp; Point(p2) = {X~{i}+Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
    p3 = newp; Point(p3) = {X~{i}+Lx~{i}/2, Y~{i}+Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
    p4 = newp; Point(p4) = {X~{i}-Lx~{i}/2, Y~{i}+Ly~{i}/2, Z~{i}-Lz~{i}/2, 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};
    e1[] = Extrude {0, 0, Lz~{i}} { Surface{s1}; };
    Magnet~{i}[] = {e1[1]};
    
    Vol~{i}= Lx~{i}*Ly~{i}*Lz~{i} ;
    Printf(" Magnet (%g) Cube ", i);
  ElseIf(M~{i} == 2) // tube
    p1 = newp; Point(p1) = {X~{i}+Rin~{i}, Y~{i}-L~{i}/2, Z~{i}, lc1};
    p2 = newp; Point(p2) = {X~{i}+R~{i}, Y~{i}-L~{i}/2, Z~{i}, lc1};
    p3 = newp; Point(p3) = {X~{i}+Rin~{i}, Y~{i}+L~{i}/2, Z~{i}, lc1};
    p4 = newp; Point(p4) = {X~{i}+R~{i}, Y~{i}+L~{i}/2, Z~{i}, lc1};
    l1 = newl; Line(l1) = {p1,p2}; l2 = newl; Line(l2) = {p2,p4};
    l3 = newl; Line(l3) = {p4,p3}; l4 = newl; Line(l4) = {p3,p1};
    ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
    s1 = news; Plane Surface(s1) = {ll1};
    e1[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{s1}; };
    e2[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e1[0]}; };
    e3[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e2[0]}; };
    e4[] = Extrude {{0, 1, 0}, {X~{i}, Y~{i}, Z~{i}}, Pi/2} { Surface{e3[0]}; };
    Magnet~{i}[] = {e1[1], e2[1], e3[1], e4[1]};
    
    Vol~{i}= L~{i}*Pi*( (R~{i}*R~{i})-(Rin~{i}*Rin~{i}) );   
    If(Vol~{i} >= 0)
	Printf(" Magnet (%g) Tube ", i);
    Else
	Printf(" Magnet (%g) Warning: The inside values cannot be larger than the outside (Tube)",i);
        Printf("   ");
    EndIf
  Else // parallelepiped with hole (Ly = Liny)
    p1 = newp; Point(p1) = {X~{i}-Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
    p2 = newp; Point(p2) = {X~{i}+Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Lz~{i}/2, lc1};
    p3 = newp; Point(p3) = {X~{i}+Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}+Lz~{i}/2, lc1};
    p4 = newp; Point(p4) = {X~{i}-Lx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}+Lz~{i}/2, 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};
    
    p5 = newp; Point(p5) = {X~{i}-Linx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Linz~{i}/2, lc1};
    p6 = newp; Point(p6) = {X~{i}+Linx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}-Linz~{i}/2, lc1};
    p7 = newp; Point(p7) = {X~{i}+Linx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}+Linz~{i}/2, lc1};
    p8 = newp; Point(p8) = {X~{i}-Linx~{i}/2, Y~{i}-Ly~{i}/2, Z~{i}+Linz~{i}/2, lc1};
    l5 = newl; Line(l5) = {p5,p6}; l6 = newl; Line(l6) = {p6,p7};
    l7 = newl; Line(l7) = {p7,p8}; l8 = newl; Line(l8) = {p8,p5};
    ll2 = newll; Line Loop(ll2) = {l5,l6,l7,l8};
    
    s1 = news; Plane Surface(s1) = {ll1,ll2};
    e1[] = Extrude {0, Ly~{i}, 0} { Surface{s1}; };
    Magnet~{i}[] = {e1[1]};
    
    VolOut~{i}= Lx~{i}*Ly~{i}*Lz~{i};
    VolIn~{i} = Linx~{i}*Ly~{i}*Linz~{i};
    Vol~{i} = VolOut~{i} - VolIn~{i};
    If(Vol~{i} >= 0)
	Printf(" Magnet (%g) SquareTube ", i);	
    Else
	Printf(" Magnet (%g) Warning: The inside values cannot be larger than the outside (SquareTube)", i);
	Printf("   ");
    EndIf
  EndIf
    
  Rotate { {0,0,1}, {X~{i},Y~{i},Z~{i}}, Rz~{i}*Pi/180  } { Volume{Magnet~{i}[]}; }
  Rotate { {0,1,0}, {X~{i},Y~{i},Z~{i}}, Ry~{i}*Pi/180  } { Volume{Magnet~{i}[]}; }
  Rotate { {1,0,0}, {X~{i},Y~{i},Z~{i}}, Rx~{i}*Pi/180  } { Volume{Magnet~{i}[]}; }

  Physical Volume(i) = Magnet~{i}[]; // magnet volume
  skin~{i}[] = CombinedBoundary{ Volume{Magnet~{i}[]}; };
  Physical Surface(100+i) = skin~{i}[]; // magnet skin
  
	
//Calculate mass from Volume and Material Density 
  If(Vol~{i} >= 0)
    Printf(" Magnet (%g) Volume = %g  m^3", i, Vol~{i});
	If(Mat~{i} == 0) //NDFeB
		Mass~{i}=Vol~{i}*dens_NdFeB;
		Printf(" Magnet (%g) Mass  = %g kg (NdFeB)", i, Mass~{i});
	ElseIf(Mat~{i} == 1) //Ferrite C5
		Mass~{i}=Vol~{i}*dens_Fe;
		Printf(" Magnet (%g) Mass = %g kg (Ferrite)", i, Mass~{i});
	Else(Mat~{i} == 2) //Alnico 5
		Mass~{i}=Vol~{i}* dens_A ;
		Printf(" Magnet (%g) Mass = %g kg (Alnico) ", i, Mass~{i});
    EndIf
  EndIf	
	
 Printf("   ");
EndFor

  
// 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;
p1 = newp; Point (p1) = {cx-lx/2, cy-ly/2, cz-lz/2, lc2};
p2 = newp; Point (p2) = {cx+lx/2, cy-ly/2, cz-lz/2, lc2};
l1 = newl; Line(l1) = {p1, p2};
e1[] = Extrude {0, ly, 0} { Line{l1}; };
e2[] = Extrude {0, 0, lz} { Surface{e1[1]}; };
Delete { Volume{e2[1]}; }
ss[] = {e1[1],e2[0],e2[2],e2[3],e2[4],e2[5]};
sl1 = newsl; Surface Loop(sl1) = {ss[]};
vv[] = {sl1};

// add magnets surface boundary
For i In {1:NumMagnets}
  sl~{i} = newsl; Surface Loop(sl~{i}) = skin~{i}[];
  vv[] += sl~{i};
EndFor
v1 = newv; Volume(v1) = {vv[]};

//create particles
For i In {0:NumParticles-1}
  For j In {0:NumParticles-1}  
    p1 = newp; Point(p1) = {-inf/2+(i)*inf/(NumParticles-1),inf,-inf/2+(j)*inf/(NumParticles-1),lc1};
    Physical Point(200+(i)*NumParticles+j) = {p1}; // Define Particles as physical point
    Point{p1} In Volume {v1};// add particles to volume so that they are included in mesh
  EndFor
EndFor

Physical Volume(NumMagnets+1) = v1; // air
Physical Surface(NumMagnets+2) = ss[]; // infinity
