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


%% dati input
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) -exp(2*x+y).*(rho*c*sin(t)+5*kcost*cos(t));       %sorgenti calore [W/m3]
u0 = @(x,y) exp(2*x+y);     %condizioni iniziali [K]
uB = @(x,y,t) exp(2*x+y)*cos(t);        %condizioni contorno [K]
t0 = 0;     %istante iniziale [s]
tf = 10;     %istante finale [s]
ue = @(x,y,t) exp(2*x+y).*cos(t);       %soluzione esatta [K]


%% generazione pointcloud
l1 = 0.1;        %lunghezza lato x [m]
l2 = 0.1;        %lunghezza lato y [m]
nEdgePoints = 10;       %numero punti in ciascun lato
%nEdgePoints = 15;       %numero punti in ciascun lato
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);
P = [boundaryNodes; P];

%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{i} = [];
for i=1:length(P)
    w2{i} = exp(-500*(h{i}.^2+k{i}.^2)).^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
dt = 5e-4;      %passo temporale [s]
%dt = 5e-5;      %passo temporale [s]
t = t0:dt:tf;       %istanti temporali [s]
u = zeros(length(P),length(t));     %soluzione numerica
u(:,1) = u0(P(:,1),P(:,2));
erru = 0*t;     %errore assoluto soluzione

params.kcost = kcost;
params.rho = rho;
params.c = c;
params.Q = Q;
params.P = P;
params.h = h;
params.k = k;
params.Psatidx = Psatidx;
params.invA = invA;
params.lbo = length(boundaryNodes);
params.w2 = w2;

for l=2:length(t)
    params.t = t(l);

    %propagazione metodo Eulero Esplicito
    u(:,l) = u(:,l-1) + dt*f(u(:,l-1),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.');
        scatter3(P(:,1),P(:,2),u(:,l),[],abs(u(:,l)),'filled');
        title(['Soluzione numerica (t=' num2str(t(l)) 's)']);
        xlabel('x');
        ylabel('y');
        zlabel('Soluzione u(x,y)');
        zlim([-2 2]);
        colorbar;

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

end


%% discretizzazione spaziale
function res = f(u,params)
kcost = params.kcost;
rho = params.rho;
c = params.c;
Q = params.Q;
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 = (kcost/(rho*c))*(d2u_dx2 + d2u_dy2) + Q(P(:,1),P(:,2),t)/(rho*c);
end

