% PROYECTO DERIVADA TOPOLOGICA CON SERBA PARALELO
function detopoviejocasa3
% Nombre del archivo ejecutable SERBAPA
ejecutable='serbapa2.exe';
% Nombre del archivo de entrada
nombre='problema';
entradan=[nombre,'.in'];
% Iteracion en la que nos encontramos
iteracion=1;
indice=num2str(iteracion);
salida=0;
pasos=0;
% Comienzo del bucle principal
%______________________________________________________________________________________
while(salida==0);
pasos++;
archivo=[nombre,indice,'.in'];
% Leo archivo Entrada.in y guardo las CC por coordenadas
fid1=fopen(archivo,'r');
% Titulo del problema
Titulo_problema=fscanf(fid1,'%s',1);
% Numero de puntos del contorno
np= fscanf(fid1,'%d',1);
% Numero de puntos internos
npi= fscanf(fid1,'%d',1);
% Numero de puntos fijos
pf= fscanf(fid1,'%d',1);
% Tipo de problema: tension plana/deformacion plana
ntip= fscanf(fid1,'%d',1);
% Constante "E"
E= fscanf(fid1,'%f',1);
% Constante de Poisson
pois= fscanf(fid1,'%f',1);
% Longitud maxima del contorno
maxlon= fscanf(fid1,'%f',1);
% Coordenada de puntos internos
for i=1:npi;
xint(i) = fscanf(fid1,'%f',1);
yint(i) = fscanf(fid1,'%f',1);
endfor
% Coordenadas de los puntos del contorno
for i=1:np;
nodo= fscanf(fid1,'%f',1);
xcon(i)=fscanf(fid1,'%f',1);
ycon(i)=fscanf(fid1,'%f',1);
endfor
% Condiciones de contorno
for i=1:np;
nodo(i)=fscanf(fid1,'%d',1);
cod1(i)=fscanf(fid1,'%d',1);
valor1(i)=fscanf(fid1,'%f',1);
cod2(i)=fscanf(fid1,'%d',1);
valor2(i)=fscanf(fid1,'%f',1);
cod3(i)=fscanf(fid1,'%d',1);
valor3(i)=fscanf(fid1,'%f',1);
cod4(i)=fscanf(fid1,'%d',1);
valor4(i)=fscanf(fid1,'%f',1);
endfor
% Matriz de conectividades
for i=1:np;
elemento(i)=fscanf(fid1,'%d',1);
nodent(i)=fscanf(fid1,'%d',1);
nodsal(i)=fscanf(fid1,'%d',1);
endfor
% Prescindibilidades
for i=1:np+npi;
nada=fscanf(fid1,'%d',1);
prescin(i)=fscanf(fid1,'%d',1);
endfor
fclose(fid1);
% Llamada a Fortran al programa de resolucion en paralelo SERBAPA
% El archivo de entrada se llamara siempre entrada.in
rename(archivo,'entrada.in');
% Numero de procesadores del cluster en los que se va a ejecutar
procesadores = '3';
comando=['mpirun.mpich -n ', procesadores,' ', ejecutable];
[p,q]=unix(comando);
% SERBAPA da como salida un archivo llamado salida.sal
% Leo archivo salida.sal con la derivada topologica
fid2= fopen('salida.sal','r');
comando=['cp ./salida.sal ./salidas/']
unix(comando)
indice
a='salida_iteracion_'
ar=[a,indice]
cd salidas
rename('salida.sal', ar)
cd ..
% Datos de entrada que ya tengo almacenados
fscanf(fid2,'%s',20);
fscanf(fid2,'%s',15);
for i=1:np;
fscanf(fid2,'%f',3);
endfor
fscanf(fid2,'%s',6);
for i=1:np;
fscanf(fid2,'%s',6);
endfor
fscanf(fid2,'%s',3);
for i=1:np;
fscanf(fid2,'%d',3);
endfor
fscanf(fid2,'%s',6);
for i=1:np;
fscanf(fid2,'%s',6);
endfor
fscanf(fid2,'%s',11);
% Almacenamiento de los desplazamientos de los puntos del contorno
for i=1:np;
fscanf(fid2,'%d',1);
desx(i)=fscanf(fid2,'%f',1);
desy(i)=fscanf(fid2,'%f',1);
endfor
fscanf(fid2,'%s',9);
% Alamcenamiento de las tensiones en el contorno
for i=1:np;
fscanf(fid2,'%d',1);
tenxx(i)=fscanf(fid2,'%f',1);
tenyy(i)=fscanf(fid2,'%f',1);
tenxy(i)=fscanf(fid2,'%f',1);
tenzz(i)=fscanf(fid2,'%f',1);
endfor
fscanf(fid2,'%s',12);
% Alamcenamiento de los desplazamientos de los puntos internos
for i=1:npi;
fscanf(fid2,'%f',1);
fscanf(fid2,'%f',1);
dpix(i)=fscanf(fid2,'%f',1);
dpiy(i)=fscanf(fid2,'%f',1);
endfor
fscanf(fid2,'%s',12);
% Almacenamiento de las tensiones de los puntos internos
for i=1:npi;
fscanf(fid2,'%f',1);
fscanf(fid2,'%f',1);
tens11(i)=fscanf(fid2,'%f',1);
tens22(i)=fscanf(fid2,'%f',1);
tens12(i)=fscanf(fid2,'%f',1);
tens33(i)=fscanf(fid2,'%f',1);
endfor
fscanf(fid2,'%s',15);
% Almacenamiento de las deformaciones del contorno
for i=1:np;
fscanf(fid2,'%d',1);
defxx(i)=fscanf(fid2,'%f',1);
defyy(i)=fscanf(fid2,'%f',1);
defxy(i)=fscanf(fid2,'%f',1);
defzz(i)=fscanf(fid2,'%f',1);
endfor
fscanf(fid2,'%s',17);
% Almacenamiento de las deformaciones de los puntos internos
for i=1:npi;
fscanf(fid2,'%f',1);
fscanf(fid2,'%f',1);
def11(i)=fscanf(fid2,'%f',1);
def22(i)=fscanf(fid2,'%f',1);
def12(i)=fscanf(fid2,'%f',1);
def33(i)=fscanf(fid2,'%f',1);
endfor
fscanf(fid2,'%s',10);
% Almacenamiento de la derivada topologica del contorno
derivada=zeros(1,np+npi);
for i=1:np;
coordenadax(i)=fscanf(fid2,'%f',1);
coordenaday(i)=fscanf(fid2,'%f',1);
derivada(i)=fscanf(fid2,'%f',1);
endfor
% Alacenamiento de la derivada topologica de los puntos internos
for i=1:npi;
coordenadax(i+np)=fscanf(fid2,'%f',1);
coordenaday(i+np)=fscanf(fid2,'%f',1);
derivada(i+np)=fscanf(fid2,'%f',1);
endfor
fscanf(fid2,'%s',4);
% Almacenamiento de la permisibilidad de np + npi
for i=1:(np+npi);
permis(i)=prescin(i);
endfor
fclose(fid2);
% Calculo del area y puntos a eliminar
if(iteracion==1);
puntos_a_eliminar=round((np+npi)*0.011 )
endif
ptosaelim=puntos_a_eliminar;
% Calculo de los puntos que sean prescindibles con DT menor
% Dichos puntos los pongo en una "lista negra"
[valor_DT,posicion]=sort(abs(derivada)); %Ordenamos derivada en orden creciente
lista_negra=zeros(1,np+npi); %Contiene los puntos "eliminables"
for i=1:(np+npi);
if(ptosaelim>0)
lista_negra(posicion(i))=1;
ptosaelim--;
endif
endfor
%-----------------------------------------------------------------------------------
% Inicializacion de bucle para quitar elementos dobles
elimina=1;
total1=0;
total2=0;
% Bucle para quitar elementos dobles y partes sueltas
while(elimina!=0)
iteracion++;
indice=num2str(iteracion);
% Genero un archivo auxiliar con la nube de puntos: nube.dat
% que es del tipo:
% num_puntos nodos x y
% coord_x coord_y
nuben=['nube',indice];
nube=[nuben,'.dat'];
fid3=fopen(nube,'w');
fprintf(fid3,"%d nodos x y \n",np+npi-puntos_a_eliminar);
for i=1:np;
if(lista_negra(i)==0);
fprintf(fid3,"%f %f\n",xcon(i),ycon(i));
total1++;
endif
endfor
for i=1:npi;
if(lista_negra(i+np)==0);
fprintf(fid3,"%f %f\n",xint(i),yint(i));
total2++;
endif
endfor
fclose(fid3);
printf("np %d, npi %d, ptsaelim %d\n",np,npi,puntos_a_eliminar);
resta=np+npi-puntos_a_eliminar-total1-total2;
printf("total1 %d, total2 %d, resta %d\n",total1,total2,resta);
resta=0;
total1=0;
total2=0;
% Llamada al meshsuite desde el suprograma call_msuite:
% que es el que crea el archivo Entrada.mshscr y llama al msuite.exe
call_msuite(nuben);
% Llamada al lee_msuite que plotea la geometria y devuelve conectividades:
[coord,conect,coordint,cint]=lee_msuite2(nuben);
[filanpi,col]=size(coordint);
[filanp,col2]=size(coord);
% Comprobacion de si hay un punto en doble frontera
elimina=0;
k=0;
for i=1:filanp
for j=1:(i-1)
if(conect(i,1)==conect(j,1))
elimina=1;
k++;
nodos(k)=conect(i,1);
printf("problemas en %d\n",conect(i,1));
endif
endfor
endfor
% Si hay problemas en un nodo, este se elimina
if(elimina==1)
for h=1:k
for j=1:np
if((xcon(j)==coord(nodos(h),1))&...
(ycon(j)==coord(nodos(h),2))&...
(lista_negra(j)==0))
lista_negra(j)=1;
puntos_a_eliminar++;
% Elimino un punto del contorno
printf("coordenadas %f %f\n",coord(nodos(h),1),coord(nodos(h),2));
endif
endfor
for j=1:npi
if((xint(j)==coord(nodos(h),1))&...
(yint(j)==coord(nodos(h),2))&...
(lista_negra(j+np)==0))
lista_negra(j+np)=1;
puntos_a_eliminar++;
% Elimino un punto interno
printf("coordenadas2 %f %f\n",coord(nodos(h),1),coord(nodos(h),2));
endif
endfor
endfor
endif
% Comprobacion de que no hay bloques desconectados
% Leo cada uno de los agujeros o trozos sueltos
fid9=0;
trozo=1;
if(elimina==0)
while fid9~=-1
% Abro el archivo que contiene al agujero
abrilo=[nuben,'_',num2str(trozo),'.con']
fid9=fopen(abrilo,'r');
if(fid9==-1)
continue
endif
eofstat = feof(fid9);
j=1;
nodor(j)= fscanf(fid9,'%d',1);
nadar(j)= fscanf(fid9,'%d',1);
% Establezco un punto q supongo dentro del dominio
% Calculo la normal al primer elemento
n1=cint(nadar(j),2)-cint(nodor(j),2);
n2=-cint(nadar(j),1)+cint(nodor(j),1);
% Y la tangente
t1=cint(nadar(j),1)-cint(nodor(j),1);
t2=cint(nadar(j),2)-cint(nodor(j),2);
% Un punto de dentro tendra coordenadas:
xcero1=cint(nodor(j),1)
ycero1=cint(nodor(j),2)
xcero=cint(nodor(j),1)+0.5*t1-0.2*n1;
ycero=cint(nodor(j),2)+0.5*t2-0.2*n2;
suma1=nodor(j);
suma2=nadar(j);
% Sumo los nodos y seran iguales cuando haya dado una vuelta completa
while(suma1-suma2!=0)
j++;
nodor(j)=fscanf(fid9,'%d',1);
nadar(j)=fscanf(fid9,'%d',1);
suma1+=nodor(j);
suma2+=nadar(j);
endwhile
% Compuebo el sentido q tienen
suma_total=0;
suma_angulo=0;
alfa1=0;
alfa2=0;
% Calculo el angulo q forman
for i=1:j
alfa1=atan2(cint(nodor(i),2)-ycero,...
cint(nodor(i),1)-xcero);
alfa2=atan2(cint(nadar(i),2)-ycero,...
cint(nadar(i),1)-xcero);
% si alguno de los dos es neativo, entonces ............
% aqui esta la clave
suma_angulo=(alfa2-alfa1)*180/pi;
if(suma_angulo<-180)
suma_angulo=360+suma_angulo;
endif
if(suma_angulo>180)
suma_angulo=-360+suma_angulo;
endif
suma_angulo;
suma_total=suma_total+suma_angulo;
endfor
k=j
suma_total;
% En el caso de ir a contrareloj
% poner en todos los nodos eliminar(i)=1
diferencia=abs(suma_total-360)
if(diferencia<10)
printf("lo quito");
elimina=1;
for i=1:k
for j=1:np
if((xcon(j)==cint(nodor(i),1))&...
(ycon(j)==cint(nodor(i),2))&...
(lista_negra(j)==0))
lista_negra(j)=1;
puntos_a_eliminar++;
endif
endfor
for j=1:npi
if((xint(j)==cint(nodor(i),1))&...
(yint(j)==cint(nodor(i),2))&...
(lista_negra(j+np)==0))
lista_negra(j+np)=1;
puntos_a_eliminar++;
endif
endfor
endfor
endif
close(fid9)
trozo++;
endwhile
elimina;
endif
% Borra archivos de nubes para q no haya tantos dando la lata
unix('rm nube*.con');
unix('rm nube*.meshscr');
unix('rm nube*.out');
endwhile
% Fin de bucle de ptos dobles y trozos sueltos
%-----------------------------------------------------------------------------------------
% Nombre de archivo salida =entradaxx
% Creo el archivo de entrada Entrada2.in y vuelvo al principio.
indice=num2str(iteracion);
entradan=[nombre,indice,'.in'];
fid4=fopen(entradan,'w');
% Escribimos el archivo de entrada de la siguiente iteracion
fprintf(fid4,"%s\n",Titulo_problema);
fprintf(fid4,"%d %d %d\n",filanp,filanpi,pf);
fprintf(fid4,"%d %f %f %f\n",ntip,E,pois,maxlon);
% Establecemos coordenadas de puntos internos
for i=1:filanpi;
fprintf(fid4,"%f %f\n",coordint(i,1),coordint(i,2));
endfor
% Establecemos coordenadas del contorno
for i=1:filanp;
fprintf(fid4,"%d %f %f\n",i,coord(i,1),coord(i,2));
endfor
% Por defecto todas las CC estan a 5678=0
for i=1:filanp;
cod1p(i)=5;
valor1p(i)=0;
cod2p(i)=6;
valor2p(i)=0;
cod3p(i)=7;
valor3p(i)=0;
cod4p(i)=8;
valor4p(i)=0;
endfor
% Inicializo el vector de permisibilidad
for i=1:(filanpi+filanp);
permisp(i)=0;
endfor
% Pongo CC por coordenadas
for i=1:np;
% printf("\n%d %d",i, permis(i))
if (permis(i)==2)
for j=1:filanp;
if((coord(j,1)==xcon(i))&(coord(j,2)==ycon(i)))
cod1p(j)=cod1(i);
valor1p(j)=valor1(i);
cod2p(j)=cod2(i);
valor2p(j)=valor2(i);
cod3p(j)=cod3(i);
valor3p(j)=valor3(i);
cod4p(j)=cod4(i);
valor4p(j)=valor4(i);
permisp(j)=2;
endif
endfor
endif
if (permis(i)==1)
for j=1:filanp;
if((coord(j,1)==xcon(i))&(coord(j,2)==ycon(i)))
permisp(j)=1;
endif
endfor
for j=1:filanpi;
if((coordint(j,1)==xcon(i))&...
(coordint(j,2)==ycon(i)))
permisp(j+filanp)=1;
endif
endfor
endif
endfor
% Establecemos CC
for i=1:filanp;
fprintf(fid4,"%d %d %f %d %f %d %f %d %f\n",...
i,cod1p(i),valor1p(i),cod2p(i),valor2p(i),...
cod3p(i),valor3p(i),cod4p(i),valor4p(i));
endfor
% Ponemos a uno los puntos no permitidos del contorno
for i=1:npi;
if(permis(i)==1)
for j=1:filanpi
if((coordint(j,1)==xint(i))&(coordint(j,2)==yint(i)))
permisp(j+filanp)=1;
endif
endfor
endif
endfor
% Escribimos matriz de conectividad
for i=1:filanp;
fprintf(fid4,"%d %d %d\n",i,conect(i,1),conect(i,2));
endfor
% Escribimos prescindibilidad
for i=1:(filanp+filanpi);
fprintf(fid4,"%d %d\n",i,permisp(i));
endfor
% Cerramos el archivo de salida
fclose(fid4);
endwhile
end
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
