Mesh.Algorithm = 6;
//Mesh.Algorithm3D = 4;
Mesh.Optimize = 1;
Mesh.SmoothNormals = 1;
Mesh.CharacteristicLengthFromCurvature = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 1;
Mesh.CharacteristicLengthFactor = 1;

dx = 0.05;
dy = 0.01;
dfar = 8.0*dx;
BL_thickness = 0.1;
lip_thickness = 0.05;
NY=Ceil(BL_thickness/dy);
Printf("NY = %g", NY);
NY_lip=Ceil(lip_thickness/dy)+1;
Printf("NY_lip = %g", NY_lip);
NZ_4th = Ceil(Pi/(4.*dx));
Printf("NZ_4th = %g", NZ_4th);

Mesh.CharacteristicLengthMin = dx;
Mesh.CharacteristicLengthMax = dfar;

//nozzle inner wall points
p1 = newp;
Point(p1) = {0.0, 0.5, 0.0, dx};
Point(p1+1) = {-1.0, 0.5, 0.0, dx};
Point(p1+2) = {-2.0, 0.5, 0.0, dx};

// BL box point outside the nozzle lip surface
Point(p1+3) = {0.5, 0.5, 0.0, dx};

//nozzle outside wall points
Point(p1+4) = {-2.0, 0.6, 0.0, dx};
Point(p1+5) = {0.0, 0.6, 0.0, dx};

// farfield points
Point(p1+6) = {-2.0, 3.0, 0.0, dfar};
Point(p1+7) = {4.0, 3.0, 0.0, dfar};

line1 = newl;
Line(line1) = {p1,p1+1};
Line(line1+1) = {p1+3,p1};

// four layers inside the BL mesh
BL_surface[] = Extrude {0, -BL_thickness, 0} { Line{line1, line1+1}; Layers{NY};};

//BL layer surfaces and volume inside nozzle
k = 0;
surf_extrude_index = BL_surface[1];
For i In {0:3}
        Nozzle_BL_volumes[] = Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} { Surface{surf_extrude_index};Layers{NZ_4th};};
        nozzle_BL_volumes_index_array[i+k*4] = Nozzle_BL_volumes[1];
        nozzle_BL_wall_surfaces_index_array[i] = Nozzle_BL_volumes[2];
        nozzle_BL_edge_surfaces_index_array[i] = Nozzle_BL_volumes[4];
        back_BL_surfaces_index_array[i] = Nozzle_BL_volumes[3];
        surf_extrude_index = Nozzle_BL_volumes[0];
EndFor

// This gives me trouble
//BL layer surfaces and BL volume outside the nozzle
k=1;
surf_extrude_index = BL_surface[5];
For i In {0:3}
        Nozzle_BL_volumes[] = Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} { Surface{surf_extrude_index};Layers{NZ_4th};};
        nozzle_BL_volumes_index_array[i+k*4] = Nozzle_BL_volumes[1];
        outside_front_BL_surfaces_index_array[i] = Nozzle_BL_volumes[5];
        outside_BL_bottom_surfaces_index_array[i] = Nozzle_BL_volumes[4];
        outside_BL_top_surfaces_index_array[i] = Nozzle_BL_volumes[2];
        surf_extrude_index = Nozzle_BL_volumes[0];
EndFor

// end of trouble section

//inlet surface
line1 = newl;
Line(line1) = {p1+1,p1+2};
line_extrude_index = line1;
For i In {0:3}
    extrude_surfaces[] = Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} { Line{line_extrude_index};};
    entry_wall_surfaces_index_array[i] = extrude_surfaces[1];
    nozzle_inlet_lines_array[i] = extrude_surfaces[2];
    line_extrude_index = extrude_surfaces[0];
EndFor

line1=newl;
Line Loop(line1) = {nozzle_inlet_lines_array[]};
nozzle_inlet_surface = news;
Plane Surface(nozzle_inlet_surface) = {line1};

//nozzle lip surface
nozzle_lip = newl;
Line(nozzle_lip) = {p1, p1+5};
Transfinite Line {nozzle_lip} = NY_lip Using Progression 1.0;

line_extrude_index = nozzle_lip;
For i In {0:3}
    extrude_surfaces[] = Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} { Line{line_extrude_index};Layers{NZ_4th};};
    lip_surfaces_index_array[i] = extrude_surfaces[1];
    line_extrude_index = extrude_surfaces[0];
EndFor

//For j In {0:#lip_surfaces_index_array[]-1}
//Printf("lip surfaces = %g", lip_surfaces_index_array[j]);
//EndFor

//nozzle outside side wall
line1 = newl;
Line(line1) = {p1+4,p1+5};

line_extrude_index = line1;
For i In {0:3}
    extrude_surfaces[] = Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} { Line{line_extrude_index};};
    nozzle_outside_wall_surfaces_index_array[i] = extrude_surfaces[1];
    nozzle_outside_inlet_lines_array[i] = extrude_surfaces[3];
    line_extrude_index = extrude_surfaces[0];
EndFor

//farfield;
line1 = newl;
Line(line1) = {p1+6,p1+7};

line_extrude_index = line1;
For i In {0:3}
    extrude_surfaces[] = Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} { Line{line_extrude_index};};
    farfield_sideline_surfaces_index_array[i] = extrude_surfaces[1];
    line_extrude_index = extrude_surfaces[0];
    farfield_outflow_lines_array[i] = extrude_surfaces[2];
    farfield_inflow_lines_array[i] = extrude_surfaces[3];
EndFor

line1=newl;
Line Loop(line1) = {nozzle_outside_inlet_lines_array[], farfield_inflow_lines_array[]};
farfield_inflow_surface = news;
Plane Surface(farfield_inflow_surface) = {line1};

line1=newl;
Line Loop(line1) = {farfield_outflow_lines_array[]};
farfield_outflow_surface = news;
Plane Surface(farfield_outflow_surface) = {line1};

unstructured_Volume_surfaces_array[] = {nozzle_outside_wall_surfaces_index_array[], 
					lip_surfaces_index_array[], 
					outside_BL_top_surfaces_index_array[], outside_front_BL_surfaces_index_array[],
					outside_BL_bottom_surfaces_index_array[], nozzle_BL_edge_surfaces_index_array[],
				 	back_BL_surfaces_index_array[], entry_wall_surfaces_index_array[], 
					nozzle_inlet_surface, farfield_inflow_surface, 
					farfield_sideline_surfaces_index_array[], farfield_outflow_surface};

surf1 = news;
Surface Loop(surf1) = {unstructured_Volume_surfaces_array[]};
volume_index = news;
Volume(volume_index) = {surf1};

Physical Surface("nozzle_inlet") = {nozzle_inlet_surface};
Physical Surface("lip_surface") = {lip_surfaces_index_array[]};
Physical Surface("farfield") = {farfield_inflow_surface, farfield_sideline_surfaces_index_array[], farfield_outflow_surface};
Physical Surface("nozzle_slip_surfaces") = {nozzle_outside_wall_surfaces_index_array[]};
Physical Surface("nozzle_BL_wall_surfaces") = {nozzle_BL_wall_surfaces_index_array[]};
Physical Surface("nozzle_entry_wall_surfaces") = {nozzle_BL_wall_surfaces_index_array[]};

Physical Volume("domain") = {nozzle_BL_volumes_index_array[], volume_index};
