% esempio di soluzione in termini di spostamenti di problema mono-assiale
% nella deformazione con condizioni al contorno cinematiche
% modificato 14/11/2006

L = 100;    % lunghezza dominio Z=[0,L]
N = 50;     % numero punti discretizzazione

E = 210000; % modulo elastico
ni = 0.3;   % coeff. di Poisson
E = E*(1-ni)/((1+ni)*(1-2*ni));

z = linspace(0,L,N)';           % ascissa discretizzata
Dz = L / (N - 1);               % intervallo di discretizzazione


w0 = 0;                         % cc in z=0
wL = 3;                         % cc in z=L
%b = zeros(size(z));             % forze di mazza nulle
%b = ones(size(z))*1000;        % forze di mazza costanti lungo z
b = linspace(1,1000,N);        % forze di mazza lineari lungo z


w1 = zeros(N,1);                % soluzione di tentativo
w1(1) = w0;
w1(N) = wL;

MaxLoop = 10000;                % numero massimo iterazioni
Toll = 1e-7;                    % tolleranza soluzione
w2 = w1;                        % inizializza soluzione aggiornata
for k=1:MaxLoop
    for j = 2:N-1               % calcola soluzione aggiornata imponendo 
                                % l'equazione Ew"+bz=0 in tutti i punti 
                                % interni del dominio
        w2(j) = w1(j+1)/2 + w1(j-1)/2 + b(j) * Dz^2 / 2 / E;
    end
    e = norm(w2 - w1);          % scarto fra sol. di tentativo e 
                                % sol. aggiornata 
    if e < Toll * norm(w2)
        break                   % se tolleranza soddisfatta esce da loop
    end
    w1 = w2;                    % nuova soluzione di tentativo
end
k
w = w2;

s = 5;                          % fattore di scala per spostamento
figure(5)
plot(z,zeros(size(z)),'.b', z + w*s, ones(size(z)),'+r')
xlabel('z')
ylabel('x')
grid on
ylim([-10 10])

figure(6)
plot(z,w)
xlabel('z')
ylabel('w(z)')
grid on

ez = zeros(size(z)) * nan;
for j=2:N-1     % deformazione lungo z ez = w'
    ez(j) = (w(j+1) - w(j-1))/2/Dz;
end
sz = E*ez;      % tensione lungo z  sz = E*ez
figure(7)
plot(z,sz)
xlabel('z')
ylabel('\sigma_z(z)')
grid on


s0 = ni/(1-ni)*sz;    % tensione trasversale   s0 = ni/(1-ni)*sz
figure(8)
plot(z,s0)
xlabel('z')
ylabel('\sigma_0(z)')
grid on
