Mesh.RecombineAll = 1;
Mesh.RecombinationAlgorithm = 0;

Li = 50.0; // distance of inflow boundary from origin
Lo = 50.0; // distance of outflow boundary from origin

n  = 50; // points on upper/lower surface of airfoil used to define airfoil
         // These points may not appear in the mesh.

lc1 = 10.0;
lc2 = 0.05;

m = 2*n - 2; // total number of points on airfoil without repetition
             // LE and TE points are common to upper/lower surface

nle = n; // point number of LE = no. of points on upper surface
         // Point(1) is trailing edge

// NACA0012 profile
// formula taken from http://turbmodels.larc.nasa.gov/naca0012_val.html
Macro NACA0012
   x2 = x * x;
   x3 = x * x2;
   x4 = x * x3;
   y = 0.594689181*(0.298222773*Sqrt(x) 
       - 0.127125232*x - 0.357907906*x2 + 0.291984971*x3 - 0.105174606*x4);
Return

// put points on upper surface of airfoil
For i In {1:n}
   theta = Pi * (i-1) / (n-1);
   x = 0.5 * (Cos(theta) + 1.0);
   Call NACA0012;
   Point(i) = {x, y, 0.0, lc2};
   xx[i] = x;
   yy[i] = y;
EndFor

// put points on lower surface of airfoil, use upper surface points and reflect
For i In {n+1:m}
   Point(i) = {xx[2*n-i], -yy[2*n-i], 0.0, lc2};
EndFor

Spline(1) = {1:n};
Spline(2) = {n:m, 1};

Point(1000) = {1.0, Li, 0.0,lc1};
Point(1001) = {0.0, Li, 0.0,lc1};
Point(1002) = {-Li, 0.0, 0.0,lc1};
Point(1003) = {0.0, -Li, 0.0,lc1};
Point(1004) = {1.0, -Li, 0.0,lc1};

Point(1005) = {Lo, 0.0, 0.0,lc1};
Point(1006) = {Lo, Li, 0.0,lc1};
Point(1007) = {Lo, -Li, 0.0,lc1};

Line(5) = {1000, 1001};
Circle(6) = {1001, nle, 1002};
Circle(7) = {1002, nle, 1003};
Line(8) = {1003, 1004};
Line(13) = {1005, 1006};
Line(14) = {1005, 1007};
Line(15) = {1006, 1000};
Line(16) = {1007, 1004};

Line Loop(1) = {1,2};
Line Loop(2) = {5,6,7,8,-16,-14,13,15};

Plane Surface(201) = {2,1};

//Define Boundary Layer
Field[1] = BoundaryLayer;
Field[1].EdgesList = {1,2};
Field[1].NodesList = {1,m};
Field[1].FanNodesList = {1,n};
Field[1].hfar = 0.05;
Field[1].hwall_n = 0.009;
Field[1].thickness = 0.2;
Field[1].ratio = 1.1;
Field[1].AnisoMax = 10;
Field[1].Quads = 1;
Field[1].IntersectMetrics = 0;
BoundaryLayer Field = 1;
