disp('%%%%%%%')
clear all; clf;

% parameters
Nx=10;           % number of grid points in x
Ny=20;            % number of grid points in y
p=-1;         % desired slope at top boundary
L=0.95;      % the length in y of the computational domain

% differentiation
[d.x,d.xx,d.wx,x]=dif1D('cheb',0,1,Nx);
[d.y,d.yy,d.wy,y]=dif1D('cheb',0,L,Ny);
[D,l,X,Y,Z,I,NN]=dif2D(d,x,y);
D0=D; Y0=Y;

ZZ=blkdiag(Z,zeros(Nx,Nx));
II=blkdiag(I,eye(Nx,Nx));

l.h=(1:NN)';
l.eta=NN+[1:Nx]';
l.top=[l.ctl; l.top; l.ctr];
l.bot=[l.cbl; l.bot; l.cbr];

%initial guess
eta0=zeros(Nx,1);
sol0=[Y(:).*(1-Y(:));eta0];
sol=sol0;

% Newton iterations
quit=0;count=0;
while ~quit

    % the present solution and its derivatives
    h=sol(l.h); hy=D.y*h;  hxx=D.xx*h; hyy=D.yy*h; 
    eta=sol(l.eta);

    % addapt the domain
    H=reshape(h,Ny,Nx); Hy=reshape(hy,Ny,Nx);
    mesh(X,Y,H,'edgecolor','b');hold on
    
    for gre=1:Nx
        ll=Y(Ny,gre);
        y=Y(:,gre);
        Y(:,gre)=Y(:,gre)*(1+eta(gre)/ll); % stretched grid
        yy=linspace(1.0001*ll,2*ll,100)'; % the grid outside the domain
        hh=H(Ny,gre)+(yy-ll)*Hy(Ny,gre); % linear extrapolation outside of the domain
        H(:,gre)=interp1([y; yy],[H(:,gre); hh],Y(:,gre),'splines'); % interpolation to the new grid
    end
    D=map2D(X,Y,D0); % the new differentiation matrices
    sol0=[Y(:).*(1-Y(:));eta0]; % update sol0 for the boundary conditions
    sol(l.h)=H(:); % update h
    sol(l.eta)=0; % set eta to zero

    mesh(X,Y,H,'edgecolor','r'); hold off; 
    xlabel('x');ylabel('y'); zlabel('h'); title('before and after interpolation')
    drawnow
    
    % the present solution and its derivatives
    h=sol(l.h); hy=D.y*h;  hxx=D.xx*h; hyy=D.yy*h; 
    eta=sol(l.eta);
    
    % nonlinear function
    f=[hxx+hyy+2; hy(l.top)+eta.*hyy(l.top)-p*ones(Nx,1)]; 

    % analytical jacobian
    A=[D.xx+D.yy, ZZ(l.h,l.eta); ...
       D.y(l.top,:)+diag(eta)*D.yy(l.top,:), diag(hyy(l.top))];

    % Boundary conditions
    loc=[l.left; l.right; l.bot; l.top];
    C=II([l.bot;l.right;l.left],:);
    
    f(loc)=[C*(sol-sol0); ... % the linear boundary conditions
            h(l.top)+diag(eta)*hy(l.top)]; % the nonlinear boundary conditions
    A(loc,:)=[C; ...
              II(l.top,:)+diag(eta)*D.y(l.top,:)*II(l.h,:)+diag(hy(l.top))*II(l.eta,:)];
    
    
    % convergence test
    res=norm(f,inf);
    disp([num2str(count) '  ' num2str(res)]);
    if count>50|res>1e5; disp('no convergence'); break; end
    if res<1e-9; quit=1; disp('converged'); continue; end

    % Newton step
    sol=sol-A\f;
    count=count+1;
end

set(gcf,'paperpositionmode','auto')
print('-dpng','-r80','domain_derivative_2D_adapt.png')

% validation
htheo=Y.*(1-Y);
err=norm(h-htheo(:),2)