nx = 10; ny = 10;
nit=200; 

E = 3;
nu = .3;

A = E/(1 - nu^2);
G = E/(2*(1+nu));


dx = 2/(nx-1); dy = 2/(ny-1);

x=0:dx:2; y=0:dy:2;
[Y,X] = meshgrid(y,x);


u=zeros(nx,ny);un=zeros(nx,ny);
v=zeros(nx,ny);vn=zeros(nx,ny);

C = 1/(2*A/dx^2 + 2*G/dy^2);
D1 = 1/(2*A/dy^2 + 2*G/dx^2);

for iit=1:nit   % iterate for solution
    for i=2:nx-1
        for j=2:ny-1
            un(i,j) = C*(A*(u(i+1,j) + u(i-1,j))/dx^2 ...
                        + G*(u(i,j+1) + u(i,j-1))/dy^2 ...
                        + (G+A*nu)*(v(i+1,j+1) - v(i-1,j+1) ...
                             - v(i+1,j-1) + v(i-1,j-1))/(4*dx*dy));
            vn(i,j) = D1*(A*(v(i,j+1) + v(i,j-1))/dy^2 ...
                        + G*(v(i+1,j) + v(i-1,j))/dx^2 ...
                        + (G+A*nu)*(u(i+1,j+1) - u(i-1,j+1) ...
                        - u(i+1,j-1) + u(i-1,j-1))/(4*dx*dy));
        end
    end

    
    
    % apply downward force across top, y stress = A*dv/dy, so dv/dy = stress/A
    vn(:,ny) = vn(:,ny-1)-1*dy/A;
    
   
    
    un(:,ny) = un(:,ny-1);  % this isn't strictly kosher as the vertical strain 
                            % produces horizontal stain ...
    
    % natural boundary conditios for sides
    
    %vn(1,2:ny-1) = vn(2,2:ny-1);
    %vn(nx,2:ny-1) = vn(nx-1,2:ny-1);
    
    % du/dy = -dv/dx
    % du/dx =-nu*dv/dy, (u(2,j) - u(1,j))/dx = -nu*(v(1,j+1)-v(1,j))/dy
    
    for k=2:ny-1
        un(1,k) = un(2,k) + nu*(vn(1,k+1)-vn(1,k-1))/2;
        un(nx,k) = un(nx-1,k) - nu*(vn(nx,k+1)-vn(nx,k-1))/2;
        vn(1,k) = vn(2,k);% + (un(1,k+1) - un(1,k-1))/2;
        vn(nx,k) = vn(nx-1,k);% - (un(10,k+1) - un(10,k-1))/2;
        
    end
    
    u=un;
    v=vn;
        
end
                            

u1=u;
v1=v;

um=zeros(nx,ny);
vm=zeros(nx,ny);

for i=1:nx
    for j=1:ny
        um(i,j) = (i-1)*dx+u(i,j);
        vm(i,j) = (j-1)*dy+v(i,j);
    end
end

plot(0,0)
hold on
        
for j=1:ny
    plot(um(:,j),vm(:,j))
end

for i=1:nx
    plot(um(i,:),vm(i,:))
end

hold off





if 1 == 0
for i=1:nx
    for j=1:ny
        um(i,j) = (i-1)*dx+u(i,j);
        vm(i,j) = (j-1)*dy+v(i,j);
    end
end



plot(0,0)
hold on
        
for i=1:nx
    for j=1:ny
        %if ~((i==6 && j== 6) ||(i==6 && j== 7) ||(i==7 && j== 6) ||(i==7 && j== 7))
           plot(um(i,j),vm(i,j))
        %end
    end
end

for i=1:nx
    for j=1:ny
        %if ~((i==6 && j== 6) ||(i==6 && j== 7) ||(i==7 && j== 6) ||(i==7 && j== 7))
           plot(um(i,j),vm(i,j))
        %end
    end
end

hold off
end

if 1 == 1
% caculate stress and strain
xstrain = zeros(nx,ny);
ystrain = zeros(nx,ny);
shearstrain = zeros(nx,ny);
xstress = zeros(nx,ny);
ystress = zeros(nx,ny);
shearstress = zeros(nx,ny);


for i=1:nx
    for j=1:ny
        
        % calculate x strain, shearstrain on right boundary
        if i==nx 
        % right boundary
            xstrain(i,j) = xstrain(i-1,j);
            % calculate shearstrain on right boundary
            if j==ny
               % top right corner
               shearstrain(i,j) = (u(i,j)-u(i,j-1))/dy + (v(i,j)-v(i-1,j))/dx;
               % right boundary not top
            else   shearstrain(i,j) = (u(i,j+1)-u(i,j))/dy + (v(i,j)-v(i-1,j))/dx; 
            end  
        % not right boundary
        else xstrain(i,j) = (u(i+1,j)-u(i,j))/dx;          
        end
        
       % calculate y strain, shearstrain on top except right boundary             
        if j==ny 
        % on top
            ystrain(i,ny) = ystrain(i,ny-1);
            % calculate shearstrain on top not right boundary
            if i~=nx shearstrain(i,ny) = (u(i,ny)-u(i,ny-1))/dy + (v(i+1,ny)-v(i,ny))/dx;
            end 
        % not on top   
        else ystrain(i,j) = (v(i,j+1)-v(i,j))/dy;
        end
        
       % calculate everywhere but top and right boundary      
        if i~=nx && j~=ny shearstrain(i,j) = (u(i,j+1)-u(i,j))/dy + (v(i+1,j)-v(i,j))/dx;
        end
            
        xstress(i,j) = A*xstrain(i,j) + A*nu*ystrain(i,j);
        ystress(i,j) = A*ystrain(i,j) + A*nu*xstrain(i,j);
        shearstress(i,j) = G*shearstrain(i,j);
    end
end

figure(2)
    surf(xstrain')
    xlabel('x')
    ylabel('y')
    title('X strain')
    figure(3)
    surf(ystrain')
    xlabel('x')
    ylabel('y')
    title('Y strain')
    figure(4)
    surf(shearstrain')
    xlabel('x')
    ylabel('y')
    title('Shear strain')
    figure(5)
    surf(xstress')
    xlabel('x')
    ylabel('y')
    title('X stress')
    figure(6)
    surf(ystress')
    xlabel('x')
    ylabel('y')
    title('Y stress')
    figure(7)
    surf(shearstress')
    xlabel('x')
    ylabel('y')
    title('Shear stress')
    

end