octave_packages/secs2d-0.0.8/Utilities/Uise2pde.m

   1 function [mesh,data_v]=Uise2pde(grid_file,pref,data_file_v,load_data,out_file)
   2
   3
   4 ## [mesh,data]=ise2pde3(grid_file,pref,data_file,load_data,out_file)
   5 ## ise2pde3
   6 ## estrae dati dal formato DF-ISE di ISE a pdetool di Matlab
   7 ## grid_file contiene il nome del file di griglia da estrarre
   8 ## pref un prefisso che verra' dato ai files temporanei creati da grep
   9 ## data_file e' un cell array delle file da estrarre
  10 ## load_data e' un cell array che contiene i nomi delle grandezze da estrarre
  11 ## out_file e' il nome del file matlab opzionale per salvare i dati estratti
  12 ##
  13 ## 17-3-2004 ver 3.1
  14 ## Marco Bellini marco_bellini_1@yahoo.it
  15 ## 14.02.2007 ver 3.2
  16 ## Octave porting and bug fixes Carlo de Falco
  17
  18 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  19 ##
  20 ##  CHANGES:
  21 ##   - riconosce i nomi dei contatti e delle regioni
  22 ##   - tratta anche dispositivi con una sola regione
  23 ##   - fornisce informazioni sulla conversione
  24 ##
  25 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  26
  27
  28 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  29 ##
  30 ##  TO DO LOG:
  31 ##
  32 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  33
  34
  35
  36 ##esempio per debug
  37 ## prova a costruire una grid per pdetool
  38 ##data_file="/home/marco/IseDB/jfet/jsdd2_mdr.dat";
  39 ##grid_file="/home/marco/IseDB/jfet/jsdd2_mdr.grd";
  40 ##pref="/home/marco/IseDB/jfet/jsdd2";
  41
  42 ## esempio per struttura onda viaggiante
  43 ## data_file={"/home/marco/IseDB/prova/spice/on1_onda1_tr_des.dat"};
  44 ## grid_file="/home/marco/IseDB/prova/spice/on_mdr.grd";
  45 ## pref="/home/marco/IseDB/prova/spice/on";
  46 ## load_data={"ElectrostaticPotential"};
  47
  48 ## dati da estrarre
  49 ##load_data={"DopingConcentration","BoronActiveConcentration","PhosphorusActiveConcentration"};
  50 ##load_data={"ElectrostaticPotential","eDensity","hDensity","eCurrentDensity","BeamGeneration"};
  51
  52 ##pdeplot(p1,e1,t1,"xydata",log10(abs(data(1,:)))',"zdata",log10(abs(data(1,:))));
  53 ##[p1,e1,t1]=ise2pde(grid_file,pref);
  54
  55 ## leggo i vertici
  56 ## i punti sono ordinati per regione
  57 ## se la prima regione e' l'ossido ed ha 269 punti
  58 ## i primi 269 valori in p1 si riferiscono ai punti dell'ossido
  59 ##
  60 ## nei file di des.cmd di simulazione i punti vengono dati per ogni regione
  61 ## quindi i punti della frontiera tra ossido e silicio vengono ripetuti due
  62 ## volte (tra l'ossido e tra il silicio): e' necessario prenderli una volta
  63 ## sola
  64 ## siccome con piu' di due materiali diversi il discorso e' piu' complicato
  65 ## si sceglie di considerare in questo programma solo un ossido ed il
  66 ## silicio
  67
  68
  69 ## p1 conterra' le coordinate x e y dei vertici
  70 [s,w]=unix(sprintf("grep Vertices %s",grid_file));
  71 n_vert_str = regexp(w,'([0-9]+)',"Tokens");
  72 n_vert=str2num(n_vert_str{1}{1});
  73 unix(sprintf("grep Vertices -A %d %s | tail -n+2 > %s_vertex.txt",n_vert,grid_file,pref));
  74 p1=load(strcat(pref,"_vertex.txt"));
  75 unix(sprintf("rm %s_vertex.txt",pref));
  76 p1(:,2)=-p1(:,2);
  77
  78
  79 fprintf("Found %d vertex\n",n_vert);
  80
  81 ##leggo gli edge
  82 ## el conterra' l'indice dei vertici degli edges degli elementi
  83 ## cioe' l'indice dei vertici dei lati dei triangoli e l'indice dei vertici
  84 ## dei segmenti
  85 [s,w]=unix(sprintf("grep Edges %s",grid_file));
  86 n_edges_str = regexp(w,'([0-9]+)',"Tokens");
  87 n_edges=str2num(n_edges_str{1}{1});
  88 unix(sprintf("grep Edges -A %d %s | tail -n+2  > %s_edges.txt",n_edges,grid_file,pref));
  89 el=load(strcat(pref,"_edges.txt"));
  90 unix(sprintf("rm %s_edges.txt",pref));
  91 fprintf("Found %d edges\n",n_edges);
  92 clear n_edges;
  93 el=el+1;
  94
  95
  96 ##leggo gli elementi triangolari
  97 ## el_tr contiene gli indici degli edge dei triangoli
  98 [s,w]=unix(sprintf("grep Elements %s",grid_file));
  99 n_els_str = regexp(w,'([0-9]+)',"Tokens");
 100 n_els=str2num(n_els_str{1}{1});
 101 ## leggo solo gli elementi che iniziano per 2 e che corrispondono ai
 102 ## triangoli, creando un file che contiene solo gli elementi triangolari
 103 unix(sprintf("grep Elements -A %d %s | head -n %d  | tail -n %d | awk '$1==2 {print $2,$3,$4}'  > %s_elements2.txt",n_els,grid_file,n_els+1,n_els,pref));
 104 el_tr=load(strcat(pref,"_elements2.txt"));
 105 unix(sprintf("rm %s_elements2.txt",pref));
 106 fprintf("Found %d triangular elements out of %d elements\n",length(el_tr),n_els);
 107
 108
 109 ## creo un file che contiene gli elementi "segmenti"
 110 unix(sprintf("grep Elements -A %d %s | head -n %d  | tail -n %d | awk '$1==1 {print $2,$3}'  > %s_elements1.txt",n_els,grid_file,n_els+1,n_els,pref));
 111 ##el_lin=load(strcat(pref,"_elements1.txt"));
 112 ##unix(sprintf("rm %s_elements1.txt",pref));
 113
 114
 115 ## creo un indice che dice se l'elemento e' un triangolo o una linea
 116 ## e lo salvo nel file pref_linee1.txt; ci sara' un 1 se l'elemento i-esimo
 117 ## e' una linea o un 2 se e' un triangolo
 118 unix(sprintf("grep Elements -A %d %s | head -n %d  | tail -n %d  | awk ' {print $1} ' > %s_linee1.txt",n_els,grid_file,n_els+1,n_els,pref));
 119 clear n_els;
 120
 121 ## leggo le regions
 122 ## cosi posso distinguere tra silicio, contatti e ossido
 123 ## al max e' possibile leggere 60 regions
 124 [s,w]=unix(sprintf("grep regions -n %s  | grep -v nb_regions | tr -d []=",grid_file));
 125 w2=w(11:length(w));
 126 n_regions=0;
 127 l2=length(w2);cont=1;num=1;str="";begin=0;
 128 while (cont<l2)
 129     c=w2(cont);
 130     cont=cont+1;
 131     if (c=="""")
 132         if begin==0
 133             begin=1;
 134             str="";
 135         else
 136             begin=0;
 137             regions{num}=str;
 138             num=num+1;  n_regions= n_regions+1;
 139         end; %begin
 140     end; ## "
 141     if ((c~="""") & (c~=" ") & (begin==1) )
 142         str=strcat(str,c);
 143     end; ## aggiungo carattere
 144 end ## while
 145
 146 fprintf("Found %d regions/contacts:\n",n_regions);
 147 for cont=1:n_regions
 148     fprintf("%d: %s\n",cont,char(regions(cont)));
 149 end
 150
 151 ## leggo le locations
 152 ## trova le righe tra cui Ãš compreso Locations
 153 [s,w]=unix(sprintf("grep Locations -n %s  | awk '/Locations/ {print $1}'",grid_file));
 154 loc_u=sscanf(w,"%d:Locations\n");
 155 [s,w]=unix(sprintf("grep Elements -n %s  | awk '/Elements/ {print $1}'",grid_file));
 156 tmp=sscanf(w,"%d:Elements\n");
 157 loc_d=tmp(1);
 158 ## generazione file
 159 unix(sprintf(" head -n %d %s | tail -n+%d | tr [:cntrl:] ' ' | tr -d }| tr ief 012 > %s_loc.txt",...
 160     loc_d-1,grid_file ,loc_u+1,pref) );
 161 loc=load(sprintf("%s_loc.txt",pref));
 162 ## loc contiene
 163 ## 0 per gli elementi interni
 164 ## 1 per gli esterni
 165 ## 2 per la frontiera
 166 unix(sprintf("rm %s_loc.txt",pref));
 167 clear loc_u;clear loc_d;
 168
 169
 170
 171
 172 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 173 ##   CONVERSIONE DELLA GRIGLIA
 174 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 175
 176 n_ed=max(size(el));
 177 n_el=max(size(el_tr));
 178
 179 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 180 ## calcolo di E1
 181 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 182
 183 ## trova gli elementi sul bordo (esterni)
 184 ## del dispositivo (quelli rossi)
 185 ind=find(loc==1);
 186 e1=el(ind,:);
 187
 188 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 189 ## calcolo di T1
 190 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 191
 192 el_si=sign(el_tr)+1;
 193 el_si=sign(el_si);
 194 ## el_si  contiene le posizioni degli elementi positivi
 195 t1=zeros(n_el,4);
 196
 197 ## elementi tutti positivi: vertice di testa
 198 ind=find( el_si(:,:)==1);
 199 t1(ind)=el(el_tr(ind)+1,1);
 200
 201 ## elementi negativi: vertice di coda
 202 ind=find( el_si(:,:)==0);
 203 t1(ind)=el(-el_tr(ind),2);
 204
 205 ## inserisco il numero di subdominio
 206
 207 ## trovo il delimitatore inferiore di una regione
 208 [s,w]=unix(sprintf("grep material -n %s | grep -v materials  | awk '/material/ {print $1 }'",grid_file));
 209 mat_pos=sscanf(w,"%d:\n");
 210 [s,w]=unix(sprintf("grep material -n %s | grep -v materials  | awk '/material/ {print $2 }' |tr -d =",grid_file));
 211 mat_type_lim=findstr(w,"\n");
 212 if isempty(mat_type_lim)
 213   mat_type=w;
 214 else
 215   for ii=1:length(mat_type_lim)-1
 216     mat_type{ii}=w(mat_type_lim(ii)+1:mat_type_lim(ii+1)-1);
 217   end
 218   mat_type{ii+1}=w(mat_type_lim(ii)+1:end);
 219 end
 220
 221 [s,w]=unix(sprintf("grep material -n %s | grep -v materials  | awk '/material/ {print $4 }' |tr -d =",grid_file));
 222 mat_name_lim=findstr(w,"\n");
 223 if isempty(mat_name_lim)
 224   mat_name=w;
 225 else
 226   for ii=1:length(mat_name_lim)-1
 227     mat_name{ii}=w(mat_name_lim(ii)+1:mat_name_lim(ii+1)-1);
 228   end
 229   mat_name{ii+1}=w(mat_name_lim(ii)+1:end);
 230 end
 231
 232
 233 ## leggo il tipo di elemento
 234 ## se e' un materiale i triangoli appartenenti a quel materiale otterranno
 235 ## indice pari all'indice della regione. t1(4,:)' perche' non ho ancora
 236 ## trasposto
 237
 238 ## se e' un contatto, i segmenti appartenenti a quel contatto saranno
 239 ## indicati con un numero in e(5,:)' (non trasposto)
 240
 241 el_type=load(strcat(pref,"_linee1.txt"));
 242 indice_tr=zeros(size(el_type));
 243 ## indice_tr e' 1 nelle posizioni che corrispondono agli elementi triangolari nella lista
 244 ## totale degli elementi nel file *.grd
 245 indice_tr(find(el_type==2))=1;
 246 ## ora faccio la somma cumulativa in modo che nella posizione i-esima del
 247 ## vettore degli elementi totali ci sia l'indice che identifica il triangolo
 248 ## i-esimo in el_tr : ho costruito una tabella di conversione dall'indice di
 249 ## ISE unico per triangoli e segmenti -> indice dei soli triangoli.
 250 indice_tr=cumsum(indice_tr);
 251 clear ind;
 252
 253 ## leggo gli elementi "segmento"
 254 try
 255 el_lin=load(strcat(pref,"_elements1.txt"));
 256 catch
 257 el_lin=-1;
 258 end_try_catch
 259
 260 ## aggiungo 1 perche' gli elementi partono da 0
 261 el_lin=el_lin+1;
 262 unix(sprintf("rm %s_elements1.txt",pref));
 263 ## costruisco una mappa di conversione anche per gli elementi lineari (indicati da 1 in el_tr).
 264 indice_lin=zeros(size(el_type));
 265 indice_lin(find(el_type==1))=1;
 266 indice_lin=cumsum(indice_lin);
 267 ##
 268 ## creo un vettore che indica il tipo di elemento di frontiera
 269 ## 1 per i lati del dispositivo (elementi esterni)
 270 ## 2 per l'interfaccia ossido-silicio
 271 ## da 3 in poi per i contatti
 272
 273 e_fron=ones(max(size(e1)),1);
 274
 275 ## AGGIUNGO LA FRONTIERA
 276 ##
 277 clear ind1;
 278 ind1=find(loc==2);
 279 name_contact{1}="External";
 280 name_contact{2}="I. Frontier";
 281
 282 if (isempty(ind1)==0)
 283 ## aggiungo gli elementi di frontiera trovati al vettore degli edges e1
 284 ## mettendo indice 2
 285     e1=[e1;el(ind1,:)];
 286     e_fron=[e_fron;2*ones(length(ind1),1)];
 287 end
 288
 289 clear el_type
 290
 291 ## numero del contatto
 292 n_contact=3;
 293 ## name_material contiene il nome del materiale i-esimo
 294 ## name_contact contiene il nome del contatto i-esimo
 295
 296
 297 for n=1:n_regions
 298
 299     if strcmp(mat_type{n},"Contact")==0
 300
 301 ## leggo gli elementi che costituiscono una regione
 302         if (n~=n_regions)
 303             unix(sprintf(" head -n %d %s | tail -n+%d |  tr -d [:cntrl:]}  > %s_tmp2.txt", mat_pos(n+1)-2,grid_file,mat_pos(n)+2,pref));
 304         else
 305 ## trattare l'ultima regione che e' diversa
 306             unix(sprintf(" tail -n+%d %s |  tr -d [:cntrl:]}  > %s_tmp2.txt",mat_pos(n)+2, grid_file,pref) );
 307         end ## e' l'ultima regione?
 308         tmp_el=load(strcat(pref,"_tmp2.txt"));
 309         tmp_el=tmp_el+1;
 310 ## trova i triangoli che fanno parte della regione n
 311         t1(indice_tr(tmp_el),4)=n;
 312         name_material{n}=mat_name{n};
 313     else
 314 ## la regione e' un contatto
 315
 316         if (n~=n_regions)
 317             unix(sprintf(" head -n %d %s | tail -n+%d |  tr -d [:cntrl:]}  > %s_tmp2.txt", mat_pos(n+1)-2,grid_file,mat_pos(n)+2,pref));
 318         else
 319 ## trattare l'ultima regione che e' diversa
 320             unix(sprintf(" tail -n+%d %s |  tr -d [:cntrl:]}  > %s_tmp2.txt", mat_pos(n)+2,grid_file,pref) );
 321         end ## e' l'ultima regione?
 322         tmp_el=load(strcat(pref,"_tmp2.txt"));
 323 ## tmp_el contiene gli indici dei vertici che appartengono al
 324 ## contatto: aggiungo 1 perche' ise parte da 0.
 325         tmp_el=tmp_el+1;
 326
 327 ## predo i vertici la cui posizione e' in tmp_el e la converto da
 328 ## dessis a pde
 329         e1=[e1;el_lin(indice_lin(tmp_el),:)];
 330
 331 ## aggiungo il numero del contatto
 332         tmp_el=tmp_el';
 333         e_fron=[e_fron;n_contact*ones(size(tmp_el))];
 334
 335         name_contact{n_contact}=regions(n);
 336         n_contact=n_contact+1;
 337
 338
 339     end ## non e' un contatto
 340 end ## fine scansione regioni
 341
 342 unix(sprintf("rm %s_linee1.txt",pref));
 343
 344
 345 ## trasposizione
 346 t1=t1';
 347 e1=e1';
 348 p1=p1';
 349
 350
 351
 352 ## individuo a che regione appartengono dei vertici
 353 ## fondamentale per leggere i set di dati prodotti dalle simulazioni
 354
 355 tmpt1=[t1(1,:),t1(2,:),t1(3,:);t1(4,:),t1(4,:),t1(4,:)];
 356 stmpt1=sortrows(tmpt1',1)';
 357 dtmpt1=diff(stmpt1(1,:));
 358 ##clear stmpt1;
 359 ind1=find(dtmpt1==1);
 360 ##%%% questa prende i valori dei vertici in ordine solo per TEST
 361 ##%%% reg_vert=[stmpt1(1,ind1),stmpt1(1,length(dtmpt1)+1)];
 362
 363 ## reg_vert contiene la regione a cui il punto appartiene
 364 reg_vert=[stmpt1(2,ind1),stmpt1(2,length(dtmpt1)+1)];
 365 clear tmpt1 stmpt1 dtmpt1 ind1;
 366
 367
 368 ## individuo i vertici appartenenti alla frontiera
 369 ##  le righe precedenti hanno assegnato questi punti
 370 ##  a una o all'altra regione. Nei file delle simulazioni questi
 371 ##  punti sono assegnati casualmente ad una delle regioni secondo
 372 ##  un criterio arbitrario. E' necessario individuare i punti alla
 373 ##  frontiera per tenerne conto.
 374
 375 ## costruisco e1 completo
 376 e1=[e1;zeros(size(e1))];
 377 e1=[e1;e_fron'];
 378
 379 ## trovo i punti di frontiera
 380 ind1=find(loc==2);
 381 ## nel caso la frontiera non esista non fa niente
 382
 383
 384
 385 if (isempty(ind1)~=1)
 386 ## controlla se c'e' una frontiera tra diversi elementi
 387
 388     el=el';
 389     tmpe1=[el(1,ind1),el(2,ind1)];
 390     stmpe1=sortrows(tmpe1',1)';
 391     dtmpe1=diff(stmpe1(1,:));
 392     clear ind1;
 393 ## non e' detto che i vertici di frontiera siano ad incrementi unitari
 394     ind1=find(dtmpe1>0);
 395
 396 ## in teoria dovrei controllare che non ci sia un solo elemento di
 397 ## frontiera
 398 ##    if (length(ind1)>1)
 399     front_vert=[stmpe1(1,ind1),stmpe1(1,length(dtmpe1)+1)];
 400 ##   else
 401 ## c'e' solo 1 punto di frontiera
 402 ##front_vert=[stmpe1(1,ind1)];
 403 ##end
 404
 405
 406
 407 ## inserisco anche nel reg_vert l'informazione se un punto e' di frontiera
 408 ## in questo modo il valore assoluto indica la regione alla quale e'
 409 ## attribuito il valore.
 410 ## i punti con valore negativo sono punti di frontiera
 411
 412     reg_vert(front_vert)=-reg_vert(front_vert);
 413
 414     clear tmpe1 stmpe1 dtmpe1 ind1;
 415 end % esiste la frontiera
 416
 417 ## % % test
 418 ##   ind1=find(reg_vert==2);
 419 ##   plot(p1(1,ind1),p1(2,ind1),'b.')
 420 ##   whos ind1
 421 ##   hold on;
 422 ##   ind1=find(reg_vert==1);
 423 ##   whos ind1
 424 ##   plot(p1(1,ind1),p1(2,ind1),'g.')
 425 ##   plot(p1(1,front_vert),p1(2,front_vert),'ro')
 426 ##   whos front_vert
 427
 428
 429
 430 ## es. testato: ok
 431 ## ind=find(reg_vert==2);
 432 ## plot(p1(1,:),p1(2,:),'b.')
 433 ## hold on;
 434 ## plot(p1(1,ind1),p1(2,ind1),'g.')
 435
 436
 437 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 438 ##         DATASETS             %
 439 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 440
 441 if (exist("data_file_v")~=0)
 442 ## guardo quanti file devo leggere
 443     n_data_file=max(size(data_file_v));
 444
 445
 446     for nf=1:n_data_file
 447
 448         data_file=char(data_file_v(nf));
 449
 450 ## leggo i nomi dei datasets
 451         [s,w]=unix(sprintf("grep datasets %s",data_file));
 452         w2=w(18:length(w)-1);
 453
 454         ndatasets=0;
 455         datasets=cell(30,1);
 456         l2=length(w2);cont=1;num=1;str="";begin=0;
 457         while (cont<l2)
 458             c=w2(cont);
 459             cont=cont+1;
 460             if (c=="""")
 461                 if begin==0
 462                     begin=1;
 463                     str="";
 464                 else
 465                     begin=0;
 466 ## controllo che non ci siano due regioni in cui e'
 467 ## presente lo stesso dato
 468 ##%  if (num==1) || (strcmp(str,datasets(num-1))==0)
 469 ##% tolgo cosi' ho le label con la molteplicita' giusta
 470                     datasets(num)={str};
 471                     num=num+1;  ndatasets= ndatasets+1;
 472 ##%  end % non e' lo stesso dataset
 473                 end; %begin
 474             end; ## "
 475             if ((c~="""") & (c~=" ") & (begin==1) )
 476                 str=strcat(str,c);
 477             end; ## aggiungo carattere
 478         end ## while
 479
 480
 481 ## estrae i valori di tutte le variabili in un unico file
 482         unix(sprintf("grep Values  -A %d %s > %s_tmp.txt",n_vert,data_file,pref));
 483
 484 ## estrae le righe che delimitano i dati da estrarre
 485
 486 ## del contiene le righe sotto le quali si trovano i valori sui vertici
 487 ## dei triangoli (Values)
 488         [s,w]=unix(sprintf("grep Values -n %s_tmp.txt  | awk '/Values/{print $1}' | tr -d : > %s_del.txt",pref,pref));
 489         del=load(sprintf("%s_del.txt",pref));
 490         unix(sprintf("rm %s_del.txt",pref));
 491
 492 ## uel contiene le righe sopra le quali si trovano i valori sui vertici
 493 ## dei triangoli (})
 494         [s,w]=unix(sprintf("grep } -n %s_tmp.txt  | awk '/}/ {print $1}' | tr -d :  > %s_del.txt",pref,pref));
 495         uel=load(sprintf("%s_del.txt",pref));
 496         unix(sprintf("rm %s_del.txt",pref));
 497
 498 ## nval contiene il numero valori sui vertici contenuto in questa
 499 ## sezione del file
 500         [s,w]=unix(sprintf("grep Values %s_tmp.txt |  tr -d  [:alpha:]{\\(\\) > %s_nval.txt",pref,pref));
 501         nval=load(sprintf("%s_nval.txt",pref));
 502         unix(sprintf("rm %s_nval.txt",pref));
 503
 504         nload=max(size(load_data));
 505         data=zeros(nload,n_vert);
 506
 507 ##calcolo gli elementi ripetuti alla frontiera
 508         n_rip=sum(loc==2)+1;
 509
 510         nl=1;nd=1;
 511         while(nl<=nload)
 512             while(nd<=ndatasets)
 513                 if strcmp(load_data(nl),datasets(nd))==1
 514
 515 ## controllo se il numero dati sui vertici e' < del
 516 ## numero totale dei vertici: accade quando c'e' dell'
 517 ## ossido
 518
 519                     if (nval(nd)<n_vert)
 520 ## caso in presenza di ossido
 521
 522 ## ho trovato il set
 523 ## controllo se e' un caso particolare (ho un set di dati per l'ossido e uno per il silicio)
 524                         if strcmp(load_data(nl),"ElectrostaticPotential")==1 | ...
 525                                 strcmp(load_data(nl),"ElectricField")==1  ...
 526                                 | strcmp(load_data(nl),"LatticeTemperature")==1
 527
 528
 529 ## leggo la prima regione
 530                             dup=del(nd)+1;
 531                             c=1; ## trova il delimitatore successivo
 532                             while uel(c)<=dup
 533                                 c=c+1;
 534                             end;
 535                             ddo=uel(c);
 536                             [s,w]=unix(sprintf(" head -n %d %s_tmp.txt | tail -n+%d |  tr -d [:cntrl:]}  > %s_tmp2.txt", ddo, pref,dup,pref));
 537                             tmp=load(strcat(pref,"_tmp2.txt")); ## tmp2 e' il vettore riga dei dati cercati
 538 ## inserisco i primi elementi: ossido, cioe'
 539 ## regione 2: gli elementi nel file .dat (var: tmp) sono
 540 ## ossido interno + frontiera.
 541 ## ind1 e' appunto dato da ossido + frontiera
 542
 543                             ind1=find(abs(reg_vert)==2 | reg_vert<0);
 544
 545                             data(nl,ind1)=tmp;
 546
 547 ## leggo la seconda regione
 548                             nd=nd+1;
 549                             dup=del(nd)+1;
 550                             c=1; ## trova il delimitatore successivo
 551                             while uel(c)<=dup
 552                                 c=c+1;
 553                             end;
 554                             ddo=uel(c);
 555                             [s,w]=unix(sprintf(" head -n %d %s_tmp.txt | tail -n+%d |  tr -d [:cntrl:]}  > %s_tmp2.txt", ddo, pref,dup,pref));
 556                             tmp=load(strcat(pref,"_tmp2.txt")); ## tmp2 e' il vettore riga dei dati cercati
 557 ## inserisco i secondi elementi
 558 ## silicio + frontiera
 559                             clear ind;
 560                             ind=find(abs(reg_vert)==1 | reg_vert<0 );
 561                             data(nl,ind)=tmp;
 562
 563                         else
 564 ## caso non particolare
 565 ## metto a zero i dati della regione mancante
 566                             ind1=find(abs(reg_vert)==2 | reg_vert<0);
 567                             data(nl,ind1)=0;
 568
 569 ## leggo la regione
 570                             dup=del(nd)+1;
 571                             c=1; ## trova il delimitatore successivo
 572                             while uel(c)<=dup
 573                                 c=c+1;
 574                             end;
 575                             ddo=uel(c);
 576                             [s,w]=unix(sprintf(" head -n %d %s_tmp.txt | tail -n+%d |  tr -d [:cntrl:]}  > %s_tmp2.txt", ddo, pref,dup,pref));
 577                             tmp=load(strcat(pref,"_tmp2.txt")); ## tmp2 il vettore riga dei dati cercati
 578
 579 ## inserisco i primi elementi
 580                             ind=find(abs(reg_vert)==1 | reg_vert<0 );
 581                             data(nl,ind)=tmp;
 582
 583 ##n_prev_region=length(tmp);
 584 ##data(nl,1:n_prev_region)=tmp;
 585
 586                         end ## trattamento casi particolari
 587
 588                     else
 589 ## il numero dei dati sui vertici e' = al numero
 590 ## dei vertici: uso lo stesso codice di ise2pde
 591                         dup=del(nd)+1;
 592                         c=1; ## trova il delimitatore successivo
 593                         while uel(c)<=dup
 594                             c=c+1;
 595                         end;
 596                         ddo=uel(c);
 597                         [s,w]=unix(sprintf(" head -n %d %s_tmp.txt | tail -n+%d |  tr -d [:cntrl:]}  > %s_tmp2.txt", ddo, pref,dup,pref));
 598                         data(nl,:)=load(strcat(pref,"_tmp2.txt")); ## tmp2 e' il vettore riga dei dati cercati
 599                     end ## numero dei vertici =
 600                 end %if trovato il set
 601                 nd=nd+1;
 602             end %nd
 603             nd=1;
 604             nl=nl+1;
 605         end %nl
 606
 607         clear nd nl;
 608
 609         data_v(:,:,nf)=data(:,:);
 610     end ## fine lettura n-esimo data_file
 611
 612 end ## if c'e' un data_file
 613 clear n_vert w2 l2 cont num str c;
 614
 615 mesh.p=p1;
 616 mesh.e=e1;
 617 mesh.t=t1;
 618 mesh.materials=name_material;
 619 mesh.contacts=name_contact;
 620
 621 for cont=1:length(name_material)
 622     fprintf("Extracted region %d : %s\n",cont,char(name_material{cont}));
 623 end;
 624 for cont=3:length(name_contact)
 625     fprintf("Extracted contact %d : %s\n",cont,char(name_contact{cont}));
 626 end;
 627
 628
 629 if exist("out_file")
 630     if (exist("data_file_v")~=0)
 631         save (out_file,"mesh","data");
 632         fprintf("mesh and data saved");
 633     else
 634         save (out_file,"mesh");
 635         fprintf("mesh saved");
 636     end;
 637 end
 638
 639 unix(sprintf("rm %s_tmp.txt",pref));
 640 unix(sprintf("rm %s_tmp2.txt",pref));
 641
 642 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 643 ## Cdf: Fix edge matrix and build output structure
 644 ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 645
 646 return
 647 mesh.e(7,end)=0;
 648 mesh = Umeshproperties(mesh);
 649
 650 alledges = [mesh.t([1,2,4],:),mesh.t([2,3,4],:),mesh.t([3,1,4],:)];
 651
 652
 653 for iedge = 1:size(mesh.e,2)
 654     whatedgeL = find((mesh.e(1,iedge)==alledges(1,:)& mesh.e(2,iedge)==alledges(2,:)));
 655     whatedgeR = find((mesh.e(1,iedge)==alledges(2,:)& mesh.e(2,iedge)==alledges(1,:)));
 656
 657     if (length(whatedgeL)==1)
 658         mesh.e(6,iedge)=alledges(3,whatedgeL);
 659     end
 660
 661     if (length(whatedgeR)==1)
 662         mesh.e(7,iedge)=alledges(3,whatedgeR);
 663     end
 664
 665 end
 666
 667 maxedge  = max(mesh.e(5,:));
 668 intedges = find((mesh.e(6,:)~=0)&((mesh.e(7,:)~=0)&((mesh.e(7,:)~=(mesh.e(6,:)) ) )));
 669 mesh.e (5,intedges) = maxedge+1;
 670 mesh.contacts{end+1} = "Interface";
 671
 672 return;