%Script prova metodo FDM
%Flessione piastra Kirchhoff monostrato isotropa
%Modifiche rispetto prova_03:
%   propagazione temporale II ordine con metodo Leapfrog
%   aggiunto termine smorzamento
%Scritto da Pavan Andrea - 01/03/2022
clear;
clc;


%% dati input
a = 0.5;        %lunghezza lato x [m]
b = 0.5;        %lunghezza lato y [m]
h = 2e-3;       %spessore [m]
E = 73e9;       %modulo Young materiale [Pa]
ni = 0.33;      %coefficiente Poisson materiale
rho = 2700;     %densità materiale [kg/m3]

dx = 0.05;     %discretizzazione lato x
x = 0:dx:a;     %griglia lato x
dy = 0.05;     %discretizzazione lato y
y = 0:dy:b;     %griglia lato y
y = y(:);
[meshX,meshY] = meshgrid(x,y);
t0 = 0;     %istante iniziale [s]
tf = 0.15;     %istante finale [s]

w0 = @(x,y) 0;      %condizioni iniziali spostamento [m]
dw0 = @(x,y) 0;     %condizioni iniziali velocità [m/s]
wB = @(x,y) 0;      %condizioni contorno Dirichlet spostamento [m]
dwBN = @(x,y) 0;    %condizioni contorno Neuman momento [Nm]

qz0 = -1000;        %intensità carico distribuito [Pa]
qz = @(x,y) qz0*sin(pi*x/a).*sin(pi*y/b);       %carico distribuito [Pa]
D = E*(h^3)/(12*(1+ni));        %rigidezza flessionale piastra [Nm]
we = @(x,y) sin(pi*x/a).*sin(pi*y/b).*qz0./(D*((pi/a)^4)+2*D*((pi/a)^2)*((pi/b)^2)+D*((pi/b)^4));       %soluzione esatta [m]

params.dx = dx;
params.dy = dy;
params.rho = rho;
params.h = h;
params.D = D;
params.qz = qz(x,y);

%% evoluzione temporale
dt = 2e-4;      %passo temporale [s]
t = t0:dt:tf;       %istanti temporali [s]
w = zeros(length(x),length(y),length(t));       %soluzione numerica
dw = 0*w;
w(:,:,1) = w0(x,y);
dw(:,:,1) = dw0(x,y);

maxw = 0*t;     %spostamento massimo ad ogni istante
for l=2:length(t)
    w(:,:,l) = 0;     %condizioni contorno Dirichlet spostamento

%     %propagazione metodo Eulero Esplicito
%     dw(:,:,l) = 0.97*dw(:,:,l-1) + dt*f(w(:,:,l-1),params);
%     w(:,:,l) = w(:,:,l-1) + dt*dw(:,:,l);

    %propagazione metodo Leapfrog
    al = f(w(:,:,l-1),params);
    w(:,:,l) = w(:,:,l-1) + dt*dw(:,:,l-1) + 0.5*dt*dt*al;
    dw(:,:,l) = dw(:,:,l-1) + 0.5*dt*(al+f(w(:,:,l),params));
    dw(:,:,l) = 0.97*dw(:,:,l);     %smorzamento

    %plot
    maxw(l) = max(max(abs(w(:,:,l))));
    figure(3);
    plot(t(1:l),maxw(1:l),'r-*');
    title('Andamento temporale spostamento massimo');
    xlabel('t');
    ylabel('max|w|');
end

%% plot
figure(2);
surf(meshX,meshY,1000*we(meshX,meshY));
title('Soluzione esatta w(x,y)');
xlabel('x');
ylabel('y');
zlabel('w [mm]');
zlim([-10 0]);
colorbar;

maxwe = max(max(abs(we(meshX,meshY))));
figure(4);
semilogy(t(1:l),abs(maxw(1:l)-maxwe)/maxwe,'k-*');
title('Andamento errore relativo su spostamento massimo');
xlabel('t');
ylabel('err(max|w|)');


%% discretizzazione spaziale
function res = f(w,params)
dx = params.dx;
dy = params.dy;
rho = params.rho;
h = params.h;
D = params.D;
qz = params.qz;

d4wdx4 = 0*w;        %derivate numeriche
d4wdy4 = 0*w;
d4wdx2dy2 = 0*w;

d4wdx4(3:end-2,:) = (w(1:end-4,:) -4*w(2:end-3,:) +6*w(3:end-2,:) -4*w(4:end-1,:) +w(5:end,:))/(dx^4);
d4wdx4(2,:) = (w(2,:)*60/11 -w(3,:)*48/11 +w(4,:)*12/11)/(dx^4);
d4wdx4(end-1,:) = (w(end-1,:)*60/11 -w(end-2,:)*48/11 +w(end-3,:)*12/11)/(dx^4);
d4wdy4(:,3:end-2) = (w(:,1:end-4) -4*w(:,2:end-3) +6*w(:,3:end-2) -4*w(:,4:end-1) +w(:,5:end))/(dy^4);
d4wdy4(:,2) = (w(:,2)*60/11 -w(:,3)*48/11 +w(:,4)*12/11)/(dy^4);
d4wdy4(:,end-1) = (w(:,end-1)*60/11 -w(:,end-2)*48/11 +w(:,end-3)*12/11)/(dy^4);
d4wdx2dy2(2:end-1,2:end-1) = (w(1:end-2,1:end-2)-2*w(1:end-2,2:end-1)+w(1:end-2,3:end) -2*(w(2:end-1,1:end-2)-2*w(2:end-1,2:end-1)+w(2:end-1,3:end)) +w(3:end,1:end-2)-2*w(3:end,2:end-1)+w(3:end,3:end))/((dy^2)*(dx^2));

res = (qz - D*(d4wdx4 + 2*d4wdx2dy2 + d4wdy4))/(rho*h);
end

