% Calculapte strain at interior points
xstrain = zeros(nx,ny);
ystrain = zeros(nx,ny);
shearstrain = zeros(nx,ny);
xstress = zeros(nx,ny);
ystress = zeros(nx,ny);
shearstress = zeros(nx,ny);

ssint =@(i,j) (u(i,j+1)-u(i,j-1))/(2*dy) + (v(i+1,j)-v(i-1,j))/(2*dx);
sstopnc =@(i,j) (u(i,j)-u(i,j-1))/dy + (v(i+1,j)-v(i-1,j))/(2*dx); 
ssrbnc =@(i,j) (u(i,j+1)-u(i,j-1))/(2*dy) + (v(i,j)-v(i-1,j))/dx;
ssbotnc =@(i,j) (u(i,j+1)-u(i,j))/dy + (v(i+1,j)-v(i-1,j))/(2*dx); 
sslbnc =@(i,j) (u(i,j+1)-u(i,j-1))/(2*dy) + (v(i+1,j)-v(i,j))/dx;

sstlc =@(i,j) (u(i,j)-u(i,j-1))/dy + (v(i+1,j)-v(i,j))/dx;
sstrc =@(i,j) (u(i,j)-u(i,j-1))/dy + (v(i,j)-v(i-1,j))/dx;
ssblc =@(i,j) (u(i,j+1)-u(i,j))/dy + (v(i+1,j)-v(i,j))/dx;
ssbrc =@(i,j) (u(i,j+1)-u(i,j))/dy + (v(i,j)-v(i-1,j))/dx;


for i=1:nx
    for j=1:ny
      
        if i==1  
            % left boundary
            if j==1
                % bottom left corner
                xstrain(1,1) = (u(2,1) - u(1,1))/dx;
                ystrain(1,1) = (v(1,2) - v(1,1))/dy;
                shearstrain(1,1) = ssblc(1,1);
                x=shearstrain(1,1)
                
            elseif j==ny              
                % top left corner
                xstrain(1,ny) = (u(2,ny) - u(1,ny))/dx;
                ystrain(1,ny) = (v(1,ny) - v(1,ny-1))/dy;
                shearstrain(1,ny) = sstlc(1,ny);
            
            else
                % left boundary, not a corner       
                xstrain(1,j) = (u(2,j) - u(1,j))/dx;
                ystrain(1,j) = (v(1,j+1) - v(1,j-1))/(2*dy);
                shearstrain(1,j) = sslbnc(1,j);
            end
            
       elseif i==nx
            % right boundary
            if j==1
                % bottom right corner
                xstrain(nx,1) = (u(nx,1) - u(nx-1,1))/dx;
                ystrain(nx,1) = (v(nx,2) - v(nx,1))/dy;
                shearstrain(nx,1) = ssbrc(nx,1);
                
            elseif j==ny
                % top right corner
                xstrain(nx,ny) = (u(nx,ny) - u(nx-1,ny))/dx;
                ystrain(nx,ny) = (v(nx,ny) - v(nx,ny-1))/dy;
                shearstrain(nx,ny) = sstrc(nx,ny);
                
            else
                % right boundary, not a corner
                xstrain(nx,j) = (u(nx,j) - u(nx-1,j))/dx;
                ystrain(nx,j) = (v(nx,j+1) - v(nx,j-1))/(2*dy);
                shearstrain(nx,j) = ssrbnc(nx,j);
            end
            
       elseif j==1
            % bottom boundary, not a corner
             xstrain(i,1) = (u(i+1,1) - u(i-1,1))/(2*dx);
             ystrain(i,1) = (v(i,2) - v(i,1))/dy;
             shearstrain(i,1) = ssbotnc(i,1);
             
       elseif j==ny
            % top boundary, not a corner
             xstrain(i,ny) = (u(i+1,ny) - u(i-1,ny))/(2*dx);
             ystrain(i,ny) = (v(i,ny) - v(i,ny-1))/dy;
             shearstrain(i,ny) = sstopnc(i,ny);
             
       else
           % interior node
           xstrain(i,j) = (u(i+1,j) - u(i-1,j))/(2*dx);
           ystrain(i,j) = (v(i,j+1) - v(i,j-1))/(2*dy);
           shearstrain(i,j) = ssint(i,j);
       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
                
               
               
      
 

if (1 == 1)
 


    
    shearstress(5,6) = G*ssrbnc(5,6);
    shearstress(5,7) = G*ssrbnc(5,7);
    shearstress(8,6) = G*sslbnc(8,6);
    shearstress(8,7) = G*sslbnc(8,7);

    shearstress(6,5) = G*ssbotnc(6,5);
    shearstress(7,5) = G*ssbotnc(7,5);
    shearstress(6,8) = G*sstopnc(6,8);
    shearstress(7,8) = G*sstopnc(6,8);


      
    shearstress(6,6)=NaN;
    shearstress(6,7)=NaN;
    shearstress(7,6)=NaN;
    shearstress(7,7)=NaN;
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')