% Schematic4.m
clear; clf; hold; N=15;h=1/N; 
for i=1:N
    x(i)=(i)*h;
end
y=x; x2=linspace(0,1,N+1);y2=x2; figure(1); hold;
for i=1:N
    for j=1:N
        plot(x(i),y(j),'.', 'Markersize', 28)
    end
end
    c = 0.7;
    cx = [x(8)-h/8 x(8)-h/8 x(9)+h/8 x(9)+h/8];
    cy = [y(4)+h/8 y(3)-h/8 y(3)-h/8 y(4)+h/8];
    patch(cx, cy, [c c c],'edgecolor','none')
for j=1:N+1
    line([x2(1) x2(end)],[y2(j) y2(j)],'color','k')
    line([x2(j) x2(j)],[y2(1) y2(end)],'color','k')
end
for i=1:N
    for j=1:N
        d(i,j)=0.27-sqrt((x(i)-0.5)^2+(y(j)-0.5)^2);
        if (abs(d(i,j))<1.01*sqrt(2)*h)
            u(i,j)=0;
        else
            u(i,j)=1;
        end
    end
end
hold on
v=0*u;
for i=2:N-1
    for j=2:N-1
      if ((abs(u(i+1,j)-u(i-1,j))+abs(u(i,j+1)-u(i,j-1)))*u(i,j)>0.5)
            v(i,j)=1000;
      else
            v(i,j)=0;
      end
    end
end
dx=0*d; dy=0*d; figure(1); hold on; s=0.5*h; dd=-abs(d);
for i=2:N-1
    for j=2:N-1
        if (v(i,j)>0.5)
            plot(x(i),y(j),'ko', 'Markersize', 20)
            dx(i,j)=(dd(i+1,j)-dd(i-1,j))/(2*h);
            dy(i,j)=(dd(i,j+1)-dd(i,j-1))/(2*h);
        end
    end
end
    quiver(x(8), y(5), -dd(8,5)*dx(8,5), -dd(8,5)*dy(8,5), 0,... 
        'linewidth',3, 'Color',[0 0 0])
    plot(x(8)-dd(8,5)*dx(8,5), y(5)-dd(8,5)*dy(8,5), 'kd',... 
        'linewidth',2, 'Markersize', 27)
    text('Interpreter','latex', 'String','$${\bf x}_{ijk}$$',...
        'position',[0.51 0.35 0],'FontSize',20)
    text('Interpreter','latex', 'String','$${\rm cp}({\bf x}_{ijk})$$','position',[0.53 0.21 0],'FontSize',20)
for i=2:N-1
    for j=2:N-1
        if (u(i,j)==0)
            plot(x(i),y(j),'k.', 'Markersize', 53)
        end
    end
end
[xx,yy]=meshgrid(x,y); contour(xx,yy,d,[0 0],'k', 'linewidth', 2)
axis image; axis([x(6) x(10) y(2) y(6)]); axis off; view(0,90);
print('-deps', 'schematic_points2')
