!This program calculates the distance from each element
!to its nearest grain boundary
Program GeomConnectivity
!
!	Parameters
	parameter(max_nodes=1000000,max_elements=1000000,max_node_in_elem=8)
	parameter(max_cells=1000,max_elem_in_cell=10000,max_neigh_per_elem=6)
	parameter(max_faces_in_cell=100)
!
!	Declare Variables
	character*256 input,input2(2)
    integer elements(max_elements,max_node_in_element)
    integer active(max_elements)
	integer num_elements(max_cells)
	integer elem_label=(max_cells,max_elem_in_cell)
	integer num_neighbours=(max_cells,max_elem_in_cell)
	integer elem_neighbours=(max_cells,max_elem_in_cell,max_neigh_per_elem)
	integer num_faces(max_cells)
	real nodes(max_nodes,3)
	real face_coords(max_cells,max_faces_in_cell,9)
	real dotprod(max_node_in_element,max_faces_in_cell)
	real average(max_faces_in_cell)
	real distance(max_cells,max_elem_in_cell)
!
!	Initialise Variables
	nodes=0.
	face_coords=0.
	dotprod=0.
	average=0.
	distance=0.
	elements=0
	active=0
	num_elements=0
	elem_label=0
	num_neighbours=0
	elem_neighbours=0
	num_faces=0
!	
!	Open Input File	and read node and element co-ordinates
    open(unit=10,file='GeomGenTemp.inp',status='unknown') 
    do while (index(input,'*Node')==0)
		read(10,'(a)')input
	end do
    nn=0
	do while(1==1)
		read(10,'(a)')input								
	 	if(index(input,'*')==0)then
			backspace(10)
			nn=nn+1	
			read(10,*)nodes(nn,1),nodes(nn,2),nodes(nn,3)
		else            	
			exit
		endif
	end do
    ne=0
	do while(1==1)
		read(10,'(a)')input								
	 	if(index(input,'*')==0)then
			backspace(10)
			ne=ne+1	
			read(10,*)dummy,elements(ne,1),elements(ne,2),elements(ne,3),elements(ne,4)&
				& ,elements(ne,5),elements(ne,6),elements(ne,7),elements(ne,8)	 		
		else
			exit
		endif
	end do
!	
!	Open GeomGen output file and read element connectivity, cell, face and
!	vertice data.
    open(unit=11,file='vertout.dat',status='old')
    read(11,*)idimension
	read(11,*)num_cells
	do i=1,num_cells
		read(11,*)num_elements(i)
		do j=1,num_elements(i)
			read(11,*)elem_label(i,j)
			read(11,*)num_neighbours(i,j)
			do k=1,num_neighbours(i,j)
				read(11,*)elem_neighbours(i,j,k)
			enddo
		enddo
		read(11,*)num_faces(i)
		do j=1,num_faces(i)
			read(11,*)face_coords(i,j,1),face_coords(i,j,2),face_coords(i,j,3)
			read(11,*)face_coords(i,j,4),face_coords(i,j,5),face_coords(i,j,6)
			read(11,*)face_coords(i,j,7),face_coords(i,j,8),face_coords(i,j,9)
		enddo
	enddo
!
!	Determine normal distance between element centroid and closest face
	do i=1,num_cells
		do j=1,num_elements(i)
			do k=1,8
				nodex=nodes(elements(elem_label(i,j),k),1)
				nodey=nodes(elements(elem_label(i,j),k),2)
				nodez=nodes(elements(elem_label(i,j),k),3)
				do m=1,num_faces(i)
					vert1x=face_coords(i,m,1)
					vert1y=face_coords(i,m,1)
					vert1z=face_coords(i,m,1)
					vert2x=face_coords(i,m,2)
					vert2y=face_coords(i,m,2)
					vert2z=face_coords(i,m,2)
					vert3x=face_coords(i,m,3)
					vert3y=face_coords(i,m,3)
					vert3z=face_coords(i,m,3)
					v1v2i=vert1x-vert2x
					v1v2j=vert1y-vert2y
					v1v2k=vert1z-vert2z
					v1v3i=vert1x-vert3x
					v1v3j=vert1y-vert3y
					v1v3k=vert1z-vert3z
					v1ni=vert1x-nodex
					v1nj=vert1y-nodey
					v1nk=vert1z-nodez
                    crossi=v1v2j*v1v3k-v1v2k*v1v3j
                    crossj=v1v2k*v1v3i-v1v2i*v1v3k
                    crossk=v1v2i*v1v3j-v1v2j*v1v3i
                    dotprod(k,m)=v1ni*crossi+v1nj*crossj+v1nk*crossk
				enddo
			enddo
			min_average=1000.
			do m=1,num_faces(i)
				average(m)=0.
				do k=1,8
					average(m)=average(m)+dotprod(k,m)
				enddo
				average(m)=average(m)/8.
				if(average(m)<min_average)min_average=average(m)
			enddo
			distance(i,j)=min_average
		enddo
	enddo
!	
!	Copy new data into updated input file	
    rewind(unit=10)
    open(unit=12,file='GeomGenINP.inp',status='unknown')
    input2(1)='**'
    do while (index(input2(1),'*End Assembly')==0)
		read(10,'(a)')input2(2)
        write(12,'(a)')input2(1)
        input2(1)=input2(2)
	enddo
    write(12,'(a)')'*End Assembly' 
    write(12,*)'*INITIAL CONDITIONS,TYPE=SOLUTION'
    open(unit=13,file='ecor.dat',status='old')
    active=0
    ierr=0
    do while (ierr==0)
       read(13,*,iostat=ierr)ielnum
       if(ierr==0)active(ielnum)=1   
    enddo  
    do i=1,num_cells
		do j=1,num_elements(i)
			write(12,'(a,8(i6,a))')'Assembly.CorPart.',elem_label(i,j),',',elem_neighbours(i,j,1),&
				& ',',elem_neighbours(i,j,2),',',elem_neighbours(i,j,3),',',elem_neighbours(i,j,4),&
				& ',',elem_neighbours(i,j,5),',',elem_neighbours(i,j,6),',',0,','
			write(12,'(2(i6,a),f18.6,a,5(i6,a))')elem_label(i,j),',',active(elem_label(i,j)),',',&
				& distance(i,j),',',0,',',0,',',0,',',0,',',0,',' 
			write(12,'(8(i6,a))')0,',',0,',',0,',',0,',',1,',',0,',',0,',',0,','
			write(12,'(6(i6,a),f18.6)')0,',',0,',',0,',',0,',',0,',',0,',',0.
		enddo
	enddo 
    ierr=0
    do while (ierr==0)
		read(10,'(a)',iostat=ierr)input2(1)
        if(ierr==0)write(12,'(a)')input2(1)          
	enddo
    close(unit=10)
    close(unit=11)
    close(unit=12)
    close(unit=13)
End Program