%Script prova metodo GFDM
%Equazione Laplace - stazionario - 2D
%Predefined connectivity
%Scritto da Pavan Andrea - 09/07/2022
clear;
clc;


%% dati input
l1 = 1;        %lunghezza lato x (m)
l2 = 1;        %lunghezza lato y (m)
nEdgePoints = 15;       %numero punti in ciascun lato
dt = 1e-3;      %passo temporale (s)

S = @(x,y,t) -6*x -6*y +4*(x.^2+y.^2-1).*exp(-x.^2-y.^2);     %sorgente
u0 = @(x,y) 1;      %condizioni iniziali
t0 = 0;     %istante iniziale (s)
tf = 0.8;       %istante finale (s)
ue = @(x,y,t) -x.^3 -y.^3 +exp(-x.^2 -y.^2);      %soluzione esatta
uB = @(x,y,t) ue(x,y,t);       %condizioni contorno


%% generazione pointcloud
boundaryNodes = [linspace(0,l1,nEdgePoints)', 0*ones(nEdgePoints,1);
    l1+0*ones(nEdgePoints,1), linspace(0,l2,nEdgePoints)';
    fliplr(linspace(0,l1,nEdgePoints))', l2+0*ones(nEdgePoints,1);
    0*ones(nEdgePoints,1), fliplr(linspace(0,l2,nEdgePoints))'];
boundaryNodes = unique(boundaryNodes,'rows','stable');      %punti contorno
P = generate_pointcloud(boundaryNodes,'nextrasteps',0,'lmax',0.2);
P = [boundaryNodes; P];
lbo = length(boundaryNodes);

%plot pointcloud
figure(1);
plot(P(:,1),P(:,2),'k.');
hold on;
plot(boundaryNodes(:,1),boundaryNodes(:,2),'r.');
title('Pointcloud iniziale')
xlabel('x (m)');
ylabel('y (m)');

axis equal;
axis square;
hold off;


%% popolamento stelle
Nminsat = 8;        %numero minimo satelliti
Psatidx = [];        %elenco punti satellite
Nsat = zeros(length(P),1);      %numero punti satellite
for i=1:length(P)
    %cerco i satelliti del punto i
    searchdist = l1/10;      %distanza ricerca punti
    while Nsat(i)<Nminsat
        %Psatidx{i} = find((abs(P(:,1)-P(i,1))<=searchdist).*(abs(P(:,2)-P(i,2))<=searchdist));
        Psatidx{i} = find((abs(P(:,1)-P(i,1))<=searchdist).*(abs(P(:,2)-P(i,2))<=searchdist).*((P(:,1)-P(i,1)).^2+(P(:,2)-P(i,2)).^2>2*eps));
        Nsat(i) = length(Psatidx{i});
        searchdist = searchdist*1.1;
    end
end

%calcolo distanze satelliti
h = [];     %distanza x satelliti
k = [];     %distanza y satelliti
for i=1:length(P)
    h{i} = zeros(Nsat(i),1);
    k{i} = zeros(Nsat(i),1);
    for j=1:Nsat(i)
        h{i}(j) = P(Psatidx{i}(j),1)-P(i,1);
        k{i}(j) = P(Psatidx{i}(j),2)-P(i,2);
    end
end

%calcolo pesi satelliti
w2 = [];
for i=1:length(P)
    R2 = max(h{i}.^2+k{i}.^2);
    w2{i} = exp(-1*(h{i}.^2+k{i}.^2)/R2).^2;
end


%% inversione matrici minimi quadrati
invA = [];     %matrici minimi quadrati invertite
for i=1:length(P)
    A = [sum(h{i}.^2.*w2{i}), sum(h{i}.*k{i}.*w2{i}), 0.5*sum(h{i}.^3.*w2{i}), 0.5*sum(h{i}.*k{i}.^2.*w2{i}), sum(h{i}.^2.*k{i}.*w2{i});
        sum(h{i}.*k{i}.*w2{i}), sum(k{i}.^2.*w2{i}), 0.5*sum(h{i}.^2.*k{i}.*w2{i}), 0.5*sum(k{i}.^3.*w2{i}), sum(h{i}.*k{i}.^2.*w2{i});
        0.5*sum(h{i}.^3.*w2{i}), 0.5*sum(h{i}.^2.*k{i}.*w2{i}), 0.25*sum(h{i}.^4.*w2{i}), 0.25*sum(h{i}.^2.*k{i}.^2.*w2{i}), 0.5*sum(h{i}.^3.*k{i}.*w2{i});
        0.5*sum(h{i}.*k{i}.^2.*w2{i}), 0.5*sum(k{i}.^3.*w2{i}), 0.25*sum(h{i}.^2.*k{i}.^2.*w2{i}), 0.25*sum(k{i}.^4.*w2{i}), 0.5*sum(h{i}.*k{i}.^3.*w2{i});
        sum(h{i}.^2.*k{i}.*w2{i}), sum(h{i}.*k{i}.^2.*w2{i}), 0.5*sum(h{i}.^3.*k{i}.*w2{i}), 0.5*sum(h{i}.*k{i}.^3.*w2{i}), sum(h{i}.^2.*k{i}.^2.*w2{i})];
    invA{i} = inv(A);
end


%% evoluzione temporale
params.S = S;
params.P = P;
params.h = h;
params.k = k;
params.Psatidx = Psatidx;
params.invA = invA;
params.lbo = lbo;
params.w2 = w2;
params.t = t0;

t = t0:dt:tf;       %istanti temporali [s]
u = zeros(length(P),length(t));     %soluzione numerica
u(:,1) = u0(P(:,1),P(:,2));
du = f(u(:,1),params);
erru = 0*t;     %errore assoluto soluzione

for l=2:length(t)
%     %propagazione metodo Eulero Esplicito
%     u(:,l) = u(:,l-1) + dt*f(u(:,l-1),params);

    %propagazione metodo RK4 Esplicito
    k1 = du;
    params.t = t(l-1)+0.5*dt;
    k2 = f(u(:,l-1)+0.5*k1*dt, params);
    params.t = t(l-1)+0.5*dt;
    k3 = f(u(:,l-1)+0.5*k2*dt, params);
    params.t = t(l-1)+dt;
    k4 = f(u(:,l-1)+k3*dt, params);
    u(:,l) = u(:,l-1) + (k1+2*k2+2*k3+k4)*dt/6;
    du = f(u(:,l), params);

    %condizioni contorno Dirichlet
    u(1:length(boundaryNodes),l) = uB(P(1:length(boundaryNodes),1),P(1:length(boundaryNodes),2),t(l));

    %calcolo errore e controllo divergenza metodo
    erru(l) = max(abs(u(:,l)-ue(P(:,1),P(:,2),t(l))));
    if max(abs(u(:,l)))>1e10
        error('Errore: la soluzione ha superato la soglia 10^10');
    end

    %grafico evoluzione
    if mod(l,100)==0
        figure(1);
        plot3(P(:,1),P(:,2),u(:,l),'k.');
        title(['Soluzione numerica (t=' num2str(t(l)) 's)']);
        xlabel('x');
        ylabel('y');
        zlabel('Soluzione u(x,y)');
        %zlim([-2 2]);

        figure(2);
        semilogy(t(1:l),erru(1:l),'k-');
        title('Errore soluzione')
        xlabel('Tempo t (s)');
        ylabel('Errore assoluto soluzione |u-ue|');
        ylim([1e-5 1]);
    end

end

%% confronto con soluzione analitica
[meshX,meshY] = meshgrid(linspace(0,l1,nEdgePoints),linspace(0,l2,nEdgePoints));
figure();
surf(meshX,meshY,ue(meshX,meshY));
hold on;
plot3(P(:,1),P(:,2),u(:,l),'k.');
title('Confronto soluzione esatta');
xlabel('x');
ylabel('y');
zlabel('u(x,y)');
colorbar;
hold off;

%distribuzione errore
figure();
plot3(P(:,1),P(:,2),abs(u(:,end)-ue(P(:,1),P(:,2),t(end))),'b.');
title('Distribuzione errore');
xlabel('x');
ylabel('y');
zlabel('Soluzione u(x,y)');


%% discretizzazione spaziale
function res = f(u,params)
S = params.S;
P = params.P;
h = params.h;
k = params.k;
invA = params.invA;
lbo = params.lbo;
Psatidx = params.Psatidx;
w2 = params.w2;
t = params.t;

d2u_dx2 = 0*u;      %derivate numeriche
d2u_dy2 = 0*u;

%calcolo derivate
for i=1+lbo:length(P)
    %calcolo derivate
    b = [sum(h{i}.*w2{i}.*(u(Psatidx{i})-u(i))); sum(k{i}.*w2{i}.*(u(Psatidx{i})-u(i))); sum(0.5*h{i}.^2.*w2{i}.*(u(Psatidx{i})-u(i))); sum(0.5*k{i}.^2.*w2{i}.*(u(Psatidx{i})-u(i))); sum(h{i}.*k{i}.*w2{i}.*(u(Psatidx{i})-u(i)))];
    D = invA{i}*b;
    d2u_dx2(i) = D(3);
    d2u_dy2(i) = D(4);
end

%du_dt = F(t,u)
res = 0*u;
res(1+lbo:end) = d2u_dx2(1+lbo:end) + d2u_dy2(1+lbo:end) - S(P(1+lbo:end,1),P(1+lbo:end,2),t);
end

