% ex06_naca_airfoil.m — NACA0012 airfoil in a circular far field.
%
% The domain is a large circle with the airfoil removed (ddiff). The size
% function refines toward the leading and trailing edges; fixed points
% pin the airfoil tip and a few far-field points. This is the heaviest
% example — give it a few seconds.

mip load --install distmesh

hlead=0.01; htrail=0.04; hmax=2; circx=2; circr=4;
a=.12/.2*[0.2969,-0.1260,-0.3516,0.2843,-0.1036];

fd=@(p) ddiff(dcircle(p,circx,0,circr),(abs(p(:,2))-polyval([a(5:-1:2),0],p(:,1))).^2-a(1)^2*p(:,1));
fh=@(p) min(min(hlead+0.3*dcircle(p,0,0,0),htrail+0.3*dcircle(p,1,0,0)),hmax);

fixx=1-htrail*cumsum(1.3.^(0:4)');
fixy=a(1)*sqrt(fixx)+polyval([a(5:-1:2),0],fixx);
fix=[[circx+[-1,1,0,0]*circr; 0,0,circr*[-1,1]]'; 0,0; 1,0; fixx,fixy; fixx,-fixy];
box=[circx-circr,-circr; circx+circr,circr];
h0=min([hlead,htrail,hmax]);

[p,t]=distmesh2d(fd,fh,h0,box,fix);

title('NACA0012 airfoil');
fprintf('mesh: %d nodes, %d triangles\n', size(p,1), size(t,1));