// Gmsh project created on Thu Feb 28 13:26:24 2019

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////// PREAMBLE /////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// using the OpenCASCADE kernel means we can't do boundary layer extrusions
// but it lets us do geometric boolean operations, without which creating 
// this mesh would be a nightmare
SetFactory("OpenCASCADE");

// perforated liner specification
linerRo = (8.625/2)*(0.0254);		// [m], liner outer radius
linerRi = linerRo - 0.009375;		// [m], liner inner radius
perforationDiameter = 0.375*0.0254;	// [m], diameter of the perforations in the liner
perforationsPerCirc = 16; 			// [-], number of perforations around a circumference (12)
perforationSpacing = 1*0.0254; 		// [m], space from center to center between rows of perforations (0.040 m)

// wire wrap specifications
wireGap = (0.14)*0.0254;			// [m], space between wires and liner
wireGauge = (0.09)*0.0254;			// [m], width of an individual wire (0.090")
wireHeight = (0.105)*0.0254;		// [m], height of an individual wire (0.105")
wireAperture = (0.03)*0.0254;		// [m], space between wires (0.008")

// mesh properties
expRate = 1.2; 							// expansion rate of node spacing away from sensitive regions
dxAperture = wireAperture/2;			// node spacing at the wire wrap apertures
dxGap =	5*dxAperture;					// node spacing at the top of the wire wrap gap layer
dxPerforation = perforationDiameter/20;	// node spacing at the perforations
dxExteriorRefine = dxGap;				// node spacing at the outer layer of the exterior refinement region
dxInteriorRefine = dxGap;				// node spacing at the inner layer of the interior refinement region





Mesh.Algorithm = 4; 				// frontal algorithm for surface meshes
Mesh.Algorithm3D = 1; 				// Delaunay algorithm for volume meshes
Geometry.LineNumbers=0;

Mesh.OptimizeThreshold = 0.9;

Mesh.MshFileVersion = 2.2;

interiorRefineR = linerRi - 2*wireGap;
exteriorRefineR = linerRo + 3*wireGap + wireHeight;

If (perforationSpacing>(linerRo*2*Pi/perforationsPerCirc))
	interiorFarFieldR = linerRi - perforationSpacing/2;
	exteriorFarFieldR = linerRo + perforationSpacing;
Else
	interiorFarFieldR = linerRi - (linerRo*2*Pi/perforationsPerCirc)/2;
	exteriorFarFieldR = linerRo + linerRo*2*Pi/perforationsPerCirc;
EndIf

// the following node spacings are calculated from prescribed parameters:

// node spacing at the exterior far field boundary
dxExteriorFar = 
	dxExteriorRefine 
	* 
	expRate^(
		Log(
			(
				(exteriorFarFieldR-exteriorRefineR)
				/
				dxExteriorRefine
			)
			*
			(expRate-1)
			+ 1
		)
		/
		Log(expRate)
		- 1
	);

// node spacing at the interior far field boundary
dxInteriorFar = 
	dxInteriorRefine 
	* 
	expRate^(
		Log(
			(
				(interiorRefineR-interiorFarFieldR)
				/
				dxExteriorRefine
			)
			*
			(expRate-1)
			+ 1
		)
		/
		Log(expRate)
		- 1
	);

dx = dxAperture;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////// LINER SURFACES AND HOLE VOLUMES /////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// radius of inner and outer surface
R = {};
R[0] = linerRi;
R[1] = linerRo;

// create the surface for the inner and outer side of the liner
For i In {0:1}
	linerC~{i} = newl; Circle(linerC~{i}) = {0,0,0,R[i],0,(2*Pi/perforationsPerCirc)};
	Rotate{{0,1,0},{0,0,0},-Pi/2}{ Curve{linerC~{i}}; }
	linerE~{i} = Extrude{perforationSpacing,0,0}{ Curve{linerC~{i}}; };
	linerS~{i} = linerE~{i}[1];
EndFor

// angle and position of each hole
X = {};
X[0] = 0;
X[1] = perforationSpacing;

Th = {};
Th[0] = 0;
Th[1] = 2*Pi/perforationsPerCirc;

For i In {0:1}
	holeC~{i} = newl; Circle(holeC~{i}) = {X[i],0,0,perforationDiameter/2};
	Rotate{{1,0,0},{0,0,0},-Th[i]}{ Curve{holeC~{i}}; }
	holeE~{i} = Extrude{0,2*linerRo*Sin(Th[i]),2*linerRo*Cos(Th[i])}{ Curve{holeC~{i}}; };
	holeS~{i} = holeE~{i}[1];
EndFor

For j In {0:1}
	For i In {0:1}
		// find the boolean difference between the surface and the hole
		BD = BooleanDifference{ Surface{linerS~{i}}; }{ Surface{holeS~{j}}; };
		// delete the original surface
		Delete{ Surface{linerS~{i}}; }
		// assign the tag for the resulting cut surface
		linerS~{i} = BD[1];
		// assign the tag for the hole in the surface
		holeS~{i}~{j} = BD[2];
		// extract the boundaries of the hole in the surface
		holeB~{i}~{j} = Boundary{ Surface{holeS~{i}~{j}}; };
	EndFor
	// delete the hole cylinder and everything related to it, not needed anymore
	Recursive Delete{ Surface{holeS~{j}}; }
EndFor

// get the labels of the lines bounding the resulting liner surfaces
For i In {0:1}
	linerB~{i} = Boundary{ Surface{linerS~{i}}; };
EndFor

// we want to move the entries in the hole boundary arrays around
// so that the 0 entry is always the arc in the center, the 1 entry
// is always the line alone the axial symmetry plane, and the 2 entry
// is the line along the angular symmetry plane
For i In {0:1}
	// the second hole is already correct, we just need to modify the
	// first hole's boundary array
	j = 0;
	tmp = holeB~{i}~{j}[];
	holeB~{i}~{j}[0] = tmp[2];
	holeB~{i}~{j}[1] = tmp[0];
	holeB~{i}~{j}[2] = tmp[1];
EndFor

// find the labels of the points that define the vertices of the holes
For j In {0:1}
	For i In {0:1}
		For k In {0:2}
			tmp = Abs( Boundary{ Curve{holeB~{i}~{j}[k]}; } );
			holeP~{i}~{j}[k] = tmp[0];
		EndFor
	EndFor
EndFor

// construct lines connecting the point pairs
For j In {0:1}
	For k In {0:2}
		holeL~{j}[k] = newl;
		Line(holeL~{j}[k]) = {holeP~{0}~{j}[k],holeP~{1}~{j}[k]};
	EndFor
EndFor

// construct the side surfaces of the holes
For j In {0:1}
	For k In {0:2}
		tmp = newll;
		Line Loop(tmp) = 
			{
				holeL~{j}[k%3],
				holeB~{0}~{j}[k%3],
				-holeL~{j}[(k+1)%3],
				-holeB~{1}~{j}[k%3]
			};
		holeSideS~{j}[k] = news;
		Surface(holeSideS~{j}[k]) = tmp;
	EndFor
EndFor

// construct the volumes for the holes
For j In {0:1}
	surfList = {};
	For k In {0:2}
		surfList[] += holeSideS~{j}[k];
	EndFor
	For i In {0:1}
		surfList[] += holeS~{i}~{j};
	EndFor
	tmp = newsl; Surface Loop(tmp) = surfList[];
	holeV[j] = newreg; Volume(holeV[j]) = tmp;
EndFor

// then delete the side surfaces because for some reason they stick around,
// and we don't want to accidentally use them for something instead of 
// correctly using Boundary{} on a volume and using the true surface bounds
For j In {0:1}
	For k In {0:2}
		Delete{ Surface{ holeSideS~{j}[k] }; }
	EndFor
EndFor

Coherence;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////// WIRE WRAP MESH LAYER //////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// The first step is to make a bunch of triangles on the axial-radial plane, one for each wire.

// It's important to calculate how many wires will fit in the domain first of all:
nWires = Floor((perforationSpacing)/(wireGauge+wireAperture));

// first let's just place all of the wire points
For i In {0:(nWires)}
	
	wireTopLP[i] = newp;
	Point(wireTopLP[i]) = {i*(wireAperture+wireGauge)-wireGauge/2,0,R[1]+wireGap+wireHeight,dx};
	
	wireTopRP[i] = newp;
	Point(wireTopRP[i]) = {(i+1)*(wireAperture+wireGauge)-wireGauge/2-wireAperture,0,R[1]+wireGap+wireHeight,dx};
	
	wireBotP[i] = newp;
	Point(wireBotP[i]) = {i*(wireAperture+wireGauge),0,R[1]+wireGap,dx};
	
EndFor

// now they need to be adjusted to make sense

// first, adjust the first top node so that it lines up with the x=0 plane
Translate{wireGauge/2,0,0}{Point{wireTopLP[0]};}

// next, do something similar to the last top node
Translate{-wireGauge/2,0,0}{Point{wireTopRP[#wireTopRP[]-1]};}

// now we're going to shift over half of the points so that the last wire
// lines up with the other bounding y-z plane
For i In {(Floor(nWires/2)+1):nWires}
	Translate{perforationSpacing-(nWires)*(wireAperture+wireGauge),0,0}{
		Point{wireTopRP[i],wireTopLP[i],wireBotP[i]};
	}
EndFor

// finally, we shift over the wire(s) in the middle half (or a third) as much,
// depending on whether there are an odd or even number of wires
If ((nWires%2)==0)
	Translate{(perforationSpacing-(nWires)*(wireAperture+wireGauge))/2,0,0}{
		Point{wireTopRP[Floor(nWires/2)],wireTopLP[Floor(nWires/2)],wireBotP[Floor(nWires/2)]};
	}
Else
	Translate{(perforationSpacing-(nWires)*(wireAperture+wireGauge))/3,0,0}{
		Point{wireTopRP[Floor(nWires/2)],wireTopLP[Floor(nWires/2)],wireBotP[Floor(nWires/2)]};
	}
	Translate{-(perforationSpacing-(nWires)*(wireAperture+wireGauge))/3,0,0}{
		Point{wireTopRP[Floor(nWires/2)+1],wireTopLP[Floor(nWires/2)+1],wireBotP[Floor(nWires/2)+1]};
	}
EndIf


// now that the triangles have been created, it's time to start connecting them with lines
For i In {0:nWires}
	
	wireTopL[2*i] = newl;
	Line(wireTopL[2*i]) = {wireTopLP[i],wireTopRP[i]};
	
	If (i<nWires)
		
		wireTopL[2*i+1] = newl;
		Line(wireTopL[2*i+1]) = {wireTopRP[i],wireTopLP[i+1]};
		
		wireBotL[i] = newl;
		Line(wireBotL[i]) = {wireBotP[i],wireBotP[i+1]};
		
	EndIf
	
	wireSideL[2*i] = newl;
	Line(wireSideL[2*i]) = {wireBotP[i],wireTopLP[i]};
	
	wireSideL[2*i+1] = newl;
	Line(wireSideL[2*i+1]) = {wireBotP[i],wireTopRP[i]};
	
EndFor

// create surfaces for the meshes between the wires
For i In {0:(#wireBotL[]-1)}
	tmp = newll;
	Line Loop(tmp) = {wireBotL[i],wireSideL[2*(i+1)],-wireTopL[2*i+1],-wireSideL[2*i+1]};
	wireRootS[i] = news;
	Surface(wireRootS[i]) = tmp;
EndFor

// extrude all of these surfaces by rotation to create the mesh between the wires
For i In {0:(#wireRootS[]-1)}
	wireE~{i} = Extrude{{1,0,0},{0,0,0},-2*Pi/perforationsPerCirc}{Surface{wireRootS[i]};};
EndFor

// also rotate the top lines of wires to eventually complete the surfaces at the top of the wires
For i In {0:nWires}
	tmp[] = Rotate{{1,0,0},{0,0,0},-2*Pi/perforationsPerCirc}
		{
			Duplicata{Curve{wireTopL[2*i]};}
		};
	wireTopRotL[i] = tmp[0];
EndFor

// we need to get label values for all of the angular top wire lines
For i In {1:(nWires-1)}
	// to find the left/right lines at the top of a wire, we need to:
	// 1. find the top boundary surface of the mesh region to the left/right of the wire
	// 2. find the right/left boundary line of said boundary surface
	tmpL = Abs( Boundary{ Surface{wireE~{i-1}[3]}; } );
	tmpR = Abs( Boundary{ Surface{wireE~{i}[3]}; } );
	// wireE~{<i-1/i>}[3] is the top surface of the left/right top surface
	// this stores the boundaries of that surface in 'tmp'
	// tmp<L/R>[0] is then the line on the 'left' side of the left/right top surface,
	// tmp<L/R>[2] is then the line on the 'right' side of the left/right top surface.
	wireTopAngLL[i] = tmpL[2];
	wireTopAngRL[i] = tmpR[0];
EndFor
// the first and last wires are special cases
// for the first wire, there is only a mesh region to the right:
i = 0;
tmpR = Abs( Boundary{ Surface{wireE~{i}[3]}; } );
wireTopAngRL[i] = tmpR[0];
// so then we create the left line by duplicating the right line:
tmp[] = Translate{-wireGauge/2,0,0}{Duplicata{Curve{tmpR[0]};}};
wireTopAngLL[i] = tmp[0];
// for the last wire, there is only a mesh to the left:
i = nWires;
tmpL = Abs( Boundary{ Surface{wireE~{i-1}[3]}; } );
wireTopAngLL[i] = tmpL[2];
// so then we create the right line by duplicating the left line:
tmp[] = Translate{wireGauge/2,0,0}{Duplicata{Curve{tmpL[2]};}};
wireTopAngRL[i] = tmp[0];

// we need to fix the node spacing on all of these lines


// now we are ready to create the surfaces on top of the wires
For i In {0:nWires}
//i=1;
	tmp = newll;
	Line Loop(tmp) = {-wireTopRotL[i],-wireTopAngLL[i],wireTopL[2*i],wireTopAngRL[i]};
	wireTopS[i] = news;
	Surface(wireTopS[i]) = tmp;
EndFor

Coherence;
// WARNING!
// using the Coherence command is necessary so that all of the curves interfacing between
// adjacent surfaces are unique and the refinement can be consistently controlled.
// HOWEVER, this means that from this point onwards, any array that held labels for curves
// or lines can't be trusted to work. The only thing that will remain useful are variables
// that contained labels for surfaces or volumes. If you need a curve somewhere after here,
// you will have to use the Boundary command on the relevant Surface.



//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// WIRE WRAP GAP MESH LAYER /////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// we need to make the surfaces on the axial and angular symmetry planes

// first step is to make some radial lines connecting the corners of the gap region
// in order to make these lines we need corner points

// we could either make new points since we know where they are in 3-space, or we can
// use the boundary command a few times to get the currently existing points.
// I choose the latter, because it is more robust if things above change

// the corner points will be stored in an array wireGapP~{i}~{j}~{k}, where
// i: 0 - bottom of layer, 1 - top of layer
// j: 0 - neg. axial sym plane, 1 - pos. axial sym plane
// k: 0 - neg. ang. sym plane, 1 - pos. ang. sym plane

//Hide{:}

// identify the points at the top of the liner mesh layer, wireGapP~{0}~{j}~{k}
For j In {0:1}
	tmp = Boundary{  Volume{ 	holeV[j] 	}; };
	tmp = Boundary{ Surface{ 	tmp[4] 		}; };
	tmp = Boundary{   Curve{ 	tmp[1] 		}; };
	wireGapP~{0}~{j}~{j} = tmp[j];
EndFor // verified correct

For j In {0:1}
	tmp = Boundary{ Surface{linerS~{1}}; };
	tmp = Boundary{ Curve{ tmp[3*j] }; };
	wireGapP~{0}~{j}~{1-j} = tmp[0];
EndFor // verified correct

// identify the corner points at the bottom of the wire mesh layer, wireGapP~{1}~{j}~{k}
For j In {0:1}
	tmp = Boundary{  Volume{ 	wireE~{(nWires-1)*j}[1] }; };
	tmp = Boundary{ Surface{ 	tmp[3] 					}; };
	tmp = Boundary{   Curve{ 	tmp[2*j] 				}; };
	For i In {0:1}
		wireGapP~{1}~{j}~{i} = tmp[(i+j)%2];
		// (i+j)%2 is basically XOR(i,j)
	EndFor
EndFor // verified correct

// next, create a line for each corner
For j In {0:1}
	For k In {0:1}
		wireGapL~{j}~{k} = newl;
		Line(wireGapL~{j}~{k}) = {wireGapP~{0}~{j}~{k},wireGapP~{1}~{j}~{k}};
	EndFor
EndFor

// now create the angular and axial plane surfaces

// axial planes
For j In {0:1}
	curveList = {};

	// the front corner curve
	curveList[] += wireGapL~{j}~{0};

	// the top curve
	tmp = Boundary{  Volume{ 	wireE~{(nWires-1)*j}[1] }; };
	tmp = Boundary{ Surface{ 	tmp[3] 					}; };
	curveList[] += Abs(tmp[2*j]);

	// the back corner curve
	curveList[] += -wireGapL~{j}~{1};

	// the bottom curve along the hole
	tmp = Boundary{  Volume{ 	holeV[j] 	}; };
	tmp = Boundary{ Surface{ 	tmp[4] 		}; };
	curveList[] += tmp[j];

	// the bottom curve along the liner
	tmp = Boundary{ Surface{linerS~{1}}; };
	curveList[] += tmp[3*j];
	
	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	wireGapAxialS~{j} = news;
	Surface(wireGapAxialS~{j}) = tmp;
	
EndFor

// angular planes
For j In {0:1}
	curveList = {};
	
	// the front corner curve
	curveList[] += wireGapL~{0}~{j};
	
	// all of the top curves from the bottom of the inter-wire meshes
	For k In {0:(nWires-1)}
		tmp = Boundary{  Volume{ 	wireE~{k}[1] 	}; };
		tmp = Boundary{ Surface{	tmp[3] 			}; };
		curveList[] += tmp[3-(2*j)];
	EndFor
	
	// the back corner curve
	curveList[] += -wireGapL~{1}~{j};
	
	// the bottom curve along the liner
	tmp = Boundary{ Surface{linerS~{1}}; };
	curveList[] += tmp[3*j+2];
	
	// the bottom curve along the hole
	tmp = Boundary{  Volume{ 	holeV[j] 	}; };
	tmp = Boundary{ Surface{ 	tmp[4] 		}; };
	curveList[] += tmp[j+1];
	
	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	wireGapAngularS~{j} = news;
	Surface(wireGapAngularS~{j}) = tmp;
	
EndFor

// now we combine all of this with existing surfaces to create the gap space mesh layer
surfList = {};

// liner surface
surfList[] += linerS~{1};

For j In {0:1}
	// hole surfaces
	tmp = Boundary{Volume{holeV[j]};};
	surfList[] += tmp[4];
	// angular and axial side surfaces
	surfList[] += wireGapAngularS~{j};
	surfList[] += wireGapAxialS~{j};
EndFor

// inter-wire mesh layer bottom surfaces
For i In {0:(nWires-1)}
	tmp = Boundary{Volume{wireE~{i}[1]};};
	surfList[] += tmp[3];
EndFor

// and make the volume
tmp = newsl;
Surface Loop(tmp) = surfList[];
wireGapV = newreg;
Volume(wireGapV) = tmp;

Coherence;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// INTERIOR REFINED REGION //////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// in this section, we will create a volume similar in size to the wire mesh gap region, but on the interior of
// the liner. this region will remain at the smallest refinement level, to resolve the flow coming out of the
// perforations and in the vicinity of the interior liner wall

// create a surface to bound this layer
interiorRefinementC = newl; Circle(interiorRefinementC) = {0,0,0,interiorRefineR,0,(2*Pi/perforationsPerCirc)};
Rotate{{0,1,0},{0,0,0},-Pi/2}{ Curve{interiorRefinementC}; }
interiorRefinementE = Extrude{perforationSpacing,0,0}{ Curve{interiorRefinementC}; };
interiorRefinementS = interiorRefinementE[1];

// identify the points at the bottom of the liner mesh layer, interiorRefinementP~{1}~{j}~{k}
For j In {0:1}
	tmp = Boundary{  Volume{ 	holeV[j] 	}; };
	tmp = Boundary{ Surface{ 	tmp[3] 		}; };
	tmp = Boundary{   Curve{ 	tmp[1] 		}; };
	interiorRefinementP~{1}~{j}~{j} = tmp[j];
EndFor // verified correct

For j In {0:1}
	tmp = Boundary{ Surface{linerS~{0}}; };
	tmp = Boundary{ Curve{ tmp[3*j] }; };
	interiorRefinementP~{1}~{j}~{1-j} = tmp[0];
EndFor // verified correct

// identify the points on the surface we just created, interiorRefinementP~{0}~{j}~{k}
tmp = Boundary{ Surface{interiorRefinementS}; };
For j In {0:1}
	For k In {0:1}
		tmpL = Abs(Boundary{ Curve{tmp[(3)-(3*k)-(j)+(2*j*k)]}; });
		// the function (3)-(3*k)-(j)+(2*j*k) was arrived at by observing the behaviour of
		// the output from the Boundary function applied to the interior refinement surface
		// and matching up the indices of the boundary curves to the desired j,k indices
		interiorRefinementP~{0}~{j}~{k} = tmpL[0];
	EndFor
EndFor

// next, create a line for each corner
For j In {0:1}
	For k In {0:1}
		interiorRefinementL~{j}~{k} = newl;
		Line(interiorRefinementL~{j}~{k}) = {interiorRefinementP~{0}~{j}~{k},interiorRefinementP~{1}~{j}~{k}};
	EndFor
EndFor

// now create the angular and axial plane surfaces

// axial planes
For j In {0:1}
	curveList = {};

	// the front corner curve
	curveList[] += interiorRefinementL~{j}~{0};
	
	// the purpose of the if statements is to make sure that the curves
	// are added to the curve list in the right order
	If (j==0)
		// the top curve along the hole
		tmp = Boundary{Volume{holeV[j]};};
		tmp = Abs(Boundary{Surface{tmp[3]};});
		curveList[] += tmp[j];
	EndIf

	// the top curve along the liner
	tmp = Abs(Boundary{ Surface{linerS~{0}};});
	curveList[] += -tmp[3*j];
	
	If (j==1)
		// the top curve along the hole
		tmp = Boundary{Volume{holeV[j]};};
		tmp = Abs(Boundary{Surface{tmp[3]};});
		curveList[] += tmp[j];
	EndIf
	
	// the back corner curve
	curveList[] += -interiorRefinementL~{j}~{1};

	// the bottom curve along the refinement surface
	tmp = Abs(Boundary{Surface{interiorRefinementS};});
	curveList[] += tmp[(2*j+3)%4];
	
	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	interiorRefinementAxialS~{j} = news;
	Surface(interiorRefinementAxialS~{j}) = tmp;
EndFor


// angular planes

For j In {0:1}
	curveList = {};

	// the front corner curve
	curveList[] += interiorRefinementL~{0}~{j};

	// the purpose of the if statements is to make sure that the curves
	// are added to the curve list in the right order
	If (j==0)
		// the top curve along the hole
		tmp = Boundary{Volume{holeV[j]};};
		tmp = Abs(Boundary{Surface{tmp[3]};});
		curveList[] += tmp[1];
	EndIf

	// the top curve along the liner
	tmp = Abs(Boundary{ Surface{linerS~{0}};});
	curveList[] += -tmp[3*j+2];

	If (j==1)
		// the top curve along the hole
		tmp = Boundary{Volume{holeV[j]};};
		tmp = Abs(Boundary{Surface{tmp[3]};});
		curveList[] += tmp[2];
	EndIf

	// the back corner curve
	curveList[] += -interiorRefinementL~{1}~{j};

	// the bottom curve along the refinement surface
	tmp = Abs(Boundary{Surface{interiorRefinementS};});
	curveList[] += tmp[(2*j+2)%4];

	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	interiorRefinementAngularS~{j} = news;
	Surface(interiorRefinementAngularS~{j}) = tmp;
EndFor

// now we combine all of this with existing surfaces to create the internal refined mesh layer
surfList = {};

// liner surface
surfList[] += linerS~{0};

For j In {0:1}
	// hole surfaces
	tmp = Boundary{Volume{holeV[j]};};
	surfList[] += tmp[3];
	// angular and axial side surfaces
	surfList[] += interiorRefinementAngularS~{j};
	surfList[] += interiorRefinementAxialS~{j};
EndFor

// bottom refinement surface
surfList[] += interiorRefinementS;

// and make the volume
tmp = newsl;
Surface Loop(tmp) = surfList[];
interiorRefinementV = newreg;
Volume(interiorRefinementV) = tmp;

Coherence;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// EXTERIOR REFINED REGION //////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// similarly to above, we are going to create a small refined zone but this time it will be exterior to the wire
// mesh layer to try and resolve the flow converging on the wires

// create a surface to bound this layer
exteriorRefinementC = newl; Circle(exteriorRefinementC) = {0,0,0,exteriorRefineR,0,(2*Pi/perforationsPerCirc)};
Rotate{{0,1,0},{0,0,0},-Pi/2}{ Curve{exteriorRefinementC}; }
exteriorRefinementE = Extrude{perforationSpacing,0,0}{ Curve{exteriorRefinementC}; };
exteriorRefinementS = exteriorRefinementE[1];


// identify the points on the surface we just created, exteriorRefinementP~{1}~{j}~{k}
tmp = Boundary{ Surface{exteriorRefinementS}; };
For j In {0:1}
	For k In {0:1}
		tmpL = Abs(Boundary{ Curve{tmp[(3)-(3*k)-(j)+(2*j*k)]}; });
		// the function (3)-(3*k)-(j)+(2*j*k) was arrived at by observing the behaviour of
		// the output from the Boundary function applied to the exterior refinement surface
		// and matching up the indices of the boundary curves to the desired j,k indices
		exteriorRefinementP~{1}~{j}~{k} = tmpL[0];
	EndFor
EndFor

// identify the corner points on the outside of the wire layer
For j In {0:1}
	tmp = Boundary{Surface{wireTopS[j*nWires]};};
	tmp = Boundary{Curve{tmp[2*(1-j)]};};
	For k In {0:1}
		exteriorRefinementP~{0}~{j}~{k} = tmp[1-((j+k)%2)];
		// 1-((j+k)%2) is basically NOT(XOR(j,k))
	EndFor
EndFor

// next, create a line for each corner
For j In {0:1}
	For k In {0:1}
		exteriorRefinementL~{j}~{k} = newl;
		Line(exteriorRefinementL~{j}~{k}) = {exteriorRefinementP~{0}~{j}~{k},exteriorRefinementP~{1}~{j}~{k}};
	EndFor
EndFor

// now create the angular and axial plane surfaces

// axial planes
For j In {0:1}
	curveList = {};
	
	// front corner curve
	curveList[] += exteriorRefinementL~{j}~{0};
	
	// top curve
	tmp = Boundary{Surface{exteriorRefinementS};};
	curveList[] += tmp[2*(1-j)+1];
	
	// the back corner curve
	curveList[] += -exteriorRefinementL~{j}~{1};
	
	// bottom curve
	tmp = Boundary{Surface{wireTopS[j*nWires]};};
	curveList[] += tmp[2*(1-j)];
	
	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	exteriorRefinementAxialS~{j} = news;
	Surface(exteriorRefinementAxialS~{j}) = tmp;
	
EndFor

// angular planes
For k In {0:1}
	curveList = {};
	
	// front corner curve
	curveList[] += exteriorRefinementL~{0}~{k};
	
	// top curve
	tmp = Boundary{Surface{exteriorRefinementS};};
	curveList[] += tmp[2*(1-k)];
	
	// back corner curve
	curveList[] += exteriorRefinementL~{1}~{k};
	
	// bottom curves
	For i In {0:(nWires-1)}
		tmp = Boundary{Surface{wireTopS[i]};}; 
		curveList[] += tmp[3-2*k];
		tmp = Boundary{Volume{wireE~{i}[1]};};
		tmp = Boundary{Surface{tmp[1]};};
		curveList[] += tmp[3-2*k];
	EndFor
	tmp = Boundary{Surface{wireTopS[nWires]};}; 
	curveList[] += tmp[3-2*k];
	
	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	exteriorRefinementAngularS~{k} = news;
	Plane Surface(exteriorRefinementAngularS~{k}) = tmp;
EndFor

// now we combine all of this with existing surfaces to create the external refined mesh layer
surfList = {};

// wire mesh surfaces
For i In {0:(nWires-1)}
	surfList[] += wireTopS[i];
	tmp = Boundary{Volume{wireE~{i}[1]};};
	surfList[] += tmp[1];
EndFor
surfList[] += wireTopS[nWires];

// axial and angular side surfaces
For i In {0:1}
	surfList[] += exteriorRefinementAngularS~{i};
	surfList[] += exteriorRefinementAxialS~{i};
EndFor

// top refinement surface
surfList[] += exteriorRefinementS;

// and make the volume
tmp = newsl;
Surface Loop(tmp) = surfList[];
exteriorRefinementV = newreg;
Volume(exteriorRefinementV) = tmp;

Delete{Surface{exteriorRefinementS};}

Coherence;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////// EXTERIOR FAR-FIELD REGION /////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// create a surface to bound this layer
exteriorFarFieldC = newl; Circle(exteriorFarFieldC) = {0,0,0,exteriorFarFieldR,0,(2*Pi/perforationsPerCirc)};
Rotate{{0,1,0},{0,0,0},-Pi/2}{ Curve{exteriorFarFieldC}; }
exteriorFarFieldE = Extrude{perforationSpacing,0,0}{ Curve{exteriorFarFieldC}; };
exteriorFarFieldS = exteriorFarFieldE[1];

// identify the points on the surface we just created, exteriorFarFieldP~{1}~{j}~{k}
tmp = Boundary{ Surface{exteriorFarFieldS}; };
For j In {0:1}
	For k In {0:1}
		tmpL = Abs(Boundary{ Curve{tmp[(3)-(3*k)-(j)+(2*j*k)]}; });
		// the function (3)-(3*k)-(j)+(2*j*k) was arrived at by observing the behaviour of
		// the output from the Boundary function applied to the exterior refinement surface
		// and matching up the indices of the boundary curves to the desired j,k indices
		exteriorFarFieldP~{1}~{j}~{k} = tmpL[0];
	EndFor
EndFor

// identify the points on the adjacent refinement surface
For j In {0:1}
	tmp = Boundary{Volume{exteriorRefinementV};};
	tmp = Boundary{Surface{tmp[2*nWires+5]};};
	tmp = Boundary{Curve{tmp[3-2*j]};};
	For k In {0:1}
		exteriorFarFieldP~{0}~{j}~{k} = tmp[(j+k)%2];
		// (j+k)%2 is basically XOR(j,k)
	EndFor
EndFor

// next, create a line for each corner
For j In {0:1}
	For k In {0:1}
		exteriorFarFieldL~{j}~{k} = newl;
		Line(exteriorFarFieldL~{j}~{k}) = {exteriorFarFieldP~{0}~{j}~{k},exteriorFarFieldP~{1}~{j}~{k}};
	EndFor
EndFor

// now create the angular and axial plane surfaces

// axial planes
For j In {0:1}
	curveList = {};
	
	// front corner curve
	curveList[] += exteriorFarFieldL~{j}~{0};
	
	// top curve
	tmp = Boundary{Surface{exteriorFarFieldS};};
	curveList[] += tmp[3-2*j];
	
	// back corner curve
	curveList[] += exteriorFarFieldL~{j}~{1};
	
	// bottom curve
	tmp = Boundary{Volume{exteriorRefinementV};};
	tmp = Boundary{Surface{tmp[2*nWires+5]};};
	curveList[] += tmp[3-2*j];
	
	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	exteriorFarFieldAxialS~{j} = news;
	Surface(exteriorFarFieldAxialS~{j}) = tmp;
EndFor

// angular planes
For k In {0:1}
	curveList = {};
	
	// front corner curve
	curveList[] += exteriorFarFieldL~{0}~{k};
	
	// top curve
	tmp = Boundary{Surface{exteriorFarFieldS};};
	curveList[] += tmp[2*(1-k)];
	
	// back corner curve
	curveList[] += exteriorFarFieldL~{1}~{k};
	
	// bottom curve
	tmp = Boundary{Volume{exteriorRefinementV};};
	tmp = Boundary{Surface{tmp[2*nWires+5]};};
	curveList[] += tmp[2*(1-k)];

	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	exteriorFarFieldAngularS~{k} = news;
	Plane Surface(exteriorFarFieldAngularS~{k}) = tmp;
EndFor

// now we combine all of this with existing surfaces to create the external far field mesh layer
surfList = {};

surfList[] += exteriorFarFieldS;

For i In {0:1}
	surfList[] += exteriorFarFieldAngularS~{i};
	surfList[] += exteriorFarFieldAxialS~{i};
EndFor

tmp = Boundary{Volume{exteriorRefinementV};};
surfList[] += tmp[2*nWires+5];

tmp = newsl;
Surface Loop(tmp) = surfList[];
exteriorFarFieldV = newreg;
Volume(exteriorFarFieldV) = tmp;

Delete{Surface{exteriorFarFieldS};}

Coherence;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////// INTERIOR FAR-FIELD REGION /////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// create a surface to bound this layer
interiorFarFieldC = newl; Circle(interiorFarFieldC) = {0,0,0,interiorFarFieldR,0,(2*Pi/perforationsPerCirc)};
Rotate{{0,1,0},{0,0,0},-Pi/2}{ Curve{interiorFarFieldC}; }
interiorFarFieldE = Extrude{perforationSpacing,0,0}{ Curve{interiorFarFieldC}; };
interiorFarFieldS = interiorFarFieldE[1];

// identify the points on the surface we just created, interiorFarFieldP~{0}~{j}~{k}
tmp = Boundary{ Surface{interiorFarFieldS}; };
For j In {0:1}
	For k In {0:1}
		tmpL = Abs(Boundary{ Curve{tmp[(3)-(3*k)-(j)+(2*j*k)]}; });
		// the function (3)-(3*k)-(j)+(2*j*k) was arrived at by observing the behaviour of
		// the output from the Boundary function applied to the exterior refinement surface
		// and matching up the indices of the boundary curves to the desired j,k indices
		interiorFarFieldP~{0}~{j}~{k} = tmpL[0];
	EndFor
EndFor

// identify the points on the adjacent refinement surface
For j In {0:1}
	tmp = Boundary{Volume{interiorRefinementV};};
	tmp = Boundary{Surface{tmp[7]};};
	tmp = Boundary{Curve{tmp[3-2*j]};};
	For k In {0:1}
		interiorFarFieldP~{1}~{j}~{k} = tmp[(j+k)%2];
		// (j+k)%2 is basically XOR(j,k)
	EndFor
EndFor

// next, create a line for each corner
For j In {0:1}
	For k In {0:1}
		interiorFarFieldL~{j}~{k} = newl;
		Line(interiorFarFieldL~{j}~{k}) = {interiorFarFieldP~{0}~{j}~{k},interiorFarFieldP~{1}~{j}~{k}};
	EndFor
EndFor

// now create the angular and axial plane surfaces

// axial planes
For j In {0:1}
	curveList = {};
	
	// front corner curve
	curveList[] += interiorFarFieldL~{j}~{0};
	
	// top curve
	tmp = Boundary{Volume{interiorRefinementV};};
	tmp = Boundary{Surface{tmp[7]};};
	curveList[] += tmp[3-2*j];
	
	// back corner curve
	curveList[] += interiorFarFieldL~{j}~{1};
	
	// bottom curve
	tmp = Boundary{Surface{interiorFarFieldS};};
	curveList[] += tmp[3-2*j];
	
	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	interiorFarFieldAxialS~{j} = news;
	Surface(interiorFarFieldAxialS~{j}) = tmp;
EndFor

// angular planes
For k In {0:1}
	curveList = {};
	
	// front corner curve
	curveList[] += interiorFarFieldL~{0}~{k};
	
	// top curve
	tmp = Boundary{Volume{interiorRefinementV};};
	tmp = Boundary{Surface{tmp[7]};};
	curveList[] += tmp[2*(1-k)];
	
	// back corner curve
	curveList[] += interiorFarFieldL~{1}~{k};
	
	// bottom curve
	tmp = Boundary{Surface{interiorFarFieldS};};
	curveList[] += tmp[2*(1-k)];
	
	// now make the surface
	tmp = newll;
	Line Loop(tmp) = curveList[];
	interiorFarFieldAngularS~{k} = news;
	Plane Surface(interiorFarFieldAngularS~{k}) = tmp;
EndFor

// now we combine all of this with existing surfaces to create the external far field mesh layer
surfList = {};

surfList[] += interiorFarFieldS;

For i In {0:1}
	surfList[] += interiorFarFieldAngularS~{i};
	surfList[] += interiorFarFieldAxialS~{i};
EndFor

tmp = Boundary{Volume{interiorRefinementV};};
surfList[] += tmp[7];

tmp = newsl;
Surface Loop(tmp) = surfList[];
interiorFarFieldV = newreg;
Volume(interiorFarFieldV) = tmp;

Delete{Surface{interiorFarFieldS};}

Coherence;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// NODE SPACING SPECIFICATION ///////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// set node spacing for all of our volumes by modifying the boundary point spacings
// the order of definition matters since subsequent commands will overwrite previous
// commands, so we will set the spacing of non-refined regions first and then overwrite
// the spacing on shared boundaries by specifying refined region spacing last


 // far fields first
Characteristic Length { PointsOf{ Volume{ interiorFarFieldV }; } } = dxInteriorFar;
Characteristic Length { PointsOf{ Volume{ exteriorFarFieldV }; } } = dxExteriorFar;

// then interior refinement regions
Characteristic Length { PointsOf{ Volume{ interiorRefinementV }; } } = dxInteriorRefine;
Characteristic Length { PointsOf{ Volume{ exteriorRefinementV }; } } = dxExteriorRefine;

// finally, space between wires
For i In {0:(nWires-1)}
	Characteristic Length { PointsOf{ Volume{ wireE~{i}[1] }; } } = dxAperture;
EndFor
// additionally, we need to correct the spacing on a couple of the surfaces
// bounding the exterior refinement region
tmp = Boundary{ Volume{ exteriorRefinementV }; };
Characteristic Length { PointsOf{ Surface{ tmp[0] 			}; } } = dxAperture;
Characteristic Length { PointsOf{ Surface{ tmp[2*nWires]	}; } } = dxAperture;

// then gap between liner and wire wrap
Characteristic Length { PointsOf{ Volume{ wireGapV }; } } = dxGap;

// then perforations
For i In {0:1}
	Characteristic Length { PointsOf{ Volume{ holeV[i] }; } } = dxPerforation;
EndFor

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// PHYSICAL VOLUME ASSIGNMENT ///////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// exterior far field
Physical Volume("outerFar",newreg) = exteriorFarFieldV;

// exterior refinement region
Physical Volume("outerRefine",newreg) = exteriorRefinementV;

// wire wrap space
volList = {};
For i In {0:(nWires-1)}
	volList[] += wireE~{i}[1];
EndFor
Physical Volume("wireWrap",newreg) = volList[];

// wire gap space
Physical Volume("wireGap",newreg) = wireGapV;

// perforations
volList = {};
For i In {0:1}
	volList[] += holeV[i];
EndFor
Physical Volume("perforations",newreg) = volList[];

// interior refinement region
Physical Volume("innerRefine",newreg) = interiorRefinementV;

// interior far field
Physical Volume("innerFar",newreg) = interiorFarFieldV;

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////// BOUNDARY SURFACE ASSIGNMENT //////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// exterior far field
tmp = Boundary{ Volume{ exteriorFarFieldV }; };
Physical Surface("exterior") = tmp[0];

// interior far field
tmp = Boundary{ Volume{ interiorFarFieldV }; };
Physical Surface("interior") = tmp[0];

// walls
surfList = {};
tmp = Boundary{ Volume{ interiorRefinementV }; };
surfList[] += tmp[0];
tmp = Boundary{ Volume{ wireGapV }; };
surfList[] += tmp[0];
For i In {0:1}
	tmp = Boundary{ Volume{ holeV[i] }; };
	surfList[] += tmp[0];
EndFor
For i In {0:(nWires-1)}
	tmp = Boundary{ Volume{ wireE~{i}[1] }; };
	surfList[] += tmp[0];
	surfList[] += tmp[2];
EndFor
tmp = Boundary{ Volume{ exteriorRefinementV }; };
For i In {0:nWires}
	surfList[] += tmp[2*i];
EndFor
Physical Surface("walls") = surfList[];

// negative axial symmetry plane
surfList = {};
tmp = Boundary{ Volume{ interiorFarFieldV }; };
surfList[] += tmp[2];
tmp = Boundary{ Volume{ interiorRefinementV }; };
surfList[] += tmp[3];
tmp = Boundary{ Volume{ holeV[0] }; };
surfList[] += tmp[1];
tmp = Boundary{ Volume{ wireGapV }; };
surfList[] += tmp[3];
tmp = Boundary{ Volume{ exteriorRefinementV }; };
surfList[] += tmp[2*nWires+2];
tmp = Boundary{ Volume{ exteriorFarFieldV }; };
surfList[] += tmp[2];
Physical Surface("axialSym1") = surfList[];

// positive axial symmetry plane
surfList = {};
tmp = Boundary{ Volume{ interiorFarFieldV }; };
surfList[] += tmp[4];
tmp = Boundary{ Volume{ interiorRefinementV }; };
surfList[] += tmp[6];
tmp = Boundary{ Volume{ holeV[1] }; };
surfList[] += tmp[1];
tmp = Boundary{ Volume{ wireGapV }; };
surfList[] += tmp[6];
tmp = Boundary{ Volume{ exteriorRefinementV }; };
surfList[] += tmp[2*nWires+4];
tmp = Boundary{ Volume{ exteriorFarFieldV }; };
surfList[] += tmp[4];
Physical Surface("axialSym2") = surfList[];

// negative angular symmetry plane
surfList = {};
tmp = Boundary{ Volume{ interiorFarFieldV }; };
surfList[] += tmp[1];
tmp = Boundary{ Volume{ interiorRefinementV }; };
surfList[] += tmp[2];
tmp = Boundary{ Volume{ holeV[0] }; };
surfList[] += tmp[2];
tmp = Boundary{ Volume{ wireGapV }; };
surfList[] += tmp[2];
For i In {0:(nWires-1)}
	tmp = Boundary{ Volume{ wireE~{i}[1] }; };
	surfList[] += tmp[4];
EndFor
tmp = Boundary{ Volume{ exteriorRefinementV }; };
surfList[] += tmp[2*nWires+1];
tmp = Boundary{ Volume{ exteriorFarFieldV }; };
surfList[] += tmp[1];
Physical Surface("angularSym1") = surfList[];


// positive angular symmetry plane
surfList = {};
tmp = Boundary{ Volume{ interiorFarFieldV }; };
surfList[] += tmp[3];
tmp = Boundary{ Volume{ interiorRefinementV }; };
surfList[] += tmp[5];
tmp = Boundary{ Volume{ holeV[1] }; };
surfList[] += tmp[2];
tmp = Boundary{ Volume{ wireGapV }; };
surfList[] += tmp[5];
For i In {0:(nWires-1)}
	tmp = Boundary{ Volume{ wireE~{i}[1] }; };
	surfList[] += tmp[5];
EndFor
tmp = Boundary{ Volume{ exteriorRefinementV }; };
surfList[] += tmp[2*nWires+3];
tmp = Boundary{ Volume{ exteriorFarFieldV }; };
surfList[] += tmp[3];
Physical Surface("angularSym2") = surfList[];





//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////// EXTRAS /////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/* // some code verifying that we still have variables defining our surfaces and volumes after the coherence commands
For j In {0:1}
	Printf('hole %0.0f is volume %0.0f',j,holeV[j]);
EndFor

For i In {0:(nWires-1)}
	Printf('wire mesh region %0.0f is volume %0.0f',i,wireE~{i}[1]);
EndFor

For i In {0:nWires}
	Printf('wire top surface mesh %0.0f is surface %0.0f',i,wireTopS[i]);
EndFor */



