function [ia,ja,nterm]=topol(nnodi,ntria,n1,triang)

%*******************************TOPOL*************************************
%
% Routine per l'analisi topologica del reticolo nel caso simmetrico
% 
%
% 
% Dati di input:
% 1) n1: massimo numero di contatti nodali
% 2) nnodi: numero nodi della griglia
% 4) ntria: numero triangoli
% 3) triang(4,ntria): tabella di incidenza (nodi dei triangoli)

% Dati di output:
% 1) JA: vettore degli indici di colonna
% 2) IA: vettore topologico
% 3) nterm: numero di elementi contenuti in JA
%    (= n. elementi non nulli della matrice)
% 
%


%
% fissa gli elementi diagonali a distanza costante n1
%
%

nterm=n1*nnodi;
for i=1:nterm
    JA(i)=0;
end
s=1;
for k=1:n1:nterm-n1+1
    JA(k)=s;
    s=s+1;
end
%
% analizza tutti i triangoli (ntria)
%
% TRIANG(4,ntria): matrice contenente la tabella della successione di nodi
%               per ciascun elemento
%
TRIANG=triang;
TRIANG=TRIANG';
imax=10^6;
for k=1:ntria
    for i=1:3
        I2(i)=TRIANG(i,k);
    end
%
% ordina i nodi dell'elemento in senso crescente in I1(3)
% 
% imax: massimo intero rappresentabile dalla macchina (definito 
%       dall'utente)
%
    for j=1:3
        kkk=1;
        kk=I2(1);
        for i=2:3
            if (I2(i)<kk)
                kk=I2(i);
                kkk=i;
            end
        end
        I1(j)=kk;
        I2(kkk)=imax;
    end
% 
% alla fine del ciclo in I1 abbiamo la terna ordinata mentre in I2 i valori
% sono tutti pari a imax
%
% genera il vettore JA
% 
    for i=1:2
        j=i+1;
        for q=j:3
            m=n1*(I1(i)-1)+q-j+2;
            mcontr=n1*I1(i);
            continua='false';
            while ((((I1(q)-JA(m))<0) | (I1(q)-JA(m))>0) & continua=='false')
                if (I1(q)-JA(m))>0
                    if JA(m)==0
                        JA(m)=I1(q);
                        continua='true';
                    else
                        m=m+1;
                        if (m-mcontr)>=0
                            iii=mcontr/n1;
                            fprintf(' RIGA=%6.0f TRIANGOLO=%6.0f AUMENTARE n1 ', iii,k);
                            break
                        else
                            continua='false';
                        end
                    end
                elseif (I1(q)-JA(m))<0
                    mm=m;
                    mm=mm+1;
                    if (m-mcontr)>=0
                        iii=mcontr/n1;
                        fprintf(' RIGA=%6.0f TRIANGOLO=%6.0f AUMENTARE n1 ', iii,k);
                        break
                    else 
                        continua2='false';
                        while ((JA(mm)>0 | JA(mm)<0) & continua2=='false')
                            mm=mm+1;
                            if (m-mcontr)>=0
                                iii=mcontr/n1;
                                fprintf(' RIGA=%6.0f TRIANGOLO=%6.0f AUMENTARE n1 ', iii,k);
                                break
                            else
                                continua2='false';
                            end
                        end
                        while JA(mm)==0
                            JA(mm)=JA(mm-1);
                            mm=mm-1;
                            while mm>m
                                JA(mm)=JA(mm-1);
                                mm=mm-1;
                            end
                        end
                    end
                    JA(m)=I1(q);
                end
            end
        end
    end
end
%
% costruisce il vettore IA
%
IA(1)=1;
m=1;
j=1;
for k=1:n1:nterm
    for i=1:n1
        if JA(k+i-1)~=0
            m=m+1;
        end
    end
    j=j+1;
    IA(j)=m;
end
%
% compatta il vettore JA eliminando gli zeri
% 
m=1;
for k=1:nterm
    if JA(k)~=0
        JA1(m)=JA(k);
        m=m+1;   
    end
end
JA=JA1;
nterm=m-1;

