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


%% dati input
l1 = 1;        %lunghezza lato x (m)
l2 = 1;        %lunghezza lato y (m)
nEdgePoints = 30;       %numero punti in ciascun lato
lmax = 0.01;        %distanza massima punti
searchdist = lmax/1000;      %distanza ricerca punti
Nminsat = 8;        %numero minimo satelliti
dt = 0.1;      %passo temporale (s)

S = @(x,y,t) 0;     %sorgente
u0 = @(x,y) 0;      %condizioni iniziali
t0 = 0;     %istante iniziale (s)
tf = 0.8;        %istante finale (s)
ue = @(x,y,t) exp(8*x).*sin(8*y)/exp(8);      %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',lmax);
P = [boundaryNodes; P];
lbo = length(boundaryNodes);


%% popolamento stelle
Psatidx = [];        %elenco punti satellite
Nsat = zeros(length(P),1);      %numero punti satellite
for i=1:length(P)
    %cerco i satelliti del punto i
    searchdisti = searchdist;
    while Nsat(i)<Nminsat
        %Psatidx{i} = find((abs(P(:,1)-P(i,1))<=searchdisti).*(abs(P(:,2)-P(i,2))<=searchdisti).*((P(:,1)-P(i,1)).^2+(P(:,2)-P(i,2)).^2>2*eps));
        Psatidx{i} = find(((P(:,1)-P(i,1)).^2+(P(:,2)-P(i,2)).^2<=searchdisti).*((P(:,1)-P(i,1)).^2+(P(:,2)-P(i,2)).^2>2*eps));
        Nsat(i) = length(Psatidx{i});
        searchdisti = searchdisti*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 = h{i}.^2+k{i}.^2;
    R2max = max(R2);
    w2{i} = exp(-1*R2/R2max).^2;
    %w2{i} = 1 - 6*(R2/R2max).^2 + 8*(R2/R2max).^3 - 3*(R2/R2max).^4;
end


%% inversione matrici minimi quadrati
invA = [];     %matrici minimi quadrati invertite
B = [];     %matrice decomposizione termine noto
C = [];     %matrice coefficienti derivate
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);
    B{i} = zeros(5,1+Nsat(i));
    B{i}(:,1) = [-sum(h{i}.*w2{i}); -sum(k{i}.*w2{i}); -0.5*sum(h{i}.^2.*w2{i}); -0.5*sum(k{i}.^2.*w2{i}); -sum(h{i}.*k{i}.*w2{i})];
    B{i}(1,2:end) = h{i}.*w2{i};
    B{i}(2,2:end) = k{i}.*w2{i};
    B{i}(3,2:end) = 0.5*h{i}.^2.*w2{i};
    B{i}(4,2:end) = 0.5*k{i}.^2.*w2{i};
    B{i}(5,2:end) = h{i}.*k{i}.*w2{i};
    C{i} = invA{i}*B{i};
end


%% evoluzione temporale
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 = [];     %errore assoluto soluzione
maxerru = 0*t;      %errore max soluzione

%matrice Eulero Implicito
M = sparse(length(P),length(P));
for i=1:lbo
    M(i,i) = 1;
end
for i=1+lbo:length(P)
    M(i,i) = 1 - dt*C{i}(3,1) - dt*C{i}(4,1);
    for j=1:Nsat(i)
        M(i,Psatidx{i}(j)) = - dt*C{i}(3,1+j) - dt*C{i}(4,1+j);
    end
end
% M = sparse(length(P)+lbo,length(P));
% for i=1+lbo:length(P)
%     M(i,i) = 1 - dt*C{i}(3,1) - dt*C{i}(4,1);
%     for j=1:Nsat(i)
%         M(i,Psatidx{i}(j)) = - dt*C{i}(3,1+j) - dt*C{i}(4,1+j);
%     end
% end
% for i=1+length(P):length(P)+lbo
%     %condizioni contorno
%     M(i,i-length(P)) = 1;
% end
%spy(M)
%condest(M)

%propagazione
for l=2:length(t)
    %propagazione metodo Eulero Implicito
    p = zeros(length(P),1);
    for i=1:lbo
        p(i) = uB(P(i,1),P(i,2),t(l));
    end
    for i=1+lbo:length(P)
        p(i) = u(i,l-1);
    end
%     p = zeros(length(P)+lbo,1);
%     for i=1:length(P)
%         p(i) = u(i,l-1);
%     end
%     for i=1+length(P):lbo+length(P)
%         p(i) = uB(P(i-length(P),1),P(i-length(P),2),t(l));
%     end
    u(:,l) = M\p;

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

    %grafico evoluzione
    if mod(l,1)==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),maxerru(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),erru{end},'b.');
title('Distribuzione errore');
xlabel('x');
ylabel('y');
zlabel('Errore assoluto err(u)');

%grafico connettività punto con errore massimo
plotidx = find(abs(erru{end})==maxerru(end));
plotidx = plotidx(1);
figure();
plot(P(:,1),P(:,2),'k.');
hold on;
plot(P(plotidx,1),P(plotidx,2),'ro');
linecolor = [rand rand rand];
for j=1:length(Psatidx{plotidx})
    plot([P(plotidx,1),P(Psatidx{plotidx}(j),1)],[P(plotidx,2),P(Psatidx{plotidx}(j),2)],'-','Color',linecolor);
end
hold off;

