% Schematic3.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; 
c = 0.7;
cx = [((x(7)+2*x(8))/3) ((x(7)+2*x(8))/3) ((x(9)+2*x(8))/3) ((x(9)+2*x(8))/3)];
cy = [y(5)+h/6 y(3)-h/6 y(3)-h/6 y(5)+h/6];
patch(cx, cy, [c c c],'edgecolor','none')
cx = [x(7)-h/6 x(7)-h/6 x(9)+h/6 x(9)+h/6];
cy = [((y(3)+2*y(4))/3) ((y(5)+2*y(4))/3) ((y(5)+2*y(4))/3) ((y(3)+2*y(4))/3)];
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
text('Interpreter','latex', 'String','$$\Delta_h \phi_{ijk}^n$$','position',[0.51 0.25 0],'FontSize',26)
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_points1')