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


%% dati input
a = 0.1;        %lunghezza lato x [m]
b = 0.1;        %lunghezza lato y [m]
nEdgePoints = 10;       %numero punti in ciascun lato
searchdist = a/10;      %distanza ricerca punti
Nminsat = 8;        %numero minimo satelliti


%% generazione pointcloud
boundaryNodes = [linspace(0,a,nEdgePoints)', 0*ones(nEdgePoints,1);
    a+0*ones(nEdgePoints,1), linspace(0,b,nEdgePoints)';
    fliplr(linspace(0,a,nEdgePoints))', b+0*ones(nEdgePoints,1);
    0*ones(nEdgePoints,1), fliplr(linspace(0,b,nEdgePoints))'];
boundaryNodes = unique(boundaryNodes,'rows','stable');      %punti contorno
P = generate_pointcloud(boundaryNodes,'nextrasteps',0,'lmax',0.1);
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
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));
        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


%% grafico connettività stelle
figure();
plot(P(:,1),P(:,2),'k.','MarkerSize',10);
hold on;
for i=1:10:length(P)
    linecolor = [rand rand rand];
    for j=1:length(Psatidx{i})
        plot([P(i,1),P(Psatidx{i}(j),1)],[P(i,2),P(Psatidx{i}(j),2)],'-','Color',linecolor);
    end
end
title('Esempio connettività pointcloud');
xlabel('x [m]');
ylabel('y [m]');
axis equal;
axis square;
hold off;


