Inverse transform for quadrilateral elements

inv_quad.m

[x y]=ginput(4);
x=100*x;
y=100*y;
subplot(1,2,1);plot(x([1:4 1]),y([1:4 1]))
hold on
subplot(1,2,2);plot([-1 -1 1 1 -1],[-1 1 1 -1 -1]);
hold on
subplot(1,2,1);title('Hit enter to exit')
while 1
    p1=[];p2=[];
    p1=ginput(1);
    
    if isempty(p1)
        break
    end
    p2=ginput(2);
    p(:,1)=p1(1,1)+(p2(1,1)-p1(1,1)).*[0:0.1:1]';
    p(:,2)=p1(1,2)+(p2(1,2)-p1(1,2)).*[0:0.1:1]';
    subplot(1,2,1);plot(p(:,1),p(:,2),'*')
    for i=1:length(p)
        u(i,1)=subs(S.u(1),[x1 x2 x3 x4 y1 y2 y3 y4 xp yp],[x' y' p(i,1) p(i,2)]);
        w(i,1)=subs(S.w(1),[x1 x2 x3 x4 y1 y2 y3 y4 xp yp],[x' y' p(i,1) p(i,2)]);
    end
    subplot(1,2,2);plot(u,w,'*')
end