Nx = 21;
Ny = 21;
bx = 100;
by = 100;

nie = 0.3;
Ee = 210000;

% ni = nie / (1-nie);
% E = Ee / (1-ni^2);
ni = nie;
E = Ee;

Dx = bx / (Nx - 1);
Dy = by / (Ny - 1);

u = zeros(Nx,Ny);
v = zeros(Nx,Ny);
x = ((1:Nx)'-1)*Dx * ones(1,Ny);
y = ones(Nx,1) * ((1:Ny)-1)*Dy;

% definizione cc
a = 0.001;
b = 0.0005;
c = 0.0002;

% lato inferiore
u(:,1) = 0;
v(:,1) = 0;

% lato sx
u(1,:) = 0;
v(1,:) = a*y(1,:);

% lato dx
u(Nx,:) = b * y(Nx,:);
v(Nx,:) = c * y(Nx,:);

% lato sup
u(:,Ny) = b*by/bx * x(:,Ny);
v(:,Ny) = a*by + (c-a) * by/bx * x(:,Ny);

u1 = u;
v1 = v;

e = 1;
toll = 1e-3;
while e > toll
    for i=2:Nx-1
        for j=2:Ny-1
            u(i,j) = (1/Dx^2 * (u1(i+1,j) + u1(i-1,j)) + ...
                      (1-ni)/2/Dy^2 * (u1(i,j+1) + u1(i,j-1)) + ...
                      (1+ni)/8/Dx/Dy * (v1(i+1,j+1) + v1(i-1,j-1) - v1(i+1,j-1) - v1(i-1,j+1))) / ...
                      (2/Dx^2 + (1-ni)/Dy^2);

            v(i,j) = (1/Dy^2 * (v1(i,j+1) + v1(i,j-1)) + ...
                      (1-ni)/2/Dx^2 * (v1(i+1,j) + v1(i-1,j)) + ...
                      (1+ni)/8/Dx/Dy * (u1(i+1,j+1) + u1(i-1,j-1) - u1(i+1,j-1) - u1(i-1,j+1))) / ...
                      (2/Dy^2 + (1-ni)/Dx^2);
        end
    end
    U = [reshape(u,[Nx*Ny, 1]); reshape(v,[Nx*Ny, 1])];
    U1 = [reshape(u1,[Nx*Ny, 1]); reshape(v1,[Nx*Ny, 1])];
    e = norm(U-U1)/norm(U);
    u1 = u;
    v1 = v;
end

[uy ux] = gradient(u);
[vy vx] = gradient(v);
ux = ux / Dx;
uy = uy / Dy;
vx = vx / Dx;
vy = vy / Dy;


sx = E/(1-ni^2) * (ux + ni*vy);
sy = E/(1-ni^2) * (vy + ni*ux);
txy = E/2/(1+ni) * (uy + vx);

dscale = 100;

figure(1)
mesh(x+u*dscale,y+v*dscale,zeros(size(x)),'edgecolor','r')
hold on
mesh(x,y,zeros(size(x)),'edgecolor','b')
hold off
title('u')
xlabel('x')
ylabel('y')
view(2)
% 
figure(2)
surf(x+u*dscale,y+v*dscale,zeros(size(x)),sx)
title('\sigma_x')
view(2)
colorbar
xlabel('x')
ylabel('y')
axis equal

figure(3)
surf(x+u*dscale,y+v*dscale,zeros(size(x)),sy)
title('\sigma_y')
view(2)
colorbar
xlabel('x')
ylabel('y')
axis equal

figure(4)
surf(x+u*dscale,y+v*dscale,zeros(size(x)),txy)
title('\tau_x')
view(2)
colorbar
xlabel('x')
ylabel('y')
axis equal

% figure(5)
% sz = ni/(1+ni) * (sx + sy);
% surf(x+u*dscale,y+v*dscale,zeros(size(x)),sz)
% title('\sigma_z')
% view(2)
% colorbar
% axis equal

figure(5)
ez = -ni/E * (sx + sy);
surf(x+u*dscale,y+v*dscale,zeros(size(x)),ez)
title('\epsilon_z')
view(2)
colorbar
xlabel('x')
ylabel('y')
axis equal

s1 = zeros(size(x));
s2 = zeros(size(x));
fi1 = zeros([size(x) 2]);
fi2 = zeros([size(x) 2]);
for i=1:Nx
    for j=1:Ny
        S = [sx(i,j), txy(i,j); txy(i,j), sy(i,j)];
        [V,L] = eig(S);
        s2(i,j) = L(1,1);
        s1(i,j) = L(2,2);
        fi2(i,j,:) = V(:,1);
        fi1(i,j,:) = V(:,2);
    end
end
figure(6)
surf(x+u*dscale,y+v*dscale,zeros(size(x)),s1)
title('\sigma_I')
view(2)
colorbar
axis equal

figure(7)
surf(x+u*dscale,y+v*dscale,zeros(size(x)),s2)
title('\sigma_{II}')
view(2)
colorbar
axis equal

sscale = .5;
figure(8)
mesh(x+u*dscale,y+v*dscale,zeros(size(x)))
hold on
quiver(x+u*dscale, y+v*dscale, fi1(:,:,1).*s1*sscale,(:,:,2).*s1*sscale, sscale, 'color', 'b')
quiver(x+u*dscale, y+v*dscale, -fi1(:,:,1).*s1*sscale,-fi1(:,:,2).*s1*sscale, sscale, 'color', 'b')
quiver(x+u*dscale, y+v*dscale, fi2(:,:,1).*s2*sscale,fi2(:,:,2).*s2*sscale, sscale, 'color', 'r')
quiver(x+u*dscale, y+v*dscale, fi2(:,:,1).*-s2*sscale,fi2(:,:,2).*-s2*sscale, sscale, 'color', 'r')
hold off
view(2)
axis equal


