Include "parameters.pro";

nm = 1.e-9;		/* conversion factor from meters to nanometers */
//bigmeshsiz = 30*nm;
//smallmeshsiz = 10*nm;
domainZsize = 1000*nm;
domainYsize = 5000*nm;
domainXsize = 1000*nm;
bigmeshsiz = domainXsize/5;
smallmeshsiz = 10*nm;

/* Front left */
frontleftoffset = 1;
xmin = -49.*nm;
xmax = -25.*nm;
zmin = -0.5*domainZsize/2;
zmax = 0.5*domainZsize/2;
y1 = 0.*nm;
y2 = 89.*nm;
y3 = 90.*nm;
y4 = 99.*nm;
y5 = 100.*nm;

nextpointnum = frontleftoffset;
Point(nextpointnum) = {xmin,y1,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y1,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y2,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y2,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y3,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y3,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y4,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y4,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y5,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y5,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;

backleftoffset = frontleftoffset+19;
nextpointnum = backleftoffset;
Point(nextpointnum) = {xmin,y1,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y1,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y2,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y2,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y3,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y3,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y4,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y4,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y5,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y5,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;

leftoffset = 1;
nextline = leftoffset;
nextarea = leftoffset;
nextvolume = leftoffset;
firstfrontpoint = frontleftoffset;
firstbackpoint = backleftoffset;
// Lines and areas on front face
Line (nextline) = {firstfrontpoint,firstfrontpoint+1}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+1,firstfrontpoint+3}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+3,firstfrontpoint+2}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+2,firstfrontpoint}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+3,firstfrontpoint+5}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+5,firstfrontpoint+4}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+4,firstfrontpoint+2}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+5,firstfrontpoint+7}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+7,firstfrontpoint+6}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+6,firstfrontpoint+4}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+7,firstfrontpoint+9}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+9,firstfrontpoint+8}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+8,firstfrontpoint+6}; nextline = nextline+1;
Line Loop (nextarea) = {leftoffset,leftoffset+1,leftoffset+2,leftoffset+3};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+2),leftoffset+4,leftoffset+5,leftoffset+6};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+5),leftoffset+7,leftoffset+8,leftoffset+9};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+8),leftoffset+10,leftoffset+11,leftoffset+12};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
//Lines and areas on back face
firstbackline = nextline;
Line (nextline) = {firstbackpoint,firstbackpoint+1}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+1,firstbackpoint+3}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+3,firstbackpoint+2}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+2,firstbackpoint}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+3,firstbackpoint+5}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+5,firstbackpoint+4}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+4,firstbackpoint+2}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+5,firstbackpoint+7}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+7,firstbackpoint+6}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+6,firstbackpoint+4}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+7,firstbackpoint+9}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+9,firstbackpoint+8}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+8,firstbackpoint+6}; nextline = nextline+1;
Line Loop (nextarea) = {firstbackline,firstbackline+1,firstbackline+2,firstbackline+3};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(firstbackline+2),firstbackline+4,firstbackline+5,firstbackline+6};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(firstbackline+5),firstbackline+7,firstbackline+8,firstbackline+9};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(firstbackline+8),firstbackline+10,firstbackline+11,firstbackline+12};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
//Lines connecting front to back
firstsideline = nextline;
Line (nextline) = {frontleftoffset,backleftoffset}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+1,backleftoffset+1}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+2,backleftoffset+2}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+3,backleftoffset+3}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+4,backleftoffset+4}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+5,backleftoffset+5}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+6,backleftoffset+6}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+7,backleftoffset+7}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+8,backleftoffset+8}; nextline = nextline+1;
Line (nextline) = {frontleftoffset+9,backleftoffset+9}; nextline = nextline+1;
//Horizontal areas, counterclockwise from above
Line Loop (nextarea) = {leftoffset,firstsideline+1,-firstbackline,-firstsideline}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+2),firstsideline+3,firstbackline+2,-(firstsideline+2)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+5),firstsideline+5,firstbackline+5,-(firstsideline+4)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+8),firstsideline+7,firstbackline+8,-(firstsideline+6)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+11),firstsideline+9,firstbackline+11,-(firstsideline+8)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
// Left side areas, outward pointing
Line Loop (nextarea) = {-(leftoffset+3),firstsideline+2,firstbackline+3,-firstsideline}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+6),firstsideline+4,firstbackline+6,-(firstsideline+2)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+9),firstsideline+6,firstbackline+9,-(firstsideline+4)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+12),firstsideline+8,firstbackline+12,-(firstsideline+6)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
//Right side areas, outward pointing
Line Loop (nextarea) = {-(leftoffset+1),firstsideline+1,firstbackline+1,-(firstsideline+3)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+4),firstsideline+3,firstbackline+4,-(firstsideline+5)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+7),firstsideline+5,firstbackline+7,-(firstsideline+7)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(leftoffset+10),firstsideline+7,firstbackline+10,-(firstsideline+9)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;

//Left line volumes, starting at the bottom, in order
// {front,bottom,right,back,top,left}
Surface Loop (nextvolume) = {leftoffset,-(leftoffset+8),leftoffset+17,-(leftoffset+4),leftoffset+9,leftoffset+13};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;
Surface Loop (nextvolume) = {leftoffset+1,-(leftoffset+9),leftoffset+18,-(leftoffset+5),leftoffset+10,leftoffset+14};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;
Surface Loop (nextvolume) = {leftoffset+2,-(leftoffset+10),leftoffset+19,-(leftoffset+6),leftoffset+11,leftoffset+15};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;
Surface Loop (nextvolume) = {leftoffset+3,-(leftoffset+11),leftoffset+20,-(leftoffset+7),leftoffset+12,leftoffset+16};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;


/* Right */
frontrightoffset = 100;
xmin = 25.*nm;
xmax = 49.*nm;

nextpointnum = frontrightoffset;
Point(nextpointnum) = {xmin,y1,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y1,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y2,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y2,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y3,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y3,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y4,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y4,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y5,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y5,zmax,smallmeshsiz}; nextpointnum = nextpointnum+1;

backrightoffset = frontrightoffset+19;
nextpointnum = backrightoffset;
Point(nextpointnum) = {xmin,y1,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y1,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y2,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y2,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y3,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y3,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y4,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y4,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {xmin,y5,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1; 
Point(nextpointnum) = {xmax,y5,zmin,smallmeshsiz}; nextpointnum = nextpointnum+1;

rightoffset = 100;
nextline = rightoffset;
nextarea = rightoffset;
nextvolume = rightoffset;
firstfrontpoint = frontrightoffset;
firstbackpoint = backrightoffset;
// Lines and areas on front face
Line (nextline) = {firstfrontpoint,firstfrontpoint+1}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+1,firstfrontpoint+3}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+3,firstfrontpoint+2}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+2,firstfrontpoint}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+3,firstfrontpoint+5}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+5,firstfrontpoint+4}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+4,firstfrontpoint+2}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+5,firstfrontpoint+7}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+7,firstfrontpoint+6}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+6,firstfrontpoint+4}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+7,firstfrontpoint+9}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+9,firstfrontpoint+8}; nextline = nextline+1;
Line (nextline) = {firstfrontpoint+8,firstfrontpoint+6}; nextline = nextline+1;
Line Loop (nextarea) = {rightoffset,rightoffset+1,rightoffset+2,rightoffset+3};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+2),rightoffset+4,rightoffset+5,rightoffset+6};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+5),rightoffset+7,rightoffset+8,rightoffset+9};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+8),rightoffset+10,rightoffset+11,rightoffset+12};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
//Lines and areas on back face
firstbackline = nextline;
Line (nextline) = {firstbackpoint,firstbackpoint+1}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+1,firstbackpoint+3}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+3,firstbackpoint+2}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+2,firstbackpoint}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+3,firstbackpoint+5}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+5,firstbackpoint+4}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+4,firstbackpoint+2}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+5,firstbackpoint+7}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+7,firstbackpoint+6}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+6,firstbackpoint+4}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+7,firstbackpoint+9}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+9,firstbackpoint+8}; nextline = nextline+1;
Line (nextline) = {firstbackpoint+8,firstbackpoint+6}; nextline = nextline+1;
Line Loop (nextarea) = {firstbackline,firstbackline+1,firstbackline+2,firstbackline+3};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(firstbackline+2),firstbackline+4,firstbackline+5,firstbackline+6};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(firstbackline+5),firstbackline+7,firstbackline+8,firstbackline+9};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(firstbackline+8),firstbackline+10,firstbackline+11,firstbackline+12};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
//Lines connecting front to back
firstsideline = nextline;
Line (nextline) = {frontrightoffset,backrightoffset}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+1,backrightoffset+1}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+2,backrightoffset+2}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+3,backrightoffset+3}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+4,backrightoffset+4}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+5,backrightoffset+5}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+6,backrightoffset+6}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+7,backrightoffset+7}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+8,backrightoffset+8}; nextline = nextline+1;
Line (nextline) = {frontrightoffset+9,backrightoffset+9}; nextline = nextline+1;
//Horizontal areas, counterclockwise from above
Line Loop (nextarea) = {rightoffset,firstsideline+1,-firstbackline,-firstsideline}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+2),firstsideline+3,firstbackline+2,-(firstsideline+2)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+5),firstsideline+5,firstbackline+5,-(firstsideline+4)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+8),firstsideline+7,firstbackline+8,-(firstsideline+6)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+11),firstsideline+9,firstbackline+11,-(firstsideline+8)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
// right side areas, outward pointing
Line Loop (nextarea) = {-(rightoffset+3),firstsideline+2,firstbackline+3,-firstsideline}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+6),firstsideline+4,firstbackline+6,-(firstsideline+2)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+9),firstsideline+6,firstbackline+9,-(firstsideline+4)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+12),firstsideline+8,firstbackline+12,-(firstsideline+6)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
//Right side areas, outward pointing
Line Loop (nextarea) = {-(rightoffset+1),firstsideline+1,firstbackline+1,-(firstsideline+3)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+4),firstsideline+3,firstbackline+4,-(firstsideline+5)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+7),firstsideline+5,firstbackline+7,-(firstsideline+7)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop (nextarea) = {-(rightoffset+10),firstsideline+7,firstbackline+10,-(firstsideline+9)}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;

//right line volumes, starting at the bottom, in order
// {front,bottom,right,back,top,right}
Surface Loop (nextvolume) = {rightoffset,-(rightoffset+8),rightoffset+17,-(rightoffset+4),rightoffset+9,rightoffset+13};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;
Surface Loop (nextvolume) = {rightoffset+1,-(rightoffset+9),rightoffset+18,-(rightoffset+5),rightoffset+10,rightoffset+14};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;
Surface Loop (nextvolume) = {rightoffset+2,-(rightoffset+10),rightoffset+19,-(rightoffset+6),rightoffset+11,rightoffset+15};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;
Surface Loop (nextvolume) = {rightoffset+3,-(rightoffset+11),rightoffset+20,-(rightoffset+7),rightoffset+12,rightoffset+16};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;


// Corners of the air space
firstdomainpoint = nextpointnum;
Point(nextpointnum) = {-domainXsize/2, 0, domainZsize/2, bigmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {domainXsize/2, 0, domainZsize/2, bigmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {domainXsize/2,domainYsize,domainZsize/2,bigmeshsiz*domainYsize/domainXsize}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {-domainXsize/2,domainYsize,domainZsize/2,bigmeshsiz*domainYsize/domainXsize}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {-domainXsize/2, 0, -domainZsize/2, bigmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {domainXsize/2, 0, -domainZsize/2, bigmeshsiz}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {domainXsize/2,domainYsize,-domainZsize/2,bigmeshsiz*domainYsize/domainXsize}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {-domainXsize/2,domainYsize,-domainZsize/2,bigmeshsiz*domainYsize/domainXsize}; nextpointnum = nextpointnum+1;
// Lines framing the air space
firstAirLine = nextline;
//front ones
Line (nextline) = {firstdomainpoint,firstdomainpoint+1}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+1,firstdomainpoint+2}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+2,firstdomainpoint+3}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+3,firstdomainpoint}; nextline = nextline+1;
//back ones
Line (nextline) = {firstdomainpoint+4,firstdomainpoint+5}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+5,firstdomainpoint+6}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+6,firstdomainpoint+7}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+7,firstdomainpoint+4}; nextline = nextline+1;
//front to back ones
Line (nextline) = {firstdomainpoint,firstdomainpoint+4}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+1,firstdomainpoint+5}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+3,firstdomainpoint+7}; nextline = nextline+1;
Line (nextline) = {firstdomainpoint+2,firstdomainpoint+6}; nextline = nextline+1;

//The easy areas framing the air space
firstAirArea = nextarea;
Printf("1st Air Area # = %f",firstAirArea);
//front plane
Line Loop(nextarea) = {firstAirLine,firstAirLine+1,firstAirLine+2,firstAirLine+3}; 
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
//back plane
Line Loop(nextarea) = {firstAirLine+4,firstAirLine+5,firstAirLine+6,firstAirLine+7};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
//left and right
Printf("Left plane loop: %f, %f, %f, %f",-(firstAirLine+3),firstAirLine+10,firstAirLine+7,-(firstAirLine+8));
Line Loop(nextarea) = {-(firstAirLine+3),firstAirLine+10,firstAirLine+7,-(firstAirLine+8)};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
Line Loop(nextarea) = {-(firstAirLine+1),firstAirLine+9,firstAirLine+5,-(firstAirLine+11)};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;
// top plane
Line Loop(nextarea) = {-(firstAirLine+2),firstAirLine+11,firstAirLine+6,-(firstAirLine+10)};
Plane Surface (nextarea) = {nextarea}; nextarea = nextarea+1;

//The bottom plane is more difficult because we must exclude holes for the lines.
Line Loop(nextarea) = {firstAirLine,firstAirLine+9,-(firstAirLine+4),-(firstAirLine+8)};
Plane Surface (nextarea) = {nextarea,leftoffset+8,rightoffset+8}; nextarea = nextarea+1;

// Define the airvolume
Printf("Air volume # = %f",nextvolume);
Surface Loop (nextvolume) = {
	-leftoffset,
	-(leftoffset+1),
	-(leftoffset+2),
	-(leftoffset+3),
	(leftoffset+4),
	(leftoffset+5),
	(leftoffset+6),
	(leftoffset+7),
	-(leftoffset+12),
	-(leftoffset+13),
	-(leftoffset+14),
	-(leftoffset+15),
	-(leftoffset+16),
	-(leftoffset+17),
	-(leftoffset+18),
	-(leftoffset+19),
	-(leftoffset+20),
		-rightoffset,
	-(rightoffset+1),
	-(rightoffset+2),
	-(rightoffset+3),
	(rightoffset+4),
	(rightoffset+5),
	(rightoffset+6),
	(rightoffset+7),
	-(rightoffset+12),
	-(rightoffset+13),
	-(rightoffset+14),
	-(rightoffset+15),
	-(rightoffset+16),
	-(rightoffset+17),
	-(rightoffset+18),
	-(rightoffset+19),
	-(rightoffset+20),
	firstAirArea,
	-(firstAirArea+1),
	firstAirArea+2,
	firstAirArea+3,
	firstAirArea+4,
	-(firstAirArea+5)
};
Volume (nextvolume) = {nextvolume}; nextvolume = nextvolume+1;

/* DEFINE CHARACTERISTIC LENGTHS

My strategy will be to use a MathEval field to make characteristic lengths that increase with 
distance from my region of interest asymptotically as r^3/2. Try this function:

fmin+Sqrt(1./r0)^3*((x-x0)^4/(sx^2+(x-x0)^2)+(y-y0)^4/(sy^2+(y-y0)^2)+(z-z0)^4/(sz^2+(z-z0)^2))^(0.75)

where x0,y0,z0 is the center of the region of interest, sx,sy,sz are the dimensions of this region, 
and r0 = rmax/maxval^(2/3) with rmax the distance to a distant point (e.g., domain edge) and maxval 
the desired characteristic length there.

Below I've substituted
r0 = 50 nm
fmin = 1 nm
[x0,y0,z0] = [37.,95.,0.]
[sx,sy,sz] = [0.5,0.5,0.5]
*/

/*
Field[1] = MathEval;
Field[1].F = "(1.e-9+Sqrt(20.)^3*((x-37.e-9)^4/(45.e-9+Fabs(x-37.e-9))^2+(y-95.e-9)^4/(30.e-9+Fabs(y-95.e-9))^2+z^4/(20.e-9+Fabs(z))^2)^0.75)";
*/


Field[1] = MathEval;
Field[1].F = "1.*(1.e-9+Sqrt(5.)*80.*((x-37.e-9)^4/(0.5e-9^2+(x-37.e-9)^2)+(y-95.e-9)^4/(0.5e-9^2+(y-95.e-9)^2)+z^4/(0.5e-9^2+z^2))^0.75)";


Mesh.CharacteristicLengthFromPoints = 0;
//Mesh.CharacteristicLengthFromCurvature = 1;
Background Field = 1;


/* Some extra, nongeometry points, to serve as markers for values checks */

firstReferencePoint = 200;
nextpointnum = firstReferencePoint;
Point(nextpointnum) = {0,20*nm,0}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,40*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,60*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,80*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,100*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,120*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,140*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,160*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,180*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {0,200*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,20*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,40*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,60*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,80*nm,0}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,100*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,120*nm,0}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,140*nm,0}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,160*nm,0}; nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,180*nm,0};  nextpointnum = nextpointnum+1;
Point(nextpointnum) = {20*nm,200*nm,0}; nextpointnum = nextpointnum+1;


/* Back to geometry */

Physical Volume(POSITIVELAYER) = {leftoffset+3,rightoffset+3};
Physical Volume(NEGATIVELAYER) = {leftoffset+1,rightoffset+1};
Physical Volume(LINES) = {leftoffset,leftoffset+2,rightoffset,rightoffset+2};
Physical Volume(VACUUM) = {rightoffset+4};

Physical Surface(GROUNDSURF) = {firstAirArea+5};
Physical Surface(INFSURF) = {firstAirArea+4};

