%Script prova metodo GFDM
%Equazione conduzione calore - stazionario - 2D
%Predefined connectivity
%Scritto da Pavan Andrea - 02/08/2022
clear;
clc;


%% dati input
l1 = 1;        %lunghezza lato x (m)
l2 = 0.1;        %lunghezza lato y (m)
nEdgePoints1 = 30;       %numero punti lato x
nEdgePoints2 = 10;       %numero punti lato y
lmax = 0.02;        %distanza massima punti
searchdist = lmax/1000;      %distanza ricerca punti
Nminsat = 8;        %numero minimo satelliti
dt = 2;      %passo temporale (s)

uL = 400;       %temperatura bordo sx (K)
uR = 300;       %temperatura bordo dx (K)
kcost = 0.1;        %conducibilità termica (W/(m*K))
rho = 8.0;      %densità (kg/m3)
c = 3.8;        %calore specifico (J/(kg*K))

Q = @(x,y,t) 0;     %sorgente
u0 = @(x,y) 0;      %condizioni iniziali
t0 = 0;     %istante iniziale (s)
tf = 120;     %istante finale (s)
ue = @(x,y,t) uL+x*(uR-uL)/l1;      %soluzione esatta


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

%indici punti con condizioni contorno
idxDirichletL = find(P(1:lbo,1)==0);
idxDirichletR = find(P(1:lbo,1)==l1);
idxNeumann = find(and(or(P(1:lbo,2)==0,P(1:lbo,2)==l2),and(P(1:lbo,1)~=0,P(1:lbo,1)~=l1)));


%% 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
x = [];     %distanza x satelliti
y = [];     %distanza y satelliti
for i=1:length(P)
    x{i} = zeros(Nsat(i),1);
    y{i} = zeros(Nsat(i),1);
    for j=1:Nsat(i)
        x{i}(j) = P(Psatidx{i}(j),1)-P(i,1);
        y{i}(j) = P(Psatidx{i}(j),2)-P(i,2);
    end
end

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


%% inversione matrici minimi quadrati
invA = [];     %matrici minimi quadrati invertite
B = [];     %matrice decomposizione termine noto
C = [];     %matrice coefficienti derivate
for j=1:length(idxDirichletL)
    %nodo contorno sx con condizione Dirichlet uL
    invA{i} = [];
    B{i} = [];
    C{i} = [];
end
for j=1:length(idxDirichletR)
    %nodo contorno dx con condizione Dirichlet uR
    invA{i} = [];
    B{i} = [];
    C{i} = [];
end
for j=1:length(idxNeumann)
    %nodo contorno con condizione Neumann dudy=0
    i = idxNeumann(j);
    A = [sum(w2{i}.*x{i}.^2), sum(w2{i}.*x{i}.^3), sum(w2{i}.*x{i}.*y{i}.^2), sum(w2{i}.*x{i}.^2.*y{i});
        sum(w2{i}.*x{i}.^3), sum(w2{i}.*x{i}.^4), sum(w2{i}.*x{i}.^2.*y{i}.^2), sum(w2{i}.*x{i}.^3.*y{i});
        sum(w2{i}.*x{i}.*y{i}.^2), sum(w2{i}.*x{i}.^2.*y{i}.^2), sum(w2{i}.*y{i}.^4), sum(w2{i}.*x{i}.*y{i}.^3);
        sum(w2{i}.*x{i}.^2.*y{i}), sum(w2{i}.*x{i}.^3.*y{i}), sum(w2{i}.*x{i}.*y{i}.^3), sum(w2{i}.*x{i}.^2.*y{i}.^2)];
    invA{i} = inv(A);
    B{i} = zeros(4, 2+Nsat(i));
    B{i}(:,1) = [-sum(w2{i}.*x{i}); -sum(w2{i}.*x{i}.^2); -sum(w2{i}.*y{i}.^2); -sum(w2{i}.*x{i}.*y{i})];
    B{i}(1,2:end-1) = w2{i}.*x{i};
    B{i}(2,2:end-1) = w2{i}.*x{i}.^2;
    B{i}(3,2:end-1) = w2{i}.*y{i}.^2;
    B{i}(4,2:end-1) = w2{i}.*x{i}.*y{i};
    B{i}(:,end) = -sum(w2{i}.*y{i})*ones(4,1);
    C{i} = invA{i}*B{i};
end
for i=lbo+1:length(P)
    %nodo interno
    A = [sum(x{i}.^2.*w2{i}), sum(x{i}.*y{i}.*w2{i}), sum(x{i}.^3.*w2{i}), sum(x{i}.*y{i}.^2.*w2{i}), sum(x{i}.^2.*y{i}.*w2{i});
        sum(x{i}.*y{i}.*w2{i}), sum(y{i}.^2.*w2{i}), sum(x{i}.^2.*y{i}.*w2{i}), sum(y{i}.^3.*w2{i}), sum(x{i}.*y{i}.^2.*w2{i});
        sum(x{i}.^3.*w2{i}), sum(x{i}.^2.*y{i}.*w2{i}), sum(x{i}.^4.*w2{i}), sum(x{i}.^2.*y{i}.^2.*w2{i}), sum(x{i}.^3.*y{i}.*w2{i});
        sum(x{i}.*y{i}.^2.*w2{i}), sum(y{i}.^3.*w2{i}), sum(x{i}.^2.*y{i}.^2.*w2{i}), sum(y{i}.^4.*w2{i}), sum(x{i}.*y{i}.^3.*w2{i});
        sum(x{i}.^2.*y{i}.*w2{i}), sum(x{i}.*y{i}.^2.*w2{i}), sum(x{i}.^3.*y{i}.*w2{i}), sum(x{i}.*y{i}.^3.*w2{i}), sum(x{i}.^2.*y{i}.^2.*w2{i})];
    invA{i} = inv(A);
    B{i} = zeros(5,1+Nsat(i));
    B{i}(:,1) = [-sum(x{i}.*w2{i}); -sum(y{i}.*w2{i}); -sum(x{i}.^2.*w2{i}); -sum(y{i}.^2.*w2{i}); -sum(x{i}.*y{i}.*w2{i})];
    B{i}(1,2:end) = x{i}.*w2{i};
    B{i}(2,2:end) = y{i}.*w2{i};
    B{i}(3,2:end) = x{i}.^2.*w2{i};
    B{i}(4,2:end) = y{i}.^2.*w2{i};
    B{i}(5,2:end) = x{i}.*y{i}.*w2{i};
    C{i} = invA{i}*B{i};
end


%% evoluzione temporale
t = t0:dt:tf;       %istanti temporali [s]
u = zeros(length(P),1);     %soluzione numerica
u(:) = u0(P(:,1),P(:,2));
uprev = u;      %soluzione al passo precedente
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 j=1:length(idxDirichletL)
    i = idxDirichletL(j);
    M(i,i) = 1;
end
for j=1:length(idxDirichletR)
    i = idxDirichletR(j);
    M(i,i) = 1;
end
for j=1:length(idxNeumann)
    i = idxNeumann(j);
    M(i,i) = 1 - (kcost*dt/rho*c)*(2*C{i}(2,1) + 2*C{i}(3,1));
    for js=1:Nsat(i)
        M(i,Psatidx{i}(js)) = -(kcost*dt/rho*c)*(2*C{i}(2,1+js) + 2*C{i}(3,1+js));
    end
end
for i=1+lbo:length(P)
    M(i,i) = 1 - (kcost*dt/rho*c)*(2*C{i}(3,1) + 2*C{i}(4,1));
    for j=1:Nsat(i)
        M(i,Psatidx{i}(j)) = -(kcost*dt/rho*c)*(2*C{i}(3,1+j) + 2*C{i}(4,1+j));
    end
end
%spy(M)
%condest(M)

%propagazione
for l=2:length(t)
    %propagazione metodo Eulero Implicito
    uprev = u;
    p = zeros(length(P),1);
%     for i=1:lbo
%         p(i) = uB(P(i,1),P(i,2),t(l));
%     end
    for j=1:length(idxDirichletL)
        i = idxDirichletL(j);
        p(i) = uL;
    end
    for j=1:length(idxDirichletR)
        i = idxDirichletR(j);
        p(i) = uR;
    end
    for j=1:length(idxNeumann)
        i = idxNeumann(j);
        p(i) = uprev(i);
    end
    for i=1+lbo:length(P)
        p(i) = uprev(i);
    end
    u = M\p;

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

    %grafico evoluzione
    if mod(l,1)==0
        figure(1);
        tiledlayout(1,2);
        nexttile;
        scatter(P(:,1),P(:,2),[],u,'filled');
        title(['Soluzione numerica (t=' num2str(t(l)) 's)']);
        xlabel('x');
        ylabel('y');
        zlabel('Soluzione u(x,y)');
        axis equal;
        colorbar;
        colormap('hot');

        nexttile;
        semilogy(t(1:l),maxerru(1:l),'k-');
        title('Errore soluzione')
        xlabel('Tempo t (s)');
        ylabel('Errore soluzione |u/ue-1|');
    end

end

