% DATI DI INPUT
%  F1 =  base del dominio
%  F2 =  altezza del dominio
%  X_0 = ascissa del punto che consideriamo origine del dominio
%  Y_0 = ascissa del punto che consideriamo origine del dominio
%
%  ___________________________________________
% |                                           |
% |                                           |
% |                                           |
% |                                           |
% |                                           |
% |                                           |  F2
% |                                           |
% |                                           |
% |                                           |
% |                                           |
% |                                           |
% |                                           |
%  ____________________________________________
%                      F1
%  DATI DI OUTPUT utili
%  la matrice TRIANG che ha le coordinate dei nodi della triangolazione
%  XY(:,1) ascisse dei nodi
%  XY(:,2) ordinate dei nodi
%  NN numero dei nodi della triangolazione
%  NT numero dei triangoli
clear
F1=input(' lato F1 ');
F2=input(' lato F2 ');
M1=input(' suddivisioni lungo x ');
M2=input(' suddivisioni lungo y ');
X_0=input(' punto X_0 ');
Y_0=input(' punto Y_0 ');
N1=M1+1;
N2=M2+1;
NN=N1*N2;
NT=2*M1*M2;
icont=0;
for i=1:M2;
    E(3)=i;
    for j=1:M1;  
        E(1)=E(3);
        E(2)=E(1)+N2+1;
        E(3)=E(2)-1;
	icont=icont+1;
	TRIANG(icont,1)=E(1);
	TRIANG(icont,2)=E(2);
	TRIANG(icont,3)=E(3);
     end
     A(3)=i+1;
     for j=1:M1;
         A(2)=A(3);
         A(1)=A(3)-1;
         A(3)=A(2)+N2;
	icont=icont+1;
	TRIANG(icont,1)=A(1);
	TRIANG(icont,2)=A(2);
	TRIANG(icont,3)=A(3);
      end
end
% CALCOLO DELLE COORDINATE DEI NODI
Z1=F1/M1;
Z2=F2/M2;
IA1=0;
for   j=1:N1;
     for  i=1:N2;
          IA1=IA1+1;
          S1=i;
          S2=j;
          XY(IA1,1)= X_0 + (S2-1.0)*Z1;
          if abs(XY(IA1,1))< 1.e-13 
                XY(IA1,1)=0.0;
          end
          XY(IA1,2)= Y_0+ F2-(S1-1.0)*Z2;
          if abs(XY(IA1,2)) < 1.e-13 
              XY(IA1,2)=0.0D+00;
          end
     end
end
clear S1 S2 IA1 i j Z1 Z2 E A N1 N2
save -ascii matrice_triang.dat TRIANG 
%  in matrice_triang.dat salviamo la matrice TRIANG come matrice
%  di dimensione NT x 3
save -ascii xycoord.dat XY 
%  in xycoord.dat salviamo la matrice XY delle ascisse e ordinate
%  dei nodi
%
%  VISUALIZZAZIONE
%  la function del MATLAB triplot permette di visualizzare la griglia
%  il comando text( ...) permette di numerare i nodi.
triplot(TRIANG,XY(:,1),XY(:,2))
if (F1==F2) 
	axis equal
end
for i=1:NN
  text(XY(i,1),XY(i,2), int2str(i), 'Color', 'blue', 'FontSize', 20)
end
% Se si vuole visualizzare la numerazione dei triangoli, conviene
% calcolare il baricentro di ciascun triangolo, e su questi
% andare a poi a scrivere il valore corrispondente nella numerazione
% Nel caso in cui sembra piu' conveniente fare due grafici, uno
% solo con la numerazione dei nodi e uno solo con la numerazione
% dei triangoli, basta modificare leggermente il codice, togliendo
% il commento alle seguenti righe
% figure(2)
% triplot(TRIANG,XY(:,1),XY(:,2))
for i=1:NT
	xbar(i) =( XY(TRIANG(i,1),1) + XY(TRIANG(i,2),1) + XY(TRIANG(i,3),1) )/3.;
	ybar(i) =( XY(TRIANG(i,1),2) + XY(TRIANG(i,2),2) + XY(TRIANG(i,3),2) )/3.;
end
for i=1:NT
	text(xbar(i), ybar(i), int2str(i),'FontAngle', 'oblique','Color', 'red')
end
