// ====================================================================
// Fonctions utiles pour Delos
// OS. 14.12.2017 -  version 1.03
//     en cours de modification.....
// Modifications : 
//     10.12.2011, OS : ajout DSWriteMesh, extension de DSShowMesh
//     15.11.2012, OS : ajout DSShowGrd ...
//     14.12.2017, OS : scilab v 6.0.0 est plus contraignant sur la syntaxe []+1 => ERROR
//
// Fonctions pour la geometrie : affichage, lecture/ecriture fichier 
// -----------------------------
// OBJET DSShowGeom : affichage d'une geometrie decrite par une liste d'aretes
//       DSShowGeom(XYcoord,laretes,reg,code,numfig,titre,colorAllReg);
// OBJET DSReadGeom : lecture d'une geometrie decrite dans le fichier "nomf"
//       [XYcoord,laretes,reg,code]=DSReadGeom(nomf);
// OBJET DSWriteGeom : ecrit la geometrie dans un fichier(au format Delos)
//       DSWriteGeom(nomf, XYcoord, laretes, reg, code )
//
// Fonctions pour le maillage : affichage, lecture/ecriture fichier
// ----------------------------
// OBJET DSShowMesh : visualise une triangulation 
//       DSShowMesh(XYcoord,TRI,REG,numfig,titre,colorAllReg)
// OBJET DSReadMesh : lecture d'un maillage dans le fichier "nomfich" 
//       (au format Delos)
//       [XYcoord,TRI,REG]=DSReadMesh(nomfich)
//       [err]=DSWriteMesh(nomfich,XY,TRI,REG)
//
// OBJET AddDENSITY : Ajout de la densite "dens" aux densites deja definies
//       (pour Delos)
//       [XYGeo,SUI,TDEN]=AddDENSITY(dens,XYGeo,SUI,TDEN);
// OBJET DSWriteDens : ecriture des densites dans le fichiers "nomf"
//       [err]=DSWriteDens(nomf,XYGeo,SUI,TDEN)
//
// TODO :
//    Procedure de verification qui fait appel a :
//     [XYM,listiM,prec]=MultiplePointsNorm(XYcoord,dig)
//
// * modifier ColorAllreg dans DSShowMesh(...REG,ColorAllreg) idem DSShowGrd
//   pour pouvoir associer des grandeurs aux triangles (qualite, taille...)
// 
// Pour utiliser ces scripts copiez : exec("DSutil.sci",-1)
// ====================================================================

function [TRIIN]=TriOnNode(ITRNOE,ITRTRI,NOETRI,noeud)
//=============================================================
// Pour l'exercice du TP 4 sur Delaunay
i=0; FIN=%F;TRIIN=[];
//itri=find(NOETRI==noeud,1);
itri=NOETRI(noeud);
if( itri < 1 ) then, printf("Pas de triangles au sommet %d\n",noeud);return;end;
itn=find(ITRNOE(itri,:)==noeud);
if( isempty(itn)) then, 
    printf("Structure non coherente au sommet %d\n",noeud);return;
end;
i=i+1;TRIIN(i)=itri;
while ( ~FIN )
    so=ITRNOE(itri,itn);
    itn=modulo(itn,3)+1;     // le sommet suivant de itn dans itri
    soms=ITRNOE(itri,itn);
    printf("Le sommet suivant de %d dans triangle %d = sommet %d\n",so,itri,soms)
    itrip=itri;
    itri=ITRTRI(itri,itn);  // le triangle oppose
    printf("Le triangle oppose au sommet %d pour le triangle %d = triangle %d\n",soms,itrip,itri);
    if ((itri == TRIIN(1)) | ( itri < 1 )) then
         FIN=%T; 
    else  
        i=i+1;TRIIN(i)=itri;  
        itn=find(ITRNOE(itri,:)==noeud);
    end
end
endfunction

// ====================================================================
// ================== POUR COMPRENDRE DELAUNAY  =======================
// ====================================================================
function [P1,CC]=Cercle3P2PD(P)
//==============================
// CC : Diametre d'un cercle defini par 3 points
P1=P(1,:);
P12=P(2,:)-P(1,:);
P13=P(3,:)-P(1,:);
A=[P12;P13];
B=[P12*P12';P13*P13'];
CC=A\B; // A*CC=B 
CC=CC';
endfunction

function PlotCercle(P1,CC)
  //==========================
  D=sqrt(CC*CC'); //prscal
  C=P1+(CC)/2; // centre
  xarc(C(:,1)-D/2,C(:,2)+D/2,D,D,0,360*64);
endfunction

function [PIC]=PinCercle(P1,CC,P)
// =======================================
// Point dans cercle ?
P1P=P-P1;
PIC= (P1P*P1P') < (CC*P1P');
endfunction
 
// ====================================================================
// ======================== SUR LA GEOMETRIE  =========================
// ====================================================================

function DSShowGeom(XYcoord,laretes,reg,code,numfig,titre,colorAllReg)
// ====================================================================
// OBJET DSShowGeom : affichage d'une geometrie decrite par une liste d'aretes
// XYcoord : coordonnees des points
// laretes : liste des aretes (chaque arete est definie par une liste de sommets (numeros))
// reg(i) : numero(s) de la (des) region(s) a l'aretes i
// regColor(i) = couleur des aretes des regions
txtleg=[];
losangeN=-4; croixC=-3; rouge=5; colorReg=rouge;
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
//printf("nbvar = %d\n",nbvarin)
if ( nbvarin < 7 ) , colorAllReg= []; end
if ( nbvarin < 6 ) , titre= []; end
if ( nbvarin < 5 ) , numfig= []; end
if ( nbvarin < 4 ) , code= []; end
if ( nbvarin < 3 ) , reg= []; end
if ( nbvarin < 2 ) , laretes= []; end
nbreg=0;
if ( ~isempty(reg) ) then
    numreg=gsort(unique(reg(find(reg>0))))
    nbreg=length(numreg)
    if ( isempty(colorAllReg) ) ,  colorAllReg=2:nbreg+1; end
end
if ( isempty(titre)) then, titre="DSShowGeom";  end
if ( isempty(numfig)) then, numfig=max(winsid())+1;  end
fig=scf(numfig); // creation d'une figure
//
XmM=[min(XYcoord(:,1)),max(XYcoord(:,1))];
YmM=[min(XYcoord(:,2)),max(XYcoord(:,2))];
rect=[XmM;YmM];
[rectm]=DSUTWinMargin(rect,0.1); // marge de 10%
plot2d(XYcoord(:,1),XYcoord(:,2),style=losangeN,rect=rectm,frameflag=3)
txtleg="points"
xtitle(titre,"X","Y")
xgrid()
nbcoord=length(XYcoord)/2;
// xstring(XYcoord(:,1),XYcoord(:,2),string([1:nbcoord]));
// remplace par (pour compatibilite avec Scicoslab 4.1.1) :
for j=1:nbcoord
    xstring(XYcoord(j,1),XYcoord(j,2),string(j));
end
// --------------- affichage des aretes ----------------------
if ( ~isempty(laretes)) then
nbare=length(laretes);
//printf("%d \n",nbare)
for i=1:nbare
    areti=laretes(i)';
    if ( isempty(code) | (code==100) ) then
        if ( ~isempty(reg)) then
            colorReg=colorAllReg(find(numreg==max(reg(i,:))));
        end
        plot2d(XYcoord(areti,1),XYcoord(areti,2),style=colorReg)
    else
        plot2d(XYcoord(areti,1),XYcoord(areti,2),style=croixC)
    end
    txtleg=[txtleg,"border_"+string(reg(i,1))+"/"+string(reg(i,2))]
end
// legend(["points","borders"],1); 
legend(txtleg,1); 
xinfo(string(nbreg)+" zones, "+string(nbare)+" borders, "+string(nbcoord)+ " coord")
else
    // seulement des points
legend(["points"],1);    
xinfo(string(nbcoord)+ " coord")
end
//
endfunction

function [XYcoord,laretes,reg,code]=DSReadGeom(nomf)
// ====================================================================
// OBJET DSReadGeom : lecture d'une geometrie decrite dans le fichier "nomf"
laretes=list(); // pas d'arete
code=[]; // pas de type d'arete
reg=[]; // pas de region separee par des aretes
XYcoord=[]; // pas de points

[fd,err]= mopen(nomf,"r");
// ----------------- lecture du bloc des coordonnees --------------------
DEBBLOC="DEBXYZ"; FINBLOC="FINXYZ"; entete=" "
while  entete ~= DEBBLOC , [entete]=mfscanf(fd,"%s");end
[dim,nb]=mfscanf(fd,"%d %d");
XYcoord=zeros(nb,dim);
printf("lecture des %d coordonnees\n",nb)
XYcoord=mfscanf(dim*nb,fd,"%f"); XYcoord=matrix(XYcoord,dim,nb)';
//for i=1:1:nb
//  for j=1:1:dim
//    XYcoord(i,j)=mfscanf(fd,"%f");
//  end
//end
// ----------------- lecture du bloc des coordonnees --------------------
DEBBLOC="DEBARE"; FINBLOC="FINARE"; entete=" "
while  entete ~= DEBBLOC , [entete]=mfscanf(fd,"%s");end
nbare=mfscanf(fd,"%d");
printf("lecture des %d aretes\n",nbare)
for i=1:nbare
    nbs=mfscanf(fd,"%d"); // nombre de sommets
//    ai=size(1,nbs)
if  (nbs==-3) then
    n=mfscanf(3,fd,"%d "); // lecture de la liste des sommets
    ai=n(1):n(3):n(2);
else
    ai=mfscanf(nbs,fd,"%d "); // lecture de la liste des sommets
end
    laretes=lstcat(laretes,list(ai));
    code(i)=mfscanf(fd,"%d"); 
    nbr=mfscanf(fd,"%d"); // nbre de regions
    reg(i,:)=[0,0];
    for j=1:1:nbr
        reg(i,j)=mfscanf(fd,"%d ");
    end
end
// nombre de regions
    numreg=unique(reg(find(reg>0)))
    nbreg=length(numreg)
    if ( exists("VERBOSE") & VERBOSE == %T ) then
        printf("Numero des (%d) regions = \n",nbreg)
        printf("\t%d",numreg)
    end
//
mclose(fd)
endfunction

function [code, reg, laretes]=DSDefaultGeom(XYcoord, laretes, reg, code)
// =====================================================================
// OBJET DSDefaultGeom : complete la definition de la geometrie dans les cas simples
code100=100; defaultReg=1;
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
//printf("nbvar = %d\n",nbvarin)
if ( nbvarin < 4 ) , code= []; end
if ( nbvarin < 3 ) , reg= []; end    
if ( nbvarin < 2 ) , laretes= list(); end    
[nb,dim]=size(XYcoord);
if (( type(laretes) == 10) & ((laretes == "AUTO") | (laretes=="auto")) ) , laretes = [1:nb,1]; end 
if ( type(laretes) == 1) , laretes = list( laretes ) ; end
if ( type(laretes) ~= 15) , printf("ERROR : no edge found\n"); return;  end
nbare=length(laretes);
if ( isempty(code)) , code=zeros(nbare,1)+code100; end
if ( isempty(reg)) , reg=zeros(nbare,2); reg(:,1)=defaultReg; end
endfunction

function DSWriteGeom(nomf, XYcoord, laretes, reg, code )
// ====================================================================
// OBJET DSWriteGeom : ecrit la geometrie dans un fichier(au format Delos)
// nomf : le nom du fichier
// points : coordonnees des points du maillage (bloc DEBXYZ)
// aretes : tableau des lignes de texte correspondant aux arêtes (bloc DEBARE)
//
[code, reg, laretes]=DSDefaultGeom(XYcoord, laretes, reg); // complete la geometrie si necessaire    
[fd,err]= mopen(nomf,"w");
[nb,dim]=size(XYcoord);
mfprintf(fd,"DEBXYZ\n");
mfprintf(fd,"%d %d\n",nb,dim);
for i=1:1:nb
  for j=1:1:dim
    mfprintf(fd,"%f ",XYcoord(i,j));
  end
  mfprintf(fd,"\n");  
end
mfprintf(fd,"FINXYZ\n");
// inutile vvv
if ( isempty(laretes)) then
    nbare=0;
    mfprintf(fd,"DEBARE\n");
    mfprintf(fd,"%d\n",nbare)
    mfprintf(fd,"FINARE\n");
    mclose(fd); 
    return ; 
end
mfprintf(fd,"DEBARE\n");
nbare=length(laretes);
mfprintf(fd,"%d\n",nbare)
for i=1:nbare
    areti=laretes(i);
    mfprintf(fd,"%d ",length(areti))
    //printf("%d ",[length(areti),areti'])
    mfprintf(fd,"%d ",areti')
    mfprintf(fd,"%d 2 %d %d\n",code(i),reg(i,1),reg(i,2))
end
mfprintf(fd,"FINARE\n");
mclose(fd);
endfunction


// ====================================================================
// ======================== SUR LES MAILLAGES =========================
// ====================================================================
function [err]=DSWriteMesh(nomfich,XY,TRI,REG)
// ====================================================================
// OBJET DSWriteMesh : ecriture d'un maillage dans le fichier "nomfich" (au format Delos)
// Lecture du fichier "nomfich" au format delos
// XY(i,2)  : coordonnees du point i
// TRI(i,:) : sommets du triangle i
// REG(i)   : numero de la region du triangle i
fd=mopen(nomfich,'w');
// ---------- elements ---------
[nb,dim]=size(XY);
DEBBLOC="DEBXYZ";ENDBLOC="FINXYZ"
mfprintf(fd,"%s\n",DEBBLOC);
mfprintf(fd,"%d\t%d\n",nb,dim);
for i=1:1:nb
  mfprintf(fd,"%f\t%f\n",XY(i,:));
end
mfprintf(fd,"%s\n",ENDBLOC);
// ---------- elements ---------
[nbe,nbnode]=size(TRI); CODETRI=3 ; 
if ( isempty(REG) ) then
  REG=zeros(nbe,1)+1; // DEFAULTREG=1
end
DEBBLOC="DEBILM";ENDBLOC="FINILM"
mfprintf(fd,"%s\n",DEBBLOC);
mfprintf(fd,"%d\n",nbe);
for i=1:1:nbe
  mfprintf(fd,"%d\t%d\t%d\t%d\t%d\t%d\n",nbnode,TRI(i,:),CODETRI,REG(i));
end
mfprintf(fd,"%s\n",ENDBLOC);
err=mclose(fd);
endfunction

function [XY,TRI,REG,GRD]=DSReadMesh(nomfich)
// ====================================================================
// OBJET DSReadMesh : lecture d'un maillage dans le fichier "nomfich" (au format Delos)
// Lecture du fichier "nomfich" au format delos
// XY(i,2)  : coordonnees du point i
// TRI(i,:) : sommets du triangle i
// REG(i)   : numero de la region du triangle i
// TODO : SANS COMMENTAIRES dans les fichiers de donnees pour l'instant
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
fd=mopen(nomfich,'r');
DEBBLOC="DEBXYZ"
ENDBLOC="FINXYZ"
entete=""
// Lecture du bloc XYZ
while  entete ~= DEBBLOC , [entete]=mfscanf(fd,"%s");end
[dim,nb]=mfscanf(fd,"%d %d");
XY=zeros(nb,dim);
i=1;
for i=1:1:nb
//while ( entete ~= ENDBLOC )
//  entete=mfscanf(fd,"%s");
//  if ( entete ~= ENDBLOC ) then
   [XY(i,:)]=mfscanf(fd," %f %f");
//   [XY(:,i)]=[x;y];
//   [XY(:,i)]=sprintf(" %f %f %f",entete);
//   i=i+1;
//  end
end
// Lecture du bloc ILM
DEBBLOC="DEBILM"
ENDBLOC="FINILM"
entete=""
// Lecture du bloc XYZ
while ( entete ~= DEBBLOC ) 
  entete=mfscanf(fd,"%s");
end
[nbt]=mfscanf(fd,"%d ");
DSTRI=zeros(nbt,6); // nbnode + 3S + region
i=1;
for i=1:1:nbt
   DSTRI(i,:)=mfscanf(fd," %d %d %d %d %d %d");  
end
TRI=DSTRI(:,2:4);
REG=DSTRI(:,6);
mclose(fd);
// ----- lecture des grandeurs ----
GRD=[];
if ( nbvarout == 4 ) then 
  [GRD]=DSReadGrd(nomfich);
  [nbg,dimG]=size(GRD);
  if ( nbg < nb ) then
  printf("ATTENTION : nb grandeurs = %d < %d (nb points)\n",nbg,nb); 
  end
  if ( dimG > 1 ) then
  printf("ATTENTION : grandeurs non-scalaires dim= %d \n",dimG); 
  end
end  
endfunction

function [GRD]=DSReadGrd(nomfich)
// ====================================================================
// Lecture des grandeurs associees a des numeros
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
fd=mopen(nomfich,'r');
DEBBLOC="DEBGRD"
ENDBLOC="FINGRD"
entete=""
// Lecture du bloc GRD --- il peut etre absent meme si donne en sortie
while ( entete ~= DEBBLOC ) 
  entete=mfscanf(fd,"%s");
  err = meof(fd);
  if ( err ~= 0 ) then
    printf("ATTENTION : pas de bloc %s \n",DEBBLOC);
    GRD=[];mclose(fd);
    return    
  end
end
A=mfscanf(fd,"%d %d");
nbg=A(1);dimG=A(2);
GRDL=mfscanf(nbg*(dimG+1),fd,"%f"); 
GRDL=matrix(GRDL,dimG+1,nbg);
GRDL=GRDL';
GRD=GRDL(GRDL(:,1),2:$); // la 1ere colonne = numero des noeuds !
mclose(fd);
endfunction

function [rectm]=DSUTWinMargin(rect,pc)
dxy=rect(:,2)-rect(:,1);
mr=dxy*pc;
mr=[-mr(1),mr(1);-mr(2),mr(2)];
rectm=rect+mr;
endfunction


function DSShowMesh(coord,itrnoe,itreg,numfig,titre,colorAllReg)
// ====================================================================
// OBJET DSShowMesh : visualise une triangulation sous scilab
// coord : coordonnees des points
// itrnoe(t,:) : les 3 sommets constituant le triangle de numero t
// itreg(i): numero de la region a laquelle appartient le triangle i
// coul(m) : couleur associee a la region de numero m
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
//printf("nbvar = %d\n",nbvarin)
if ( nbvarin < 6 ) , colorAllReg= []; end
if ( nbvarin < 5 ) , titre = []; end
if ( nbvarin < 4 ) , numfig = []; end
if ( nbvarin < 3 ) , itreg = []; end
if ( nbvarin < 2 ) , itrnoe = []; end
//
RANDOMCOLOR=%F
[m,n]=size(itrnoe);
nbreg=0; numreg=1;
// A REVOIR : itreg====> la couleur pour afficher chaque triangle (indice dans la palette)
//            colorAllReg =====> la palette !!! matrice ou string !!!!
// ------------ gestion des couleurs ------------------------- appel DSColorTRI
if ( ~isempty(itreg) ) then
  numreg=gsort(unique(itreg(find(itreg>0))))
else
  if ( isempty(colorAllReg) ) then
    numreg=1; itreg=zeros(m,1)+1; // ni region, ni couleur =1
  else
    // pas region mais des couleurs => tirage aleatoire des couleurs :
    RANDOMCOLOR=%T
    numreg=length(colorAllReg); itreg=int(rand(m,1)*(numreg-1)+1);
  end
end
// on remplace le numero de la region par la couleur associee
nbreg=length(numreg)
if ( isempty(colorAllReg) ) ,  colorAllReg=2:nbreg+1; end
itrcol=itreg;
for i=1 : nbreg
    itrcol(find(itreg==numreg(i)))=i; // on remplace par la couleur
  end
// --------------------------- pause;
//
if ( isempty(titre)) then, titre="DSShowMesh";  end
numfig_max = max(winsid()); // Added since v 6.0.0
if ( isempty(numfig_max) ) then numfig_max = 0; end
if ( isempty(numfig)) then, numfig= numfig_max +1;  end                        
//if ( isempty(numfig)) then, numfig=max(winsid())+1;  end                        // ERROR since v 6.0.0
fig=scf(numfig); // creation d'une figure
//
XmM=[min(coord(:,1)),max(coord(:,1))];
YmM=[min(coord(:,2)),max(coord(:,2))];
rect=[XmM;YmM];
//rect=[min(coord(:,1)),max(coord(:,1));min(coord(:,2)),max(coord(:,2))];
[rectm]=DSUTWinMargin(rect,0.1); // marge de 10%
//plot2d(-mr(1),-mr(2),style=-1,rect=rectm,frameflag=3);
plot2d(-rectm(1,1),-rectm(2,1),style=-1,rect=rectm,frameflag=3);
//plot2d(XmM,YmM,style=-1,rect=rectm);
//plot2d(coord(:,1),coord(:,2),style=-3,rect=rectm);
[m,n]=size(itrnoe);
xpols=matrix(coord(itrnoe,1),m,n);
ypols=matrix(coord(itrnoe,2),m,n);
//xset("colormap",rand(m,3)); // RVB pour chaque materiau
//fill=zeros();
//xfpolys(xpols',ypols',[m/4:m/4+m-1]); // pour tracer des polygones
//xfpolys(xpols',ypols',[2:m+1]); // pour tracer des polygones
//couleur=zeros(1,m)+5;
//couleur=coul(itreg);
couleur=colorAllReg(itrcol);
xfpolys(xpols',ypols',couleur); // pour tracer des polygones
xtitle(titre,"X","Y")
xgrid()
// colorbar : a revoir les numeros des regions ne se suivent pas forcement
// if (( nbreg > 1 ) & ( RANDOMCOLOR == %T ) )then
//  colorbar(1,2,itrcol)
// end
//nbnodes = length(coord(itrnoe,2)); // Les aretes ?
nbnodes = length(coord(:,2));
xinfo(string(nbreg)+" zones, "+string(m)+" triangles, "+string(nbnodes)+ " nodes")
endfunction

// ====================================================================
// ======================== SUR LES DENSITES  =========================
// ====================================================================

function [err]=DSWriteDens(nomf,XYGeo,SUI,TDEN)
// ====================================================================
// OBJET DSWriteDens : ecriture des densites dans le fichiers "nomf"
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
//printf("nbvar = %d\n",nbvarin)
if ( nbvarin < 4 ) , printf("ERROR : incomplete definition"); return; end
//
[fd,err]= mopen(nomf,"w");
[nb,dim]=size(XYGeo);
mfprintf(fd,"* generation automatique par DSWriteDens\n"); // ajouter la date...
mfprintf(fd,"DEBGEO\n");
mfprintf(fd,"%d %d\n",nb,dim);
for i=1:1:nb
  mfprintf(fd,"%d ",i);
  for j=1:1:dim
    mfprintf(fd,"%f ",XYGeo(i,j));
  end
  mfprintf(fd,"\n");  
end
mfprintf(fd,"FINGEO\n");
//
nbsui=length(SUI);
mfprintf(fd,"DEBSUI\n");
mfprintf(fd,"%d\n",nbsui);
for isui=1:1:nbsui
  sui=SUI(isui);
  // numero, type, taille, raison
  mfprintf(fd,"%d\t%d\t%f\t%f\n",isui,sui(1),sui(2),sui(3));
end
mfprintf(fd,"FINSUI\n");
//
nbden=length(TDEN);
mfprintf(fd,"DEBDEN\n");
mfprintf(fd,"%d  1\n",nbden); // 1 = generation
for iden=1:1:nbden
  den=TDEN(iden)
  if ( den(2) == 2 ) , mfprintf(fd,"%d\t%d\t%d\t%d\t%d\n",iden,den(1),den(2),den(3),den(4)); end
  if ( den(2) == 3 ) , mfprintf(fd,"%d\t%d\t%d\t%d\t%d\t%d\n",iden,den(1),den(2),den(3),den(4),den(5)); end  
end
mfprintf(fd,"FINDEN\n");
mclose(fd);
err=0; //OK
endfunction


function [XYGeo,SUI,TDEN]=AddDENSITY(dens,XYGeo,SUI,TDEN)
// ====================================================================
// OBJET AddDENSITY : Ajout de la densite "dens" aux densites deja definies (pour Delos)
// conversion d'une DENSITE en objets Delos, cette procedure a pour objectif
// de faciliter la definition des densites : elle complete les definitions...
suite_geometrique=1;
suite_arithmetique=2;
concentration_ponctuelle=1; // 'POIN'
concentration_axiale=2; // 'DROI'
concentration_segment=3;// 'SEGM'
//
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
//printf("nbvar = %d\n",nbvarin)
if ( nbvarin < 4 ) , TDEN= list(); end
if ( nbvarin < 3 ) , SUI= list(); end
if ( nbvarin < 2 ) , XYGeo= []; end
if (( nbvarin < 1 ) | (isempty(dens))) , return; end
    // C'est une initialisation :-)
    //  printf("ERROR : no density given\n"); return; end
if ( ( type(dens) == 1 ) & (length(dens) == 1 ) ) then
    XYcoord=[];    Taille=dens;    RaisonSUI=[];    TypeDEN=[];
  else
    XYcoord=dens(1);    Taille=dens(2);    
    RaisonSUI=dens(3);    TypeDEN=dens(4);
end
[nbp,dim]=size(XYGeo);[nbsui]=length(SUI);[nbden]=length(TDEN)
[nbcoord,dim]=size(XYcoord);[nbt]=length(Taille);[nbr]=length(RaisonSUI);
if ( nbcoord == 0 ) then
  // --------------------- Taille max du maillage --------------------------
  if ( nbt == 1 ) then
    // Taille minimale imposee (en n'importe quel points)
    if ( nbp == 0 ) then
      nbp=nbp+1; XYGeo(nbp,:)=[0,0] ;
  end 
    SUI= lstcat( SUI, list([ suite_geometrique , Taille, 1 ]))
    nbsui=length(SUI);
    TDEN = lstcat( TDEN, list([ concentration_ponctuelle, 2 , nbsui, nbp ]));  // ERROR since v 6.0.0
  else
    printf("Error : no point for %d size of mesh\n",nbt); return;
  end
end
if ( nbcoord == 1 ) then
  // --------------------- Concentration ponctuelle --------------------------
  if ( nbt == 1 ) then
    nbp=nbp+1;    XYGeo(nbp,:)=XYcoord;
    SUI= lstcat( SUI, list([ suite_geometrique , Taille, RaisonSUI]));
    nbsui=length(SUI);
    TDEN = lstcat( TDEN, list([ concentration_ponctuelle, 2 , nbsui, nbp]));     
  else
    printf("Error : one point => one size (%d found)\n",nbt); return;
  end  
end
if ( nbcoord > 1 ) then
  // --------------------- Concentration axiale --------------------------
  nbp=nbp+nbcoord;   XYGeo=[XYGeo;XYcoord];
  if ( nbt == 1 ) then
    if ( TypeDEN == 'POIN' ) then
      // tous les points sont a la meme concentration ponctuelle
      if ( nbr == 1 ) then
        for i=1:nbt
          SUI= lstcat( SUI, list([ suite_geometrique , Taille, RaisonSUI]));
          nbsui=length(SUI);
          TDEN = lstcat( TDEN, list([ concentration_ponctuelle,2,nbsui,nbp-nbt+i ]));
        end
      else
        printf("Error : unknown error \n"); return;
      end            
    end
    //
    if ((TypeDEN == 'DROI') | (TypeDEN == 'SEGM')) then
      if (TypeDEN == 'DROI') , tyDEN=concentration_axiale; end
      if (TypeDEN == 'SEGM') , tyDEN=concentration_segment; end
      if ( (nbcoord == 2) & ( nbt == 1 ) & ( nbr == 1 )) then
        SUI= lstcat( SUI, list([ suite_geometrique , Taille, RaisonSUI]));
        nbsui=length(SUI);
        TDEN = lstcat( TDEN, list([ tyDEN,3,nbsui,nbp-1,nbp]));     
      else
        printf("Error : unknown error \n"); return;
      end
    end
  end
end
endfunction
  
// ====================================================================
// ======================== APPELS SYSTEME/DELOS ======================
// ====================================================================

function [err,MESH,RES]=DSdelos(GEO,DEN)
// ====================================================================
// OBJET DSdelos : Lancement de DELOS sous Linux
// pour Windows (TODO)
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
//printf("nbvar = %d\n",nbvarin)
if ( nbvarin < 6 ) , colorAllReg= []; end
//
[OS,vers]=getos();
if (( nomden == "" ) | (isempty(nomden))) , nomden="-" ; end
if (( nomres == "" ) | (isempty(nomres))) , nomres="-" ; end
commande="delos v "+nomgeo+"  "+nomden+"  "+nommesh+"  "+nomres+" - ";
err=unix_g(commande);
endfunction

function [message,XY,TRI,REG]=DSdelosfile(nomgeo,nomden,nommesh,nomres)
// ====================================================================
// OBJET DSdelos : Lancement de DELOS sur des fichiers
// pour Windows (TODO)
tempfile=""
XY=[];TRI=[];REG=[]; message="";
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
if ( nbvarin < 4 ) then, nomres=""; end
if ( nbvarin < 3 ) then, nommesh=""; end
if ( nbvarin < 2 ) then, nomden=""; end
//[OS,vers]=getos();
if (isempty(nomden)) then, nomden="-" ; end
if (isempty(nomres)) then, nomres="-" ; end
if (isempty(nommesh)) then
    repTemp=getenv("TMPDIR");
    tempfile=repTemp+"\"+"tmp_delos.mesh"
    nommesh=tempfile ; 
    pause;
end
commande="delos v "+nomgeo+"  "+nomden+"  "+nommesh+"  "+nomres+" - ";
message=unix_g(commande);
if (nbvarout>1) then
  [XY,TRI,REG,GRD]=DSReadMesh(nommesh);    
end
if (~isempty(tempfile)) then, unix("rm "+tempfile);end
endfunction

// ====================================================================
// ========================  POST TRAITEMENT  =========================
// ====================================================================

function DSShowGrd(coord,itrnoe,Grd,numfig,titre,colorAll,Fmesh)
// ====================================================================
// OBJET DSShowGrd : visualise une grandeur sur triangulation sous scilab
// coord : coordonnees des points
// itrnoe(t,:) : les 3 sommets constituant le triangle de numero t
// itreg(i,:): valeur(s) de la grandeur(vecteur/tenseur) au noeud i
// 
[nbvarout,nbvarin]=argn(); //disp(argn(1),argn(2))
//printf("nbvar = %d\n",nbvarin)
if ( nbvarin < 7 ) , Fmesh= []; end
if ( nbvarin < 6 ) , colorAll= []; end
if ( nbvarin < 5 ) , titre = []; end
if ( nbvarin < 4 ) , numfig = []; end
if ( nbvarin < 3 ) , Grd = []; end
if ( nbvarin < 2 ) , itrnoe = []; end
//
if ( isempty(Fmesh)) then, Fmesh=%T;  end
if ( isempty(titre)) then, titre="DSShowGrd";  end
numfig_max = max(winsid()); // Added since v 6.0.0
if ( isempty(numfig_max) ) then numfig_max = 0; end
if ( isempty(numfig)) then, numfig= numfig_max +1;  end                        
//if ( isempty(numfig)) then, numfig=max(winsid())+1;  end                        // ERROR since v 6.0.0
hfig=scf(numfig); // creation d'une figure
// -----------------------------------
if ( ~isempty(colorAll) ) then
  if ( typeof(colorAll) == "string" ) then colorAll=evstr(colorAll); end;
  if ( typeof(colorAll) == "constant" ) then
    [nbcolor,RVB]=size(colorAll);
    if ( RVB == 3 ) then, hfig.color_map=colorAll; end
  end
else
  nbcolor=128;
  xset("colormap",jetcolormap(nbcolor)); 
end
// -----------------------------------
[nbnodes,dim]=size(coord);
[nbt,trois]=size(itrnoe);
num=1:nbt; flag=zeros(1,nbt); 
//
XmM=[min(coord(:,1)),max(coord(:,1))];
YmM=[min(coord(:,2)),max(coord(:,2))];
rect=[XmM;YmM];
//rect=[min(coord(:,1)),max(coord(:,1));min(coord(:,2)),max(coord(:,2))];
[rectm]=DSUTWinMargin(rect,0.1); // marge de 10%
//
triangles=[num', itrnoe, flag'];
// On visualise :
fec(coord(:,1),coord(:,2),triangles,Grd,rect=rectm,strf="040", mesh=Fmesh) ; 
// strf est un code qui signifie ici "sans axes" et "vue isometrique"
// voir : help plol2d
colorbar(min(Grd),max(Grd)); // 
// Il faut trouver le moyen de mettre une legende !!!!!!!!
//
xinfo(string(nbt)+" triangles, "+string(nbnodes)+ " nodes,"+string(nbcolor)+" colors")
endfunction


function [COL,numREG,numCOL]=DSColorTRI(REG,ColorsREG)
// ===========================================================
// Affecte une couleur a chaque triangle en fonction de REG et ColorsREG 
// ? et de la palette courante ?
//
// REG(i) = "region" a laquelle appartient le triangle "i" 
// ou nbt = alors l'affectation des couleurs sera aleatoire.
//
// ColorsREG(1:nbreg) = INDICE des couleurs associees a chaque region
//              les indices des couleurs correspondent au numero dans la palette 
//                     [i,j,k,l...] : indice des couleurs dans la palette
// NON...                    ["red","blue"...] : nom des couleurs dans la palette  
// ATTENTION : il faut que la figure ai ete cree ?!
// NON...                     "rand" : couleurs aleatoires pour chaque region/triangles
//
// COL(i) = couleur du triangle numero "i" 
// Exemples :
//   [COL]=DSColorTRI(nbt) : couleur uniforme pour tous les triangles
//   [COL]=DSColorTRI(nbt,ColorsREG) : couleur aleatoire tire dans ColorsREG
//                                     couleur uniforme si length(ColorsREG)==1
//   [COL]=DSColorTRI(REG,ColorsREG) : numero des regions croissant=> N premieres coul.
// ------------ gestion des couleurs -------------------------
if ( argn(2) < 2 ) then, ColorsREG= []; end
COL=[];numREG=[];numCOL=[];
if ( isempty(REG) ) then return; end;
//
m=length(REG); // m = le nombre de triangles
// 1. Les regions des triangles ne sont pas donnees :
if ( m == 1 ) then 
  m=REG; // m = le nombre de triangles
  if ( isempty(ColorsREG) ) then
    numREG=1; REG=zeros(m,1)+numREG; // ni region, ni couleur => couleur uniforme =1
    numCOL=3; // =color("cyan"); // Mais il faut la palette...
  else
    // pas region mais des couleurs => tirage aleatoire des regions 
    //                                 (meme nombre que les couleurs) :
    numREG=length(ColorsREG); REG=int(rand(m,1)*(numREG-1)+1);
  end
else
  // Recherche des regions
   numREG=gsort(unique(REG(find(REG>0))),'g','i') 
end
// on remplace le numero de la region par la couleur associee
nbreg=length(numREG)
if ( isempty(ColorsREG) ) ,  numCOL=3:nbreg+2; end
COL=REG; // juste pour la taille
for i=1 : nbreg
    COL(find(REG==numREG(i)))=numCOL(i); // on remplace par la couleur
end
// --------------------------- pause;
endfunction