function [p,t]=AdvancingFront(p,f,Nf,DistanzaInterna,DivisioneDistanzaSup,TempoVis,SezioneX)

global TOLLERANZA;
TOLLERANZAPICCOLA=1e-10;
ORDINEM=2;
Tempo=clock;
% FORMULETTE SGUAPPE PER IL CALCOLO DEL CIRCOCENTRO E CIRCORAGGIO DI UN
% TRIANGONO IN 3D DATI I SUOI 3 VERTICI a, b, c
% 
%     |                                                           |
%     | |c-a|^2 (b-a)x(c-a)]x(b-a) + |b-a|^2 (c-a)x[(b-a)x(c-a)   |
%     |                                                           |
% r = -------------------------------------------------------------,
%                          2 | (b-a)x(c-a) |^2
% 
%         |c-a|^2 (b-a)x(c-a)]x(b-a) + |b-a|^2 (c-a)x[(b-a)x(c-a)
% m = a + ---------------------------------------------------------.
%                            2 | (b-a)x(c-a) |^2

% PREPROCESSING: CALCOLO DELLA DISTANZA MEDIA TRA I NODI SUPERFICIALI (NSI= node spacing information)
fprintf('Esecuzione pre-processing...');
NNodiSup=size(p,1);
NSI=zeros(size(p,1),1);
edges=sort([f(:,[1 2]); f(:,[1 2]); f(:,[1 2])],2);
edges=unique(edges,'rows');
L=sqrt(sum((  p(edges(:,2),:)-p(edges(:,1),:) ).^2,2));
for i=1:size(p,1)
    NSI(i)=mean(L(any(edges-repmat([i i],size(edges,1),1)==0,2)));
end
PassoMaxGriglia=min(min(NSI)/DivisioneDistanzaSup,DistanzaInterna);
%calcolo i circocentri dei tringoli di superficie
T1=f(:,1); T2=f(:,2); T3=f(:,3); 
v1=[p(T2,1)-p(T1,1) p(T2,2)-p(T1,2) p(T2,3)-p(T1,3)]; %b-a delle formuline
v2=[p(T3,1)-p(T1,1) p(T3,2)-p(T1,2) p(T3,3)-p(T1,3)]; %c-a delle formuline
v_cross=cross(v1,v2,2);                                 %(b-a)x(c-a) delle formuline
CC=(repmat(sum(v2.^2,2),1,3).*cross(v_cross,v1,2)+repmat(sum(v1.^2,2),1,3).*cross(v2,v_cross,2))./(2*repmat(sum(v_cross.^2,2),1,3));
CircoCentri=CC+[p(T1,1) p(T1,2) p(T1,3)]; 
CircoRaggi=sqrt(sum(CC.^2,2));

fprintf('fatto\n');

% PRIMA PARTE: GENERAZIONE DEI VERTICI INTERNI AL SOLIDO

%Genero una griglia molto fitta che includa il volume da meshare
fprintf('Generazione nodi interni molto densi...');
m=min(p);
M=max(p);
[x,y,z] = ndgrid(m(1):PassoMaxGriglia:M(1),m(2):PassoMaxGriglia:M(2),m(3):PassoMaxGriglia:M(3));
%La ordino in modo che una scansione lineare corrisponda ad una scansione
%spaziale "a strati" partendo dal basso verso l'alto
Griglia=[z(:) y(:) x(:)];
Griglia=sortrows(Griglia);
Griglia(:,[1 2 3])=Griglia(:,[3 2 1]);

%Per ogni punto della griglia calcolo quale tringolo della superficie è il
%+ vicino (distanza punto-triangolo vedere articolo riportato... è un casino vettorizzarlo ma ce l'ho fatta ;) )
App1=p(f(:,1),:)-p(f(:,2),:);
App2=p(f(:,1),:)-p(f(:,3),:);
V1=App1./repmat(sqrt(sum(App1.^2,2)),1,3)+App2./repmat(sqrt(sum(App2.^2,2)),1,3);
App1=p(f(:,2),:)-p(f(:,1),:);
App2=p(f(:,2),:)-p(f(:,3),:);
V2=App1./repmat(sqrt(sum(App1.^2,2)),1,3)+App2./repmat(sqrt(sum(App2.^2,2)),1,3);
App1=p(f(:,3),:)-p(f(:,1),:);
App2=p(f(:,3),:)-p(f(:,2),:);
V3=App1./repmat(sqrt(sum(App1.^2,2)),1,3)+App2./repmat(sqrt(sum(App2.^2,2)),1,3);
GrigliaMod=repmat(Griglia,size(f,1),1);
Proiezioni=GrigliaMod;
App=fliplr(sortrows(repmat(eye(size(f,1)),size(Griglia,1),1)));
%calcolo le proiezioni di ciascun punto della griglia sui piani
%identificati da ciascun triangolo
Proiezioni=Proiezioni-repmat(dot((Proiezioni-App*p(f(:,1),:)),App*Nf,2),1,3).*(App*Nf);
App1=App*p(f(:,1),:);
App2=App*p(f(:,2),:);
App3=App*p(f(:,3),:);
V12=App2-App1;
V23=App3-App2;
V31=App1-App3;
V1X=Proiezioni-App1;
V2X=Proiezioni-App2;
V3X=Proiezioni-App3;
V13=App3-App1;
V21=App1-App2;
V32=App2-App3;
DentroAlTriangolo=zeros(size(Proiezioni,1),1);
%Devo capire quli proiezioni sono cadute all'interno dei corrispondenti triangoli e quali no
DentroAlTriangolo(find( dot(cross(V1X,V12,2),cross(V13,V12,2),2)>=0  & dot(cross(V2X,V23,2),cross(V21,V23,2),2)>=0 & dot(cross(V3X,V31,2),cross(V32,V31,2),2)>=0 ))=1;
%per i punti la cui proiezione cade dentro al triangolo la distanza dal
%triangolo coincide con la distanza tra il punto e la sua proiezione (in realtà calcolo anche le altre perché mi servono per dopo)
Distanze=sqrt(sum((Proiezioni-GrigliaMod).^2,2));
%per i punti la cui proiezione cade fuori dal triangolo, calcolo la
%distanza di tale proiezione dal triangolo (è la minore tra le distanze con ciascun lato)
%e la distanza complessiva si calcola col mitico teorema di pitagora
IndiciFuori=find(DentroAlTriangolo==0);
R1=cross(cross(App2(IndiciFuori,:)-Proiezioni(IndiciFuori,:),App1(IndiciFuori,:)-Proiezioni(IndiciFuori,:),2),App2(IndiciFuori,:)-App1(IndiciFuori,:),2);
Norma=sqrt(sum(R1.^2,2));
Bug=find(Norma==0);
warning off;
R1=R1./repmat(Norma,1,3);
warning on;
R1(Bug,:)=0;
Proiezione1=Proiezioni(IndiciFuori,:)+repmat((dot(App1(IndiciFuori,:)-Proiezioni(IndiciFuori,:),R1,2)),1,3).*R1;
t=sign(dot(Proiezione1-App1(IndiciFuori,:),App2(IndiciFuori,:)-App1(IndiciFuori,:),2)).*sqrt(sum((Proiezione1-App1(IndiciFuori,:)).^2,2))./sqrt(sum((App2(IndiciFuori,:)-App1(IndiciFuori,:)).^2,2));
Distanze1=zeros(size(IndiciFuori,1),1);
I1=find(t>=0 & t<=1); I2=find(t<0); I3=find(t>1);
Distanze1(I1)=sqrt(sum((Proiezione1(I1,:)-Proiezioni(IndiciFuori(I1),:)).^2,2));
Distanze1(I2)=sqrt(sum((App1(IndiciFuori(I2),:)-Proiezioni(IndiciFuori(I2,:),:)).^2,2));
Distanze1(I3)=sqrt(sum((App2(IndiciFuori(I3),:)-Proiezioni(IndiciFuori(I3,:),:)).^2,2)); %Distanze1 contiene la distanza tra la proiezione sul piano ed il segmento 1-2
R2=cross(cross(App3(IndiciFuori,:)-Proiezioni(IndiciFuori,:),App1(IndiciFuori,:)-Proiezioni(IndiciFuori,:),2),App3(IndiciFuori,:)-App1(IndiciFuori,:),2);
Norma=sqrt(sum(R2.^2,2));
Bug=find(Norma==0);
warning off;
R2=R2./repmat(sqrt(sum(R2.^2,2)),1,3);
warning on;
R2(Bug,:)=0;
Proiezione2=Proiezioni(IndiciFuori,:)+repmat((dot(App1(IndiciFuori,:)-Proiezioni(IndiciFuori,:),R2,2)),1,3).*R2;
t=sign(dot(Proiezione2-App1(IndiciFuori,:),App3(IndiciFuori,:)-App1(IndiciFuori,:),2)).*sqrt(sum((Proiezione2-App1(IndiciFuori,:)).^2,2))./sqrt(sum((App3(IndiciFuori,:)-App1(IndiciFuori,:)).^2,2));
Distanze2=zeros(size(IndiciFuori,1),1);
I1=find(t>=0 & t<=1); I2=find(t<0); I3=find(t>1);
Distanze2(I1)=sqrt(sum((Proiezione2(I1,:)-Proiezioni(IndiciFuori(I1),:)).^2,2));
Distanze2(I2)=sqrt(sum((App1(IndiciFuori(I2),:)-Proiezioni(IndiciFuori(I2),:)).^2,2));
Distanze2(I3)=sqrt(sum((App3(IndiciFuori(I3),:)-Proiezioni(IndiciFuori(I3),:)).^2,2)); %Distanze2 contiene la distanza tra la proiezione sul piano ed il segmento 1-3
R3=cross(cross(App3(IndiciFuori,:)-Proiezioni(IndiciFuori,:),App2(IndiciFuori,:)-Proiezioni(IndiciFuori,:),2),App3(IndiciFuori,:)-App2(IndiciFuori,:),2);
Norma=sqrt(sum(R3.^2,2));
Bug=find(Norma==0);
warning off;
R3=R3./repmat(sqrt(sum(R3.^2,2)),1,3);
warning on;
R3(Bug,:)=0;
Proiezione3=Proiezioni(IndiciFuori,:)+repmat((dot(App2(IndiciFuori,:)-Proiezioni(IndiciFuori,:),R3,2)),1,3).*R3;
t=sign(dot(Proiezione3-App2(IndiciFuori,:),App3(IndiciFuori,:)-App2(IndiciFuori,:),2)).*sqrt(sum((Proiezione3-App2(IndiciFuori,:)).^2,2))./sqrt(sum((App3(IndiciFuori,:)-App2(IndiciFuori,:)).^2,2));
Distanze3=zeros(size(IndiciFuori,1),1);
I1=find(t>=0 & t<=1); I2=find(t<0); I3=find(t>1);
Distanze3(I1)=sqrt(sum((Proiezione3(I1,:)-Proiezioni(IndiciFuori(I1),:)).^2,2));
Distanze3(I2)=sqrt(sum((App2(IndiciFuori(I2),:)-Proiezioni(IndiciFuori(I2),:)).^2,2));
Distanze3(I3)=sqrt(sum((App3(IndiciFuori(I3),:)-Proiezioni(IndiciFuori(I3),:)).^2,2)); %Distanze3 contiene la distanza tra la proiezione sul piano ed il segmento 2-3
Distanze(IndiciFuori)=sqrt(Distanze(IndiciFuori).^2+min(min(Distanze1.^2,Distanze2.^2),Distanze3.^2));
%riordino le distanze in una matrice
Distanze=reshape(Distanze,size(Griglia,1),size(f,1));
%Ora devo controllare per TUTTI i triangoli a distanza minima che il punto
%si trovi nella direzione della normale (cioè, appunto, all'interno del volume)
[DistMin TriangoliPiuVicini]= min(Distanze,[],2);
[Point Tri]=find(Distanze<repmat(DistMin+TOLLERANZA,1,size(f,1)));%ipotesi: se le due distanze dei punti dai triangoli differiscono per meno di TOLLERANZA allora sono uguali
[Point m, n]=unique(Point);
IndiciBuoni=DistMin(n)>TOLLERANZA & dot(CircoCentri(Tri,:)-Griglia(Point(n),:),Nf(Tri,:),2)<0;
App=(1:size(Point,1))';
App=histc(App(n).*IndiciBuoni,0:size(Point,1)');
vec=histc(n,1:max(n));
%Cancello i punti della griglia che sono esterni o sulla superficie
Griglia=Griglia(find(App(2:size(App,1))==vec),:);
DistMin=DistMin(find(App(2:size(App,1))==vec));
NNodi=size(Griglia,1);
fprintf('fatto\n');
%Ora devo generare le NSI per ogni punto della griglia così ottenuta
[DistMin Ix]=sort(DistMin);
Griglia=[Griglia(Ix,:); p];
NSI=[zeros(NNodi,1); NSI];
NSI(NNodi)=DistanzaInterna; %il nodo più interno al volume ha NSI=DistanzaInterna cioè il parametro passato alla funzione
fprintf('Generazione delle Node spacing information...');
App2=(NNodi:(NNodi+NNodiSup))';
hProgBar=progressbar( [],0,'Generazione delle Node spacing information, Avanzamento' );
for i=1:(NNodi-1) %si può vettorizzare... ma non ho voglia!
    Trovato=0;
    %Uso solo le NSI sulla superficie e centrale
    %trova il nodo con NSI più vicino a quello corrente (i)
    %[m I1]=min(sum((repmat(Griglia(i,:),NNodiSup+1,1)-Griglia(NNodi:(NNodi+NNodiSup),:)).^2,2));
    [m VettScelte1]=sort(sum((repmat(Griglia(i,:),NNodiSup+1,1)-Griglia(NNodi:(NNodi+NNodiSup),:)).^2,2));
    Scelta1=0;
    while (~Trovato & Scelta1<size(VettScelte1,1))
        Scelta1=Scelta1+1;    
        I1=App2(VettScelte1(Scelta1));
        App=find(dot(repmat(Griglia(I1,:)-Griglia(i,:),NNodiSup+1,1), Griglia(NNodi:(NNodi+NNodiSup),:)-repmat(Griglia(i,:),NNodiSup+1,1),2)<=0);
        App=App2(App);
        App_2=setdiff(App,[I1]);
        App=App_2;
        [m VettScelte2]=sort(sum((repmat(Griglia(i,:),size(App,1),1)-Griglia(App,:)).^2,2));
        Scelta2=0;
        while (~Trovato & Scelta2<size(VettScelte2,1))
            Scelta2=Scelta2+1;
            I2=App_2(VettScelte2(Scelta2));
            A01=Griglia(I1,:)-Griglia(i,:);
            A01=A01./repmat(norm(A01),1,3);
            A02=Griglia(I2,:)-Griglia(i,:);
            A02=A02./repmat(norm(A02),1,3);
            A03=Griglia(NNodi:(NNodi+NNodiSup),:)-repmat(Griglia(i,:),NNodiSup+1,1);
            A03=A03./repmat(sqrt(sum(A03.^2,2)),1,3);
            App=find(sum((A03+repmat(A01+A02,NNodiSup+1,1)).^2,2)<=1);
            App=App2(App);
            App_3=setdiff(App,[I1; I2]);
            App=App_3;
            [m VettScelte3]=sort(sum((repmat(Griglia(i,:),size(App,1),1)-Griglia(App,:)).^2,2));
            Scelta3=0;
            while(~Trovato & Scelta3<size(VettScelte3,1))
                Scelta3=Scelta3+1;
                I3=App_3(VettScelte3(Scelta3));
                V0=dot(cross(repmat(Griglia(I2,:),NNodiSup+1,1)-repmat(Griglia(I1,:),NNodiSup+1,1),repmat(Griglia(I3,:),NNodiSup+1,1)-repmat(Griglia(I1,:),NNodiSup+1,1),2),Griglia(NNodi:(NNodi+NNodiSup),:)-repmat(Griglia(I1,:),NNodiSup+1,1),2);
                V1=dot(cross(repmat(Griglia(I2,:),NNodiSup+1,1)-repmat(Griglia(i,:),NNodiSup+1,1),repmat(Griglia(I3,:),NNodiSup+1,1)-repmat(Griglia(i,:),NNodiSup+1,1),2),Griglia(NNodi:(NNodi+NNodiSup),:)-repmat(Griglia(i,:),NNodiSup+1,1),2);
                V2=dot(cross(repmat(Griglia(i,:),NNodiSup+1,1)-repmat(Griglia(I1,:),NNodiSup+1,1),repmat(Griglia(I3,:),NNodiSup+1,1)-repmat(Griglia(I1,:),NNodiSup+1,1),2),Griglia(NNodi:(NNodi+NNodiSup),:)-repmat(Griglia(I1,:),NNodiSup+1,1),2);
                V3=dot(cross(repmat(Griglia(I2,:),NNodiSup+1,1)-repmat(Griglia(I1,:),NNodiSup+1,1),repmat(Griglia(i,:),NNodiSup+1,1)-repmat(Griglia(I1,:),NNodiSup+1,1),2),Griglia(NNodi:(NNodi+NNodiSup),:)-repmat(Griglia(I1,:),NNodiSup+1,1),2);
                V4=dot(cross(repmat(Griglia(I2,:),NNodiSup+1,1)-repmat(Griglia(I1,:),NNodiSup+1,1),repmat(Griglia(I3,:),NNodiSup+1,1)-repmat(Griglia(I1,:),NNodiSup+1,1),2),repmat(Griglia(i,:)-Griglia(I1,:),NNodiSup+1,1),2);
                App=find(min([V0.*V1, V1.*V2, V2.*V3, V3.*V4],[],2)>=0 & ~V0==0);
                App=App2(App);
                App=setdiff(App,[I1; I2; I3]);
                [m I4]=min(sum((repmat(Griglia(i,:),size(App,1),1)-Griglia(App,:)).^2,2));
                I4=App(I4);
                if (length(I4)>0)
                    Trovato=1;
                end
            end
        end  
    end
    %Interpoliamo le information spacing dai quattro vertici selezionati
    %per ottenere quelle del nodo in questione
    Distanze=sqrt(sum((repmat(Griglia(i,:),4,1)-Griglia([I1;I2;I3;I4],:)).^2,2));
    NSI(i)=sum(NSI([I1;I2;I3;I4])./(Distanze.^ORDINEM) )/sum(1./(Distanze.^ORDINEM));
    hProgBar = progressbar( hProgBar,1/(NNodi-1) );
end
progressbar(hProgBar,-1 );
fprintf('fatto\n');
fprintf('Scrematura nodi interni in base alle node spacing information...');
hProgBar=progressbar( [],0,'Scrematura nodi interni in base alle node spacing information, Avanzamento' );
%Selezioniamo alcuni nodi interni in base alle information spacing
VettSelezionati=[zeros(NNodi,1); ones(NNodiSup,1)];
for i=1:NNodi
   DaConfrontare=find(VettSelezionati);
   if (min(sum((repmat(Griglia(i,:),size(DaConfrontare,1),1)- Griglia(DaConfrontare,:)).^2,2))> NSI(i)^2) 
       VettSelezionati(i)=1;
   end
    hProgBar = progressbar( hProgBar,1/(NNodi) );
end
progressbar(hProgBar,-1 );
Griglia=Griglia(find(VettSelezionati),:);
NSI=NSI(find(VettSelezionati));
fprintf('fatto\n');

% SECONDA PARTE: GENERAZIONE DEGLI EDGE CHE COSTITUISCONO LA MESH

NNodiTot=size(Griglia,1);
Griglia=[Griglia((NNodiTot-NNodiSup+1):NNodiTot,:); Griglia(1:(NNodiTot-NNodiSup),:)];
SuperficieGenerazione=f;
NormaliSupGenerazione=Nf;
ListaNodiSuperficie=(1:NNodiSup)';
ListaNodiInterni=((NNodiSup+1):NNodiTot)';
t=zeros(0,4);
QualitaTet=zeros(0,0);
fprintf('Generazione degli edge della mesh...');
hProgBar=progressbar( [],0,'Costruzione degli edge, Avanzamento');
CostanteProgBar=size(ListaNodiInterni,1);
while (size(SuperficieGenerazione,1)>0)
    %qui dovrei ordinare le superfici in base alla loro area! in modo da
    %usare come base dei tetraedri le facce + piccole...  ma non c'ho voglia!
    FacciaScelta=0;
    TrovatoCandidato=0;
    while(~TrovatoCandidato)
        FacciaScelta=FacciaScelta+1;
        if (FacciaScelta>size(SuperficieGenerazione,1))
            fprintf('Errore fatale... impossibile inserire un nuovo tetraedro');
            return;
        end
        %Seleziono tutti i nodi che possono essere utili a formare il nuovo
        %tetraedro: tutti i nodi che sono ancora interni + tutti i nodi di
        %superficie
        NNodiSup=size(ListaNodiSuperficie,1);
        NNodiInt=size(ListaNodiInterni,1);
        ListaNodiCandidati=[ListaNodiSuperficie; ListaNodiInterni];
        %Diamo una prima scrematura ai candidati selezionando solo quelli
        %che sono dalla parte giusta del triangolo di base (cioè dalla parte della normale)
        %e che, se selezionati, danno luogo ad un tetraedro che non
        %contiene nessun altro punto (nemmeno sulla loro superficie!!)
        V0=squeeze(dot(cross(repmat(Griglia(SuperficieGenerazione(FacciaScelta,2),:),[NNodiSup+NNodiInt,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),repmat(Griglia(SuperficieGenerazione(FacciaScelta,3),:),[NNodiSup+NNodiInt,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),2),repmat(Griglia(ListaNodiCandidati,:),[1,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),2));
        V1=squeeze(dot(cross(repmat(Griglia(SuperficieGenerazione(FacciaScelta,2),:),[NNodiSup+NNodiInt,1,NNodiTot])-repmat(reshape(Griglia',[1,3,NNodiTot]),[NNodiSup+NNodiInt,1,1]),repmat(Griglia(SuperficieGenerazione(FacciaScelta,3),:),[NNodiSup+NNodiInt,1,NNodiTot])-repmat(reshape(Griglia',[1,3,NNodiTot]),[NNodiSup+NNodiInt,1,1]),2),repmat(Griglia(ListaNodiCandidati,:),[1,1,NNodiTot])-repmat(reshape(Griglia',[1,3,NNodiTot]),[NNodiSup+NNodiInt,1,1]),2));
        V2=squeeze(dot(cross(repmat(reshape(Griglia',[1,3,NNodiTot]),[NNodiSup+NNodiInt,1,1])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),repmat(Griglia(SuperficieGenerazione(FacciaScelta,3),:),[NNodiSup+NNodiInt,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),2),repmat(Griglia(ListaNodiCandidati,:),[1,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),2));
        V3=squeeze(dot(cross(repmat(Griglia(SuperficieGenerazione(FacciaScelta,2),:),[NNodiSup+NNodiInt,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),repmat(reshape(Griglia',[1,3,NNodiTot]),[NNodiSup+NNodiInt,1,1])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),2),repmat(Griglia(ListaNodiCandidati,:),[1,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),2));
        V4=squeeze(dot(cross(repmat(Griglia(SuperficieGenerazione(FacciaScelta,2),:),[NNodiSup+NNodiInt,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),repmat(Griglia(SuperficieGenerazione(FacciaScelta,3),:),[NNodiSup+NNodiInt,1,NNodiTot])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),2),repmat(reshape(Griglia',[1,3,NNodiTot]),[NNodiSup+NNodiInt,1,1])-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),[NNodiSup+NNodiInt,1,NNodiTot]),2)); 
        SS=dot( Griglia(ListaNodiCandidati,:)-repmat(Griglia(SuperficieGenerazione(FacciaScelta,1),:),NNodiSup+NNodiInt,1),repmat(NormaliSupGenerazione(FacciaScelta,:),NNodiSup+NNodiInt,1),2);
        Confronto=V0.*V1; Confronto(:,:,2)=V1.*V2; Confronto(:,:,3)=V2.*V3; Confronto(:,:,4)=V3.*V4; 
        Confronto=min(Confronto,[],3);
        App=find( max(Confronto,[],2)<=0 & sum(Confronto==0,2)==4 & SS>0);
        ListaNodiCandidati=ListaNodiCandidati(App);
        %I candidati rimasti vengono ulteriormante selezionati imponendo
        %che il tetraedro che se ne ottiene non intersechi nessun'altra
        %faccia
        CandidatoBuono=zeros(size(ListaNodiCandidati,1),1);
        FattoreDiForma=zeros(size(ListaNodiCandidati,1),1);
        for i=1:size(ListaNodiCandidati,1)
            FacceAggiunte=[SuperficieGenerazione(FacciaScelta,1) SuperficieGenerazione(FacciaScelta,2) ListaNodiCandidati(i); SuperficieGenerazione(FacciaScelta,1) SuperficieGenerazione(FacciaScelta,3) ListaNodiCandidati(i); SuperficieGenerazione(FacciaScelta,2) SuperficieGenerazione(FacciaScelta,3) ListaNodiCandidati(i)];
            FacceTot=sort(SuperficieGenerazione,2);
            %FacceTot=union(sort(SuperficieGenerazione,2),sort([t(:,[1,2,3]);t(:,[1,2,4]); t(:,[1,3,4]); t(:,[2,3,4])],2),'rows'); %usare questa riga solo per sicurezza... in realtà basta fare il controllo solo sulle facce di superficie
            FacceAggiunte=FacceAggiunte(find(~ismember(sort(FacceAggiunte,2),FacceTot,'rows')),:);
            NFA=size(FacceAggiunte,1);
            NFT=size(FacceTot,1);
            Intersezione=zeros(NFA,NFT);
            if (NFA>0)
                N1=cross(Griglia(FacceAggiunte(:,2),:)-Griglia(FacceAggiunte(:,1),:),Griglia(FacceAggiunte(:,3),:)-Griglia(FacceAggiunte(:,1),:),2);
                N2=cross(Griglia(FacceTot(:,2),:)-Griglia(FacceTot(:,1),:),Griglia(FacceTot(:,3),:)-Griglia(FacceTot(:,1),:),2);
                d1=-dot(N1,Griglia(FacceAggiunte(:,1),:),2);
                d2=-dot(N2,Griglia(FacceTot(:,1),:),2);
                dV1_1=reshape(dot( repmat(N2,[1 1 NFA]), repmat(reshape(Griglia(FacceAggiunte(:,1),:)',[1 3 NFA]),[NFT,1,1]),2) + repmat(d2,[1 1 NFA]),[NFT NFA])'; 
                dV1_2=reshape(dot( repmat(N2,[1 1 NFA]), repmat(reshape(Griglia(FacceAggiunte(:,2),:)',[1 3 NFA]),[NFT,1,1]),2) + repmat(d2,[1 1 NFA]),[NFT NFA])';
                dV1_3=reshape(dot( repmat(N2,[1 1 NFA]), repmat(reshape(Griglia(FacceAggiunte(:,3),:)',[1 3 NFA]),[NFT,1,1]),2) + repmat(d2,[1 1 NFA]),[NFT NFA])';
                dV2_1=reshape(dot( repmat(N1,[1 1 NFT]), repmat(reshape(Griglia(FacceTot(:,1),:)',[1 3 NFT]),[NFA,1,1]),2) + repmat(d1,[1 1 NFT]),[NFA NFT]);
                dV2_2=reshape(dot( repmat(N1,[1 1 NFT]), repmat(reshape(Griglia(FacceTot(:,2),:)',[1 3 NFT]),[NFA,1,1]),2) + repmat(d1,[1 1 NFT]),[NFA NFT]);
                dV2_3=reshape(dot( repmat(N1,[1 1 NFT]), repmat(reshape(Griglia(FacceTot(:,3),:)',[1 3 NFT]),[NFA,1,1]),2) + repmat(d1,[1 1 NFT]),[NFA NFT]);
                CoppieDisgiunte=find((dV1_1.*dV1_2>0 & dV1_2.*dV1_3>0) | (dV2_1.*dV2_2>0 & dV2_2.*dV2_3>0));
                CoppieComplanari=find( dV1_1==0 & dV1_2==0 & dV1_3==0);
                CoppieIntrippate=setdiff(setdiff(find(ones(NFA,NFT)>0),CoppieComplanari),CoppieDisgiunte);
                %Esaminiamo le coppie di triangoli intrippate (si veda articolo sull'intersezione di triangoli)
                if (length(CoppieIntrippate)>0)
                    D=cross(repmat(N1,[1 1 NFT]),repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]),2);
                    det=reshape(dot(repmat(N1,[1 1 NFT]),repmat(N1,[1 1 NFT]),2).*dot(repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]),repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]),2)-dot(repmat(N1,[1 1 NFT]),repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]),2).^2,[NFA NFT]);
                    warning off;
                    c1=(-repmat(d1,1,NFT).*reshape(dot(repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]),repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]),2),[NFA NFT])+repmat(reshape(d2,[1 NFT]),NFA,1).*reshape(dot(repmat(N1,[1 1 NFT]),repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]),2),[NFA NFT]))./det;
                    c2=(-repmat( reshape(d2,[1 NFT]),NFA,1).*reshape(dot(repmat(N1,[1 1 NFT]),repmat(N1,[1 1 NFT]),2),[NFA NFT])+repmat(d1,1,NFT).*reshape(dot(repmat(N1,[1 1 NFT]),repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]),2),[NFA NFT]))./det;
                    O=repmat(reshape(c1,[NFA,1,NFT]),[1 3 1]).*repmat(N1,[1 1 NFT])+repmat(reshape(c2,[NFA,1,NFT]),[1 3 1]).*repmat(reshape(N2',[1 3 NFT]),[NFA 1 1]);
                    D=D./repmat(sqrt(sum(D.^2,2)),[1 3 1]);
                    warning on;
                    %calcolo l'intersezione dei due triangoli con la retta O+tD
                    pV1_1=reshape(dot(D,repmat(Griglia(FacceAggiunte(:,1),:),[1,1,NFT])-O,2),[NFA, NFT]);
                    pV1_2=reshape(dot(D,repmat(Griglia(FacceAggiunte(:,2),:),[1,1,NFT])-O,2),[NFA, NFT]);
                    pV1_3=reshape(dot(D,repmat(Griglia(FacceAggiunte(:,3),:),[1,1,NFT])-O,2),[NFA, NFT]);
                    I12=intersect(find(dV1_1.*dV1_2>0 | (dV1_1==0 & dV1_2==0) ),CoppieIntrippate);
                    I13=setdiff(intersect(find(dV1_1.*dV1_3>0  | (dV1_1==0 & dV1_3==0)),CoppieIntrippate),I12);
                    I23=setdiff(intersect(find(dV1_2.*dV1_3>0 | (dV1_2==0 & dV1_3==0)),CoppieIntrippate),union(I12,I13));
                    I12=union(I12, setdiff(intersect(find(dV1_1.*dV1_2>=0),CoppieIntrippate),union(I12,union(I13,I23))));
                    I13=union(I13, setdiff(intersect(find(dV1_1.*dV1_3>=0),CoppieIntrippate),union(I12,union(I13,I23))));
                    I23=union(I23, setdiff(intersect(find(dV1_2.*dV1_3>=0),CoppieIntrippate),union(I12,union(I13,I23))));                    
                    t1_1=zeros(NFA, NFT);
                    t1_2=zeros(NFA, NFT);
                    t1_1(I12)=pV1_2(I12)+(pV1_3(I12)-pV1_2(I12)).*dV1_2(I12)./(dV1_2(I12)-dV1_3(I12));
                    t1_1(I13)=pV1_1(I13)+(pV1_2(I13)-pV1_1(I13)).*dV1_1(I13)./(dV1_1(I13)-dV1_2(I13));
                    t1_1(I23)=pV1_3(I23)+(pV1_1(I23)-pV1_3(I23)).*dV1_3(I23)./(dV1_3(I23)-dV1_1(I23));
                    t1_2(I12)=pV1_1(I12)+(pV1_3(I12)-pV1_1(I12)).*dV1_1(I12)./(dV1_1(I12)-dV1_3(I12));
                    t1_2(I13)=pV1_3(I13)+(pV1_2(I13)-pV1_3(I13)).*dV1_3(I13)./(dV1_3(I13)-dV1_2(I13));
                    t1_2(I23)=pV1_2(I23)+(pV1_1(I23)-pV1_2(I23)).*dV1_2(I23)./(dV1_2(I23)-dV1_1(I23));
                    t1=t1_1;
                    t1(:,:,2)=t1_2;
                    t1=sort(t1,3);
                    pV2_1=reshape(dot(D,repmat(reshape(Griglia(FacceTot(:,1),:)',[1 3 NFT]),[NFA,1,1])-O,2),[NFA, NFT]);
                    pV2_2=reshape(dot(D,repmat(reshape(Griglia(FacceTot(:,2),:)',[1 3 NFT]),[NFA,1,1])-O,2),[NFA, NFT]);
                    pV2_3=reshape(dot(D,repmat(reshape(Griglia(FacceTot(:,3),:)',[1 3 NFT]),[NFA,1,1])-O,2),[NFA, NFT]);
                    I12=intersect(find(dV2_1.*dV2_2>0 | (dV2_1==0 & dV2_2==0) ),CoppieIntrippate);
                    I13=setdiff(intersect(find(dV2_1.*dV2_3>0  | (dV2_1==0 & dV2_3==0)),CoppieIntrippate),I12);
                    I23=setdiff(setdiff(intersect(find(dV2_2.*dV2_3>0 | (dV2_2==0 & dV2_3==0)),CoppieIntrippate),I12),I13);
                    I12=union(I12, setdiff(intersect(find(dV2_1.*dV2_2>=0),CoppieIntrippate),union(I12,union(I13,I23))));
                    I13=union(I13, setdiff(intersect(find(dV2_1.*dV2_3>=0),CoppieIntrippate),union(I12,union(I13,I23))));
                    I23=union(I23, setdiff(intersect(find(dV2_2.*dV2_3>=0),CoppieIntrippate),union(I12,union(I13,I23))));
                    t2_1=zeros(NFA, NFT);
                    t2_2=zeros(NFA, NFT);
                    t2_1(I12)=pV2_2(I12)+(pV2_3(I12)-pV2_2(I12)).*dV2_2(I12)./(dV2_2(I12)-dV2_3(I12));
                    t2_1(I13)=pV2_1(I13)+(pV2_2(I13)-pV2_1(I13)).*dV2_1(I13)./(dV2_1(I13)-dV2_2(I13));
                    t2_1(I23)=pV2_3(I23)+(pV2_1(I23)-pV2_3(I23)).*dV2_3(I23)./(dV2_3(I23)-dV2_1(I23));
                    t2_2(I12)=pV2_1(I12)+(pV2_3(I12)-pV2_1(I12)).*dV2_1(I12)./(dV2_1(I12)-dV2_3(I12));
                    t2_2(I13)=pV2_3(I13)+(pV2_2(I13)-pV2_3(I13)).*dV2_3(I13)./(dV2_3(I13)-dV2_2(I13));
                    t2_2(I23)=pV2_2(I23)+(pV2_1(I23)-pV2_2(I23)).*dV2_2(I23)./(dV2_2(I23)-dV2_1(I23));
                    t2=t2_1;
                    t2(:,:,2)=t2_2;
                    t2=sort(t2,3);
                    Inf=max(t1(:,:,1),t2(:,:,1));
                    Sup=min(t1(:,:,2),t2(:,:,2));
                    CoppieIntrippateIntersecate=intersect(find(Inf<=Sup),CoppieIntrippate);
                    [I J]=ind2sub([NFA NFT],CoppieIntrippateIntersecate);
                    %Le uniche intersezioni accettate per le coppie di triangoli intrippate sono quelle in cui i 2 triangoli
                    %hanno un lato in comune (e quindi l'intersezione coincide col lato) oppure hanno un vertice in comune e si
                    %intersecano SOLO in quel punto
                    NodiSecondo=[FacceTot(J,1),FacceTot(J,2),FacceTot(J,3)];
                    NumNodiComuni=sum([NodiSecondo-repmat(FacceAggiunte(I,1),1,3) NodiSecondo-repmat(FacceAggiunte(I,2),1,3) NodiSecondo-repmat(FacceAggiunte(I,3),1,3)]==0,2);
                    if (sum(NumNodiComuni==2 | (NumNodiComuni==1 & reshape(Sup(CoppieIntrippateIntersecate),length(CoppieIntrippateIntersecate),1)-reshape(Inf(CoppieIntrippateIntersecate),length(CoppieIntrippateIntersecate),1)<TOLLERANZAPICCOLA ))<size(CoppieIntrippateIntersecate,1))
                        CandidatoBuono(i)=0;
                        continue;
                    end
                end
                %Esaminiamo ora le coppie complanari di triangoli
                if (length(CoppieComplanari)>0)
                    [I J]=ind2sub([NFA NFT],CoppieComplanari);
                    O1=Griglia(FacceAggiunte(I,1),:); O2=Griglia(FacceTot(J,1),:);
                    O1(:,:,2)=Griglia(FacceAggiunte(I,2),:); O2(:,:,2)=Griglia(FacceTot(J,1),:);
                    O1(:,:,3)=Griglia(FacceAggiunte(I,3),:); O2(:,:,3)=Griglia(FacceTot(J,1),:);
                    O1(:,:,4)=Griglia(FacceAggiunte(I,1),:); O2(:,:,4)=Griglia(FacceTot(J,2),:);
                    O1(:,:,5)=Griglia(FacceAggiunte(I,2),:); O2(:,:,5)=Griglia(FacceTot(J,2),:);
                    O1(:,:,6)=Griglia(FacceAggiunte(I,3),:); O2(:,:,6)=Griglia(FacceTot(J,2),:);
                    O1(:,:,7)=Griglia(FacceAggiunte(I,1),:); O2(:,:,7)=Griglia(FacceTot(J,3),:);
                    O1(:,:,8)=Griglia(FacceAggiunte(I,2),:); O2(:,:,8)=Griglia(FacceTot(J,3),:);
                    O1(:,:,9)=Griglia(FacceAggiunte(I,3),:); O2(:,:,9)=Griglia(FacceTot(J,3),:);
                    D1=Griglia(FacceAggiunte(I,2),:)-Griglia(FacceAggiunte(I,1),:); D2=Griglia(FacceTot(J,2),:)-Griglia(FacceTot(J,1),:);
                    D1(:,:,2)=Griglia(FacceAggiunte(I,3),:)-Griglia(FacceAggiunte(I,2),:); D2(:,:,2)=Griglia(FacceTot(J,2),:)-Griglia(FacceTot(J,1),:);
                    D1(:,:,3)=Griglia(FacceAggiunte(I,1),:)-Griglia(FacceAggiunte(I,3),:); D2(:,:,3)=Griglia(FacceTot(J,2),:)-Griglia(FacceTot(J,1),:);
                    D1(:,:,4)=Griglia(FacceAggiunte(I,2),:)-Griglia(FacceAggiunte(I,1),:); D2(:,:,4)=Griglia(FacceTot(J,3),:)-Griglia(FacceTot(J,2),:);
                    D1(:,:,5)=Griglia(FacceAggiunte(I,3),:)-Griglia(FacceAggiunte(I,2),:); D2(:,:,5)=Griglia(FacceTot(J,3),:)-Griglia(FacceTot(J,2),:);
                    D1(:,:,6)=Griglia(FacceAggiunte(I,1),:)-Griglia(FacceAggiunte(I,3),:); D2(:,:,6)=Griglia(FacceTot(J,3),:)-Griglia(FacceTot(J,2),:);
                    D1(:,:,7)=Griglia(FacceAggiunte(I,2),:)-Griglia(FacceAggiunte(I,1),:); D2(:,:,7)=Griglia(FacceTot(J,1),:)-Griglia(FacceTot(J,3),:);
                    D1(:,:,8)=Griglia(FacceAggiunte(I,3),:)-Griglia(FacceAggiunte(I,2),:); D2(:,:,8)=Griglia(FacceTot(J,1),:)-Griglia(FacceTot(J,3),:);
                    D1(:,:,9)=Griglia(FacceAggiunte(I,1),:)-Griglia(FacceAggiunte(I,3),:); D2(:,:,9)=Griglia(FacceTot(J,1),:)-Griglia(FacceTot(J,3),:);
                    D2_bis=Griglia(FacceTot(J,3),:)-Griglia(FacceTot(J,1),:);
                    D2_bis(:,:,2)=Griglia(FacceTot(J,3),:)-Griglia(FacceTot(J,1),:);
                    D2_bis(:,:,3)=Griglia(FacceTot(J,3),:)-Griglia(FacceTot(J,1),:);
                    D2_bis(:,:,4)=Griglia(FacceTot(J,1),:)-Griglia(FacceTot(J,2),:);
                    D2_bis(:,:,5)=Griglia(FacceTot(J,1),:)-Griglia(FacceTot(J,2),:);
                    D2_bis(:,:,6)=Griglia(FacceTot(J,1),:)-Griglia(FacceTot(J,2),:);
                    D2_bis(:,:,7)=Griglia(FacceTot(J,2),:)-Griglia(FacceTot(J,3),:);
                    D2_bis(:,:,8)=Griglia(FacceTot(J,2),:)-Griglia(FacceTot(J,3),:);
                    D2_bis(:,:,9)=Griglia(FacceTot(J,2),:)-Griglia(FacceTot(J,3),:);
                    %Controlliamo che la nuova faccia non sia interamente
                    %contenuta in una faccia già esistente
                    App=reshape(dot(cross(D2,O1-O2,2),cross(D2,D2_bis,2),2)>0,[size(CoppieComplanari,1),9]);
                    if (size(find(sum(App(:,[1 4 7]),2)==3 | sum(App(:,[2 5 8]),2)==3 | sum(App(:,[3 6 9]),2)==3 ),1)>0)
                        CandidatoBuono(i)=0;
                        continue;
                    end
                    Det=dot(D1,D1,2).*dot(D2,D2,2)-dot(D1,D2,2).^2;
                    Det=reshape(Det,[size(CoppieComplanari,1),9]);
                    LatiIntersecosi=zeros(size(CoppieComplanari,1),9);
                    LatiParalleli=find(Det==0);
                    warning off;
                    t1=reshape(-dot(D2,D2,2).*dot(D1,O1-O2,2)+dot(D1,D2,2).*dot(D2,O1-O2,2),[size(CoppieComplanari,1),9])./Det;
                    t2=reshape(-dot(D1,D2,2).*dot(D1,O1-O2,2)+dot(D1,D1,2).*dot(D2,O1-O2,2),[size(CoppieComplanari,1),9])./Det;
                    warning on;
                    LatiIntersecosi(find( t1>=0 & t1<=1 & t2>=0 & t2<=1))=1;
                    LatiIntersecosi(LatiParalleli)=0;
                    DaControllare=find(LatiIntersecosi);
                    if (length(find((t1(DaControllare)==0 & t2(DaControllare)==0) | (t1(DaControllare)==1 & t2(DaControllare)==1) | (t1(DaControllare)==1 & t2(DaControllare)==0) | (t1(DaControllare)==0 & t2(DaControllare)==1)))<size(DaControllare,1))
                        CandidatoBuono(i)=0;
                        continue;
                    end
                    LatiIsoRettici=intersect(find( reshape(sum(cross(D1,O1-O2,2).^2,2),[size(CoppieComplanari,1),9])==0 ),LatiParalleli);
                    t_=reshape(dot(D1,O2-O1,2)./sum(D1.^2,2),[size(CoppieComplanari,1),9]);
                    t_(:,:,2)=reshape(dot(D1,O2+D2-O1,2)./sum(D1.^2,2),[size(CoppieComplanari,1),9]);
                    t_=sort(t_,3);
                    App=max(t_(:,:,1),0)<=min(t_(:,:,2),1);
                    LatiIntersecosi(LatiIsoRettici)=App(LatiIsoRettici);
                    LatiIsoRetticiIntersecanti=intersect(find(App>0),LatiIsoRettici);
                    t1=t_(:,:,1);
                    t2=t_(:,:,2);
                    if (length(find(   (t1(LatiIsoRetticiIntersecanti)>1-TOLLERANZAPICCOLA & t1(LatiIsoRetticiIntersecanti)<1+TOLLERANZAPICCOLA) | abs(t2(LatiIsoRetticiIntersecanti))<TOLLERANZAPICCOLA | ( abs(t1(LatiIsoRetticiIntersecanti))<TOLLERANZAPICCOLA & t2(LatiIsoRetticiIntersecanti)>1-TOLLERANZAPICCOLA & t2(LatiIsoRetticiIntersecanti)<1+TOLLERANZAPICCOLA )))<size(LatiIsoRetticiIntersecanti,1))
                        CandidatoBuono(i)=0;
                        continue;
                    end
                end
                CandidatoBuono(i)=1;
                TrovatoCandidato=1;
            else
                CandidatoBuono(i)=1;
                TrovatoCandidato=1;
            end
            %Calcolo del fattore di forma del tetraedro
            %Calcolo del circoraggio
            MatriceApp=[Griglia(SuperficieGenerazione(FacciaScelta,2),:)-Griglia(SuperficieGenerazione(FacciaScelta,1),:); Griglia(SuperficieGenerazione(FacciaScelta,3),:)-Griglia(SuperficieGenerazione(FacciaScelta,1),:); Griglia(ListaNodiCandidati(i),:)-Griglia(SuperficieGenerazione(FacciaScelta,1),:)];
            CircoRaggio=sqrt(sum((inv(2*MatriceApp)*sum(MatriceApp.^2,2)).^2));
            %Il fattore di forma è il rapporto tra il volume del tetraedro
            %ed il volume del tetraedro equilatero con lo stesso cicoraggio (normalizzato a [0 1])
            V_opt=(8*CircoRaggio^3)/(9*sqrt(3));
            FattoreDiForma(i)=(V_opt-abs(dot(MatriceApp(1,:),cross(MatriceApp(2,:),MatriceApp(3,:))))/6)/V_opt;
        end
        
        if (TrovatoCandidato)
            IndiciBuoni=find(CandidatoBuono);
            [App Perm]=sort(FattoreDiForma(IndiciBuoni));
            IndiciBuoni=IndiciBuoni(Perm);
            L_eletto=ListaNodiCandidati(IndiciBuoni(1));
            QualitaTet=[QualitaTet; App(1)];
            %creo il tetraedro costituito dai tre nodi della faccia scelta
            %e dal nodo eletto ed aggiorno tutte le strutture dati
            NuovoTetraedro=[SuperficieGenerazione(FacciaScelta,1) SuperficieGenerazione(FacciaScelta,2) SuperficieGenerazione(FacciaScelta,3) L_eletto];
            FacceTetraedro=sort([SuperficieGenerazione(FacciaScelta,1) SuperficieGenerazione(FacciaScelta,2) L_eletto; SuperficieGenerazione(FacciaScelta,1) SuperficieGenerazione(FacciaScelta,3) L_eletto; SuperficieGenerazione(FacciaScelta,2) SuperficieGenerazione(FacciaScelta,3) L_eletto; SuperficieGenerazione(FacciaScelta,1) SuperficieGenerazione(FacciaScelta,2) SuperficieGenerazione(FacciaScelta,3)],2);
            VerticiOpposti=[setdiff(NuovoTetraedro,FacceTetraedro(1,:)); setdiff(NuovoTetraedro,FacceTetraedro(2,:)); setdiff(NuovoTetraedro,FacceTetraedro(3,:)); setdiff(NuovoTetraedro,FacceTetraedro(4,:))];
            IndiciDaAggiungere=find(~ismember(FacceTetraedro,sort(SuperficieGenerazione,2),'rows'));
            FacceDaAggiungere=FacceTetraedro(IndiciDaAggiungere,:);
            VerticiOpposti=VerticiOpposti(IndiciDaAggiungere,:);
            NomaliDaAggiungere=cross(Griglia(FacceDaAggiungere(:,2),:)-Griglia(FacceDaAggiungere(:,1),:),Griglia(FacceDaAggiungere(:,3),:)-Griglia(FacceDaAggiungere(:,1),:),2);
            DaInvertire=find(dot(Griglia(VerticiOpposti,:)-Griglia(FacceDaAggiungere(:,1),:),NomaliDaAggiungere,2)>0);
            NomaliDaAggiungere(DaInvertire,:)=-NomaliDaAggiungere(DaInvertire,:);
            NomaliDaAggiungere=NomaliDaAggiungere./repmat(sqrt(sum(NomaliDaAggiungere.^2,2)),1,3);
            IndiciDaEliminare=find(ismember(sort(SuperficieGenerazione,2),FacceTetraedro,'rows'));
            SuperficieGenerazione=SuperficieGenerazione(setdiff(1:size(SuperficieGenerazione,1),IndiciDaEliminare),:);
            NormaliSupGenerazione=NormaliSupGenerazione(setdiff(1:size(NormaliSupGenerazione,1),IndiciDaEliminare),:);
            t=[t; NuovoTetraedro];
            if (ismember(L_eletto,ListaNodiInterni)) 
                % Se il nuovo tetraedro è stato costruito su un nodo nuovo
                % allora i 3 vertici di base restano ancora vertici di superficie
                hProgBar = progressbar( hProgBar,1/(CostanteProgBar) );
                ListaNodiInterni=setdiff(ListaNodiInterni,L_eletto);
                ListaNodiSuperficie=[ListaNodiSuperficie; L_eletto];
            else
                % In caso contrario i nodi interni non variano... devo solo togliere dai nodi di superficie quelli
                % che appartengono a 3 facce del tetraedro che non sono state aggiunte
                if (size(FacceDaAggiungere,1)==1) 
                    %In questo caso tolgo solo il nodo che non appartiene
                    %alla faccia aggiunta solo se non serve a descrivere altre facce superficiali
                    if (~ismember(VerticiOpposti,SuperficieGenerazione))
                        ListaNodiSuperficie=setdiff(ListaNodiSuperficie,VerticiOpposti);
                    end
                end
                if (size(FacceDaAggiungere,1)==0)
                    %In questo caso vanno tolti tutti sempre se non servono più
                    App=find(~ismember(NuovoTetraedro,SuperficieGenerazione));
                    if (length(App)>0)
                        ListaNodiSuperficie=setdiff(ListaNodiSuperficie,NuovoTetraedro(App)');
                    end
                end
            end
            SuperficieGenerazione=[SuperficieGenerazione; FacceDaAggiungere];
            NormaliSupGenerazione=[NormaliSupGenerazione; NomaliDaAggiungere];
            if (etime(clock,Tempo)>TempoVis)
                Tempo=clock;
                hold off;
                incl=find(Griglia(:,1)>SezioneX);
                tri1=SuperficieGenerazione(any(ismember(SuperficieGenerazione,incl),2),:);
                h=trimesh(tri1,Griglia(:,1),Griglia(:,2),Griglia(:,3));
                hold on;
                set(h,'edgecolor',[0.5058 0.7882 0.0352]);
                axis equal;
                drawnow
            end
        else
           % fprintf('\nCacchio non l''ho trovato!!!\n');
        end
    end
end
progressbar(hProgBar,-1 );
fprintf('fatto\n');
p=Griglia;
h=figure(100);
hist(QualitaTet,50);
pause;
close(h);




