Lx = 500;   % dimensioni dominio rettangolare
Ly = 100;

Nx = 101;   % discretizzazione
Ny = 21;

ni = 0.3;   % moduli elastici
E = 210000;

q = 10;     % carico



Dx = Lx / (Nx - 1);
Dy = Ly / (Ny - 1);

x = ((1:Nx)'-1)*Dx * ones(1,Ny);    % coordinate nodi dominio
y = ones(Nx,1) * ((1:Ny)-1)*Dy;

psi = zeros(size(x));               % inizializza funzione di airy


% imposta condizioni al contorno
% AB

% AD
psi(2,:) = Dx * (Lx/4 * q - Lx/4 * q * cos(pi/Ly * y(2,:)));

% BC
psi(Nx-1,:) = Dx * (Lx/4 * q - Lx/4 * q * cos(pi/Ly * y(2,:)));

% CD
psi(:,Ny) = Lx/2 * q * x(:,Ny) - 1/2 * q * x(:,Ny).^2;
psi(:,Ny-1) = psi(:,Ny);
% psi([1 2 Nx-1 Nx],:) = 1;
% psi(:,[1 2 Ny-1 Ny]) = 1;


f = psi;    % funzione tentativo

e = 1;    % errore (distanza fra soluzione tentativo e soluzione corretta) 
toll = 1e-6;    % tolleranza    
while e > toll
    for i=3:Nx-2
        for j=3:Ny-2
            psi(i,j) = -1/(6*Dy^4 + 6*Dx^4 + 8*Dx^2*Dy^2) * ...
                       (Dy^4 * (psi(i+2,j) - 4*psi(i+1,j) ...
                          - 4*psi(i-1,j) + psi(i-2,j)) ...
                        + Dx^4 * (psi(i,j+2) - 4*psi(i,j+1) ...
                          - 4*psi(i,j-1) + psi(i,j-2)) ...
                        + 2*Dx^2*Dy^2 * (psi(i+1,j+1) - 2*psi(i,j+1) ...
                          + psi(i-1,j+1) - 2*psi(i+1,j) - 2*psi(i-1,j) ...
                          + psi(i+1,j-1) - 2*psi(i,j-1) + psi(i-1,j-1)));
        end
    end
    F = reshape(psi,[Nx*Ny, 1]);
    F1 = reshape(f,[Nx*Ny, 1]);
    e = norm(F-F1)/norm(F)
    f = psi;
end

% calcola derivate parziali di psi
[py px] = gradient(psi);
[pxy pxx] = gradient(px);
[pyy pyx] = gradient(py);

% tensioni
sx = pyy/Dy^2;
sy = pxx/Dx^2;
txy = -pxy/Dx/Dy;


figure(1)
surf(x,y,psi)
title('\psi')
xlabel('x')
ylabel('y')

figure(2)
surf(x,y,zeros(size(x)),sx)
title('\sigma_x')
view(2)
axis equal
colorbar
xlim([0 Lx])
ylim([0 Ly])

figure(3)
surf(x,y,zeros(size(x)),sy)
title('\sigma_y')
view(2)
axis equal
colorbar
xlim([0 Lx])
ylim([0 Ly])

figure(4)
surf(x,y,zeros(size(x)),txy)
title('\tau_{xy}')
view(2)
axis equal
colorbar
xlim([0 Lx])
ylim([0 Ly])


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(5)
surf(x,y,zeros(size(x)),s1)
title('\sigma_I')
view(2)
colorbar
axis equal
xlim([0 Lx])
ylim([0 Ly])

figure(6)
surf(x,y,zeros(size(x)),s2)
title('\sigma_{II}')
view(2)
colorbar
axis equal
xlim([0 Lx])
ylim([0 Ly])

sscale = .1;
figure(7)
mesh(x,y,zeros(size(x)))
hold on
quiver(x, y, fi1(:,:,1).*s1*sscale,fi1(:,:,2).*s1*sscale, 0, 'color', 'b')
quiver(x, y, -fi1(:,:,1).*s1*sscale,-fi1(:,:,2).*s1*sscale, 0, 'color', 'b')
quiver(x, y, fi2(:,:,1).*s2*sscale,fi2(:,:,2).*s2*sscale, 0, 'color', 'r')
quiver(x, y, fi2(:,:,1).*-s2*sscale,fi2(:,:,2).*-s2*sscale, 0, 'color', 'r')
hold off
view(2)
axis equal
xlim([0 Lx])
ylim([0 Ly])
title('Direzioni principali tensione')

figure(8)
Lz = 10;
J = (Ly^3 * Lz) / 12;
M = 1/8 * q*Lz * Lx^2;
sx_sv = M/J * (Ly/2 - y(1,:));
plot(sx(floor(Nx/2),:),y(1,:), sx_sv,y(1,:))
xlabel('\sigma_x(L_x/2, y)')
ylabel('y')
title('confronto con soluzione Saint Venant')

