// Gmsh geometry for a rectangular coil surrounded by by two tori.

Function Brick
  // Creates a rectangular brick.

  // Define points.
  cp1 = newp;
  Point(cp1) = {xmax,ymax,zmax,lcar};
  cp2 = newp;
  Point(cp2) = {xmax,ymax,zmin,lcar};
  cp3 = newp;
  Point(cp3) = {xmax,ymin,zmax,lcar};
  cp4 = newp;
  Point(cp4) = {xmax,ymin,zmin,lcar};
  cp5 = newp;
  Point(cp5) = {xmin,ymax,zmax,lcar};
  cp6 = newp;
  Point(cp6) = {xmin,ymax,zmin,lcar};
  cp7 = newp;
  Point(cp7) = {xmin,ymin,zmax,lcar};
  cp8 = newp;
  Point(cp8) = {xmin,ymin,zmin,lcar};

  // Lines.
  cl1 = newreg;
  Line(cl1) = {cp1,cp2};
  cl2 = newreg;
  Line(cl2) = {cp1,cp3};
  cl3 = newreg;
  Line(cl3) = {cp1,cp5};
  cl4 = newreg;
  Line(cl4) = {cp2,cp4};
  cl5 = newreg;
  Line(cl5) = {cp2,cp6};
  cl6 = newreg;
  Line(cl6) = {cp3,cp4};
  cl7 = newreg;
  Line(cl7) = {cp3,cp7};
  cl8 = newreg;
  Line(cl8) = {cp4,cp8};
  cl9 = newreg;
  Line(cl9) = {cp5,cp6};
  cl10 = newreg;
  Line(cl10) = {cp5,cp7};
  cl11 = newreg;
  Line(cl11) = {cp6,cp8};
  cl12 = newreg;
  Line(cl12) = {cp7,cp8};

  // Surfaces.
  cll1 = newreg; Line Loop(cll1) = {cl1,cl4,-cl6,-cl2};    Plane Surface(newreg) = {cll1};
  cll2 = newreg; Line Loop(cll2) = {cl1,cl5,-cl9,-cl3};    Plane Surface(newreg) = {cll2};
  cll3 = newreg; Line Loop(cll3) = {cl2,cl7,-cl10,-cl3};   Plane Surface(newreg) = {cll3};
  cll4 = newreg; Line Loop(cll4) = {cl4,cl8,-cl11,-cl5};   Plane Surface(newreg) = {cll4};
  cll5 = newreg; Line Loop(cll5) = {cl6,cl8,-cl12,-cl7};   Plane Surface(newreg) = {cll5};
  cll6 = newreg; Line Loop(cll6) = {cl9,cl11,-cl12,-cl10}; Plane Surface(newreg) = {cll6};

  // We then store the surface loops identification numbers in a list
  // for later reference (we will need these to define the final
  // volume).
  theloops[t] = newreg;
  // Define outer surface.
  Surface Loop(theloops[t]) = {cll1+1,cll2+1,cll3+1,cll4+1,cll5+1,cll6+1};

Return

Function Torus
  // Given a zshift and two radii irad and orad, and a zshift, this
  // creates a torus parallel to the x-y-plane.
  // The points:
  tp1 = newp;
  Point(tp1) = {0,orad,zshift,lcar};
  tp2 = newp;
  Point(tp2) = {0,irad+orad,zshift,lcar};
  tp3 = newp;
  Point(tp3) = {0,orad,zshift+irad,lcar};
  tp4 = newp;
  Point(tp4) = {0,orad,zshift-irad,lcar};
  tp5 = newp;
  Point(tp5) = {0,-irad+orad,zshift,lcar};
  // One cirlce:
  tc1 = newreg;
  Circle(tc1) = {tp2,tp1,tp3};
  tc2 = newreg;
  Circle(tc2) = {tp3,tp1,tp5};
  tc3 = newreg;
  Circle(tc3) = {tp5,tp1,tp4};
  tc4 = newreg;
  Circle(tc4) = {tp4,tp1,tp2};
  // The extrusion to the torus:
  tll1 = newreg;
  Line Loop(tll1) = {tc1,tc2,tc3,tc4};
  ts1 = newreg;
  Plane Surface(ts1) = {tll1};
  // Gmsh cannot rotate beyond PI, so split the extrusion up in three parts.
  // numLayers = 2*Pi*orad / lcar
  //numLayers = 2*Pi*orad / lcar;
  Printf("numLayers = %g", numLayers);
  ts2[] = Extrude {{0,0,1}, {0,0,0}, 2*Pi/3}{Surface{ts1}; Layers{numLayers};};
  ts3[] = Extrude {{0,0,1}, {0,0,0}, 2*Pi/3}{Surface{ts2[0]}; Layers{numLayers};};
  ts4[] = Extrude {{0,0,1}, {0,0,0}, 2*Pi/3}{Surface{ts3[0]}; Layers{numLayers};};

  // We then store the surface loops identification numbers in a list
  // for later reference (we will need these to define the final
  // volume).
  theloops[t] = newreg;
  // Define outer surface.
  //Surface Loop(theloops[t]) = {ts1, ts2[0], ts2[2], ts2[3], ts2[4], ts2[5]};
  Surface Loop(theloops[t]) = {ts2[2], ts2[3], ts2[4], ts2[5],
                               ts3[2], ts3[3], ts3[4], ts3[5],
                               ts4[2], ts4[3], ts4[4], ts4[5]};

  thetorus = newreg;
  Volume(thetorus) = theloops[t];
  Physical Volume (t) = thetorus;
Return

// Actual creation starts here.

//// Coil.
//t = 1;
//lcar = 0.01;
//xmin = -0.05;
//xmax = 0.05;
//ymin = -0.05;
//ymax = 0.05;
//zmin = -0.15;
//zmax = 0.15;
//Call Brick;
//thebrick = newreg;
//Volume(thebrick) = theloops[t];
//Physical Volume (t) = thebrick;


// Tori.
t = 1;
lcar = 0.05;
orad = 0.5;
irad = 0.2;
zshift = 0.04;
// 2*pi*orad / lcar
numLayers = 63;
Call Torus;
//t = 2;
//zshift = -0.04;
//Call Torus;

// Outer cube.
t = 0;
lcar = 0.5;
xmin = -0.75;
xmax = 0.75;
ymin = -0.75;
ymax = 0.75;
zmin = -0.75;
zmax = 0.75;
Call Brick;
thebrick = newreg;
// The volume of the cube, without the cavities, is now defined by a number of
// surface loops: the first surface loop (t=0) defines the exterior surface;
// the surface loops other than the first one define holes.  (Again,
// to reference an array of variables, its identifier is followed by
// square brackets):
Volume(thebrick) = {theloops[]};
// We finally define a physical volume for the elements discretizing
// the cube, without the holes (whose elements were already tagged
// with numbers 1 to 5 in the `For' loop):
Physical Volume (newreg) = thebrick;
Physical Volume(58) = {1};
Physical Volume(59) = {1};
