Function MyCircle
  // Given a zshift and two radii rad and xshift, and a zshift, this
  // creates a torus parallel to the x-y-plane.
  // The points:
  tp1 = newp;
  Point(tp1) = {xshift,yshift,z,lcar};
  tp2 = newp;
  Point(tp2) = {rad+xshift,yshift,z,lcar};
  tp3 = newp;
  Point(tp3) = {xshift,yshift+rad,z,lcar};
  tp4 = newp;
  Point(tp4) = {xshift,yshift-rad,z,lcar};
  tp5 = newp;
  Point(tp5) = {-rad+xshift,yshift,z,lcar};

  // Build circle from arcs.
  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};

  theloops[t] = newreg;
  Line Loop(theloops[t]) = {tc1,tc2,tc3,tc4};
  surf1 = newreg;
  Plane Surface(surf1) = {theloops[t]};

  // Name the circles.
  Physical Surface(Sprintf("coil %g", t)) = surf1;

  // TODO don't hardcode field names
  Field[3] = Attractor;
  Field[3].NNodesByEdge = 100;
  Field[3].EdgesList = {tc1, tc2, tc3, tc4};
  Field[4] = Threshold;                                                           
  Field[4].IField = 3;
  Field[4].LcMin = lcar_fine;
  Field[4].LcMax = lcar;
  Field[4].DistMin = width;
  Field[4].DistMax = 2 * width;
Return

z = 0;
lcar = 0.2;
lcar_fine = lcar / 10.0;
width = 0.1;

// Circle.
t = 1;
rad = 1.0;
xshift = 0;
yshift = 0.0;
Call MyCircle;

z = 0;
lcar_air = lcar;
t = 2;
xmin = -2;
xmax = 2.0;
ymin = -2.0;
ymax= 2.0;

// Points.
cp1 = newp;
Point(cp1) = {xmin,ymin,z,lcar_air};
cp2 = newp;
Point(cp2) = {xmax,ymin,z,lcar_air};
cp3 = newp;
Point(cp3) = {xmax,ymax,z,lcar_air};
cp4 = newp;
Point(cp4) = {xmin,ymax,z,lcar_air};

// Lines.
cl1 = newreg;
Line(cl1) = {cp1,cp2};
cl2 = newreg;
Line(cl2) = {cp2,cp3};
cl3 = newreg;
Line(cl3) = {cp3,cp4};
cl4 = newreg;
Line(cl4) = {cp4,cp1};

// Stitch it all together.
t = 0;
theloops[t] = newreg;
Line Loop(theloops[t]) = {cl1,cl2,cl3,cl4};
air = newreg;
Plane Surface(air) = {theloops[]};
Physical Surface(Sprintf("air")) = air;

// Make elements match at interface boundaries.
// Without 'Coherence', the elements may still match, but the
// edges may be registered twice. This causes trouble with
// Extrude.
Coherence;

Field[7] = Min;                                                                 
Field[7].FieldsList = {4};
Background Field = 7;
