lc=50;
length=1500;
dec=1000;

Merge "cylindercone.stl";

s1 = news ;
s2 = s1 -1 ;
sl1 = newsl ; Surface Loop (sl1) = {s2} ;
//v1=newv ; Volume(v1) = {sl1};

Transfinite Surface {s2}; Recombine Surface {s2};



// ############ boundary layer creation ############
BLlayers = 8;
BLgrowth = 1.2;
cells[] = {1,1};
heights = 1/10;
heights[1] = heights * (BLgrowth+1);
For i In {2:BLlayers-1}
// Printf("new %g",i);
cells[i] = 1;
heights[i] = heights[i-1] + BLgrowth * (heights[i-1] -heights[i-2]);
EndFor

v2[] = Extrude{Surface{s2};Layers{cells[],heights[]};Recombine;};

// ############ boundary layer creation ############



sl2 = newsl ; Surface Loop (sl2) = {-v2[0]} ;
Transfinite Surface "*" ; 
Transfinite Volume "*" ; 


// External Domain
xmin = -length + dec;
xmax =  length ; 
ymin = -length + dec;
ymax =  length - dec;
zmin = -length + dec;
zmax =  length - dec; 

p1 = newp ; Point (p1) = { xmin, ymin, zmin ,lc }; p2 = newp ; Point (p2) = { xmax, ymin, zmin ,lc };
p3 = newp ; Point (p3) = { xmax, ymax, zmin ,lc }; p4 = newp ; Point (p4) = { xmin, ymax, zmin ,lc };
p5 = newp ; Point (p5) = { xmin, ymin, zmax ,lc }; p6 = newp ; Point (p6) = { xmax, ymin, zmax ,lc };
p7 = newp ; Point (p7) = { xmax, ymax, zmax ,lc }; p8 = newp ; Point (p8) = { xmin, ymax, zmax ,lc };
l1 = newl ; Line(l1)   = {p1, p4}; l2 = newl ; Line(l2)   = {p4, p8}; l3 = newl ; Line(l3)   = {p8, p5}; l4 = newl ; Line(l4)   = {p5, p1};
l5 = newl ; Line(l5)   = {p2, p3}; l6 = newl ; Line(l6)   = {p3, p7}; l7 = newl ; Line(l7)   = {p7, p6}; l8 = newl ; Line(l8)   = {p6, p2};
l9 = newl ; Line(l9)   = {p1, p2}; l10 = newl ; Line(l10) = {p4, p3}; l11 = newl ; Line(l11) = {p5, p6}; l12 = newl ; Line(l12) = {p8, p7};
ll1 = newll ; Line Loop (ll1) = { l1,l2,l3,l4};     s3 = news ; Plane Surface (s3) ={ll1} ;
ll2 = newll ; Line Loop (ll2) = { l5,l6,l7,l8};     s4 = news ; Plane Surface (s4) ={ll2} ;
ll3 = newll ; Line Loop (ll3) = { -l3,l12,l7,-l11}; s5 = news ; Plane Surface (s5) ={ll3} ;
ll4 = newll ; Line Loop (ll4) = { l1,l10,-l5,-l9};  s6 = news ; Plane Surface (s6) ={ll4} ;
ll5 = newll ; Line Loop (ll5) = { l4,l9,-l8,-l11};  s7 = news ; Plane Surface (s7) ={ll5} ;
ll6 = newll ; Line Loop (ll6) = { -l2,l10,l6,-l12}; s8 = news ; Plane Surface (s8) ={ll6} ;
sl3 = newsl ; Surface Loop (sl3) = {s3,s4,s5,s6,s7,s8} ;

v3=newv ; Volume(v3) = {sl3,sl2};

Physical Surface("body") = {s2};
Physical Surface("xmin") = {s3};
Physical Surface("xmax") = {s4};
Physical Surface("ymin") = {s5};
Physical Surface("ymax") = {s6};
Physical Surface("zmin") = {s7};
Physical Surface("zmax") = {s8};
Physical Volume (v4) = { v3, v2[1]} ;
