
%	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

