%Script prova metodo GFDM - flusso a potenziale attorno NACA-0012
%Equazione Laplace - stazionario - 2D
%Predefined connectivity
%Scritto da Pavan Andrea - 16/07/2022
clear;
clc;

%% dati input
airfoilNpoints = 200;
boundaryRadius = 20;     %raggio bordo esterno
boundaryNpoints = 50;      %numero punti bordo esterno
lmax = 3;        %distanza massima punti
searchdist = 1e-9;      %distanza ricerca punti
Nminsat = 8;        %numero minimo satelliti

u0 = @(x,y) 1*x;      %condizioni iniziali
t0 = 0;     %istante iniziale (s)
tf = 1000;        %istante finale (s)
dt = 100;      %passo temporale (s)


%% generazione geometria
%nodi profilo
boundaryNodes{1} = naca4_generator('4412',airfoilNpoints);
airfoilNormal = 0*boundaryNodes{1};     %versori normali entranti
airfoilNormal(1,:) = [-1 0];
for i=2:airfoilNpoints
    t = boundaryNodes{1}(i,:)-boundaryNodes{1}(i-1,:);
    li = norm(t);
    t = t/li;
    airfoilNormal(i,:) = [-t(2), t(1)];
end

%nodi contorno
boundaryNodes{2}(:,1) = 0.25 + boundaryRadius*cos(0:(2*pi/boundaryNpoints):2*pi);
boundaryNodes{2}(:,2) = boundaryRadius*sin(0:(2*pi/boundaryNpoints):2*pi);
boundaryNodes{2}(end,:) = [];


%% generazione pointcloud
P = generate_pointcloud_2(boundaryNodes,'boolInternal',[false true],'nextrasteps',1,'lmax',lmax);
P = [boundaryNodes{1}; boundaryNodes{2}; P];
lbo = length(boundaryNodes{1})+length(boundaryNodes{2});


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

%matrice Eulero Implicito
M = sparse(length(P),length(P));
for i=1:airfoilNpoints      %condizione impermeabilità parete profilo
    M(i,i) = C{i}(1,1)*airfoilNormal(i,1) + C{i}(2,1)*airfoilNormal(i,2);
    for j=1:Nsat(i)
        M(i,Psatidx{i}(j)) = C{i}(1,1+j)*airfoilNormal(i,1) + C{i}(2,1+j)*airfoilNormal(i,2);
    end
end
for i=1+airfoilNpoints:lbo      %condizione farfield
    M(i,i) = C{i}(1,1);
    for j=1:Nsat(i)
        M(i,Psatidx{i}(j)) = C{i}(1,1+j);
    end
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
%spy(M)
%condest(M)

%propagazione
for l=2:length(t)
    %propagazione metodo Eulero Implicito
    p = zeros(length(P),1);
    for i=1:airfoilNpoints
        p(i) = 0;
    end
    for i=1+airfoilNpoints:lbo
        p(i) = 1;
    end
    for i=1+lbo:length(P)
        p(i) = u(i,l-1);
    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)');
    end

end

%% confronto con metodo pannelli
cp = zeros(length(P),1);
Vx = zeros(length(P),1);
Vy = zeros(length(P),1);
for i=1:length(P)
    Vx(i) = C{i}(1,:)*u([i;Psatidx{i}],end);
    Vy(i) = C{i}(2,:)*u([i;Psatidx{i}],end);
    cp(i) = 1-(Vx(i)^2+Vy(i)^2)^2;
end

[~,~,cp_panel,~,~,~] = panel1(boundaryNodes{1},0);
figure();
plot(cp_panel(:,1),cp_panel(:,2),'b-');
hold on;
plot(P(1:airfoilNpoints,1),cp(1:airfoilNpoints),'ro');
title('Confronto distribuzione pressione');
legend('Metodo a pannelli','GFDM');
xlabel('x/c');
ylabel('cp');
set(gca,'Ydir','reverse');
ylim([-2 1]);
hold off;

