phd-scripts/Unpublished/3D_Voxel_Assign/GrainMask_VUSD_Nodes.for

205 lines
5.9 KiB
Fortran

! This program generates a micro-strucutre mask and applies it to
! an existing mesh.
! J.Grogan 05/08/11
program Corrosion_Preprocessor
!
! Parameters
parameter(max_elements=600000,max_cells=20000,max_fc=100,max_faces=10000)
parameter(max_elem_in_cell=50000,max_neigh_per_elem=6,max_neighbours=50)
!
! Variables
character(len=256)input
character(len=256)input2(2)
integer,allocatable,dimension(:)::num_faces
integer,dimension(max_elements,max_neighbours)::neighbour
integer mesh_type
integer(1) iseed
double precision,allocatable,dimension(:,:)::distance
double precision,allocatable,dimension(:,:)::nbr_dist
double precision,allocatable,dimension(:,:)::cor_faces
double precision,dimension(max_elements)::cor_dist
double precision,allocatable,dimension(:,:)::ele_centroid
!
! Mesh Type: 2 = 2D, 3 = 3D
mesh_type=3
iseed=3
call random_seed(iseed)
!
allocate(ele_centroid(max_elements,3))
! Get element centroids
call elem_centroids(num_elements,ele_centroid,mesh_type)
!
! Get neighbouring elements
allocate(nbr_dist(max_elements,max_neighbours))
do i=1,num_elements
icount=1
do j=1,num_elements
if(icount>max_neighbours)then
rmaxdist=0.
do k=1,max_neighbours
if(nbr_dist(i,k)>rmaxdist)then
rmaxdist=nbr_dist(i,k)
index_max_dist=k
endif
enddo
endif
cent1x=ele_centroid(i,1)
cent1y=ele_centroid(i,2)
cent1z=ele_centroid(i,3)
cent2x=ele_centroid(j,1)
cent2y=ele_centroid(j,2)
cent2z=ele_centroid(j,3)
dist=sqrt((cent1x-cent2x)*(cent1x-cent2x)+(cent1y-cent2y)*(cent1y-cent2y)&
&+(cent1z-cent2z)*(cent1z-cent2z))
if(icount>max_neighbours)then
if(dist<rmaxdist)then
nbr_dist(i,index_max_dist)=dist
neighbour(i,index_max_dist)=j
endif
else
nbr_dist(i,j)=dist
neighbour(i,j)=j
endif
icount=icount+1
enddo
print *,float(i)/float(num_elements)
enddo
deallocate(nbr_dist)
!
! Get corrosion surface distances
open(unit=22,file='CorSurf.dat',status='unknown')
allocate(cor_faces(max_faces,3))
ierr=0
inum_faces=1
do while (ierr==0)
read(22,*,iostat=ierr)cor_faces(inum_faces,1),cor_faces(inum_faces,2),&
&cor_faces(inum_faces,3)
if(ierr==0)inum_faces=inum_faces+1
enddo
close(unit=22)
do i=1,num_elements
centx=ele_centroid(i,1)
centy=ele_centroid(i,2)
centz=ele_centroid(i,3)
distmin=1000.
do j=1,inum_faces-1
facex=cor_faces(j,1)
facey=cor_faces(j,2)
facez=cor_faces(j,3)
dist=sqrt((centx-facex)*(centx-facex)+(centy-facey)*(centy-facey)&
&+(centz-facez)*(centz-facez))
if(dist<distmin)distmin=dist
enddo
cor_dist(i)=distmin
enddo
deallocate(cor_faces)
!
rewind(10)
open(unit=13,file='GeomGenINP.inp',status='unknown')
input2(1)='**'
do while (index(input2(1),'*End Assembly')==0)
read(10,'(a)')input2(2)
write(13,'(a)')input2(1)
input2(1)=input2(2)
enddo
write(13,'(a)')'*End Assembly'
write(13,*)'*INITIAL CONDITIONS,TYPE=SOLUTION'
do i=1,num_elements
call random_number(randnum)
randnum=0
write(13,'(3(a,i6),2(a,f18.6),a,2(i6,a),f18.6,a)')'Assembly.CorPart.',i,',' &
& ,i,',',max_neighbours,',',0.,',',cor_dist(i),',',0,',',1,',',randnum,','
do k=1, max_neighbours
if(mod(k,8)/=0)then
if(k/=max_neighbours)then
write(13,'(i6,a)',advance='no')neighbour(i,k),','
else
write(13,'(i6,a,f18.6,a,i6,a,i6,a,i6)')neighbour(i,k),',',cor_dist(i),',',0,',',0,',',0
endif
else
if(k/=max_neighbours)then
write(13,'(i6,a)')neighbour(i,k),','
else
write(13,'(i6,a,f18.6,a,i6,a,i6,a,i6)')neighbour(i,k),',',cor_dist(i),',',0,',',0,',',0
endif
endif
enddo
enddo
ierr=0
do while (ierr==0)
read(10,'(a)',iostat=ierr)input2(1)
if(ierr==0)write(13,'(a)')input2(1)
enddo
end program
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Get Element Centroids
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
subroutine elem_centroids(ne,ele_centroid,mesh_type)
!
! Parameters
parameter(max_nodes=6000000,max_elements=600000,max_node_in_elem=8)
!
! Variables
character(len=256)input
integer,dimension(max_elements,max_node_in_elem)::elements
double precision,dimension(max_nodes,3)::nodes
double precision,dimension(max_elements,3)::ele_centroid
!
! 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
if(mesh_type==3)then
read(10,*)dummy,nodes(nn,1),nodes(nn,2),nodes(nn,3)
else
read(10,*)dummy,nodes(nn,1),nodes(nn,2)
endif
else
if(mesh_type==2)then
nele_type=4
else
if(index(input,'C3D8')/=0)nele_type=8
if(index(input,'C3D4')/=0)nele_type=4
endif
exit
endif
end do
ne=0
do while(1==1)
read(10,'(a)')input
if(index(input,'*')==0)then
backspace(10)
ne=ne+1
if(nele_type==8)then
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)
elseif(nele_type==4)then
read(10,*)dummy,elements(ne,1),elements(ne,2),elements(ne,3),elements(ne,4)
endif
else
exit
endif
end do
!
! Get element centroid
do i=1,ne
centroidx=0.d0
centroidy=0.d0
centroidz=0.d0
do j=1,nele_type
centroidx=centroidx+nodes(elements(i,j),1)
centroidy=centroidy+nodes(elements(i,j),2)
if(mesh_type==3)centroidz=centroidz+nodes(elements(i,j),3)
enddo
ele_centroid(i,1)=centroidx/float(nele_type)
ele_centroid(i,2)=centroidy/float(nele_type)
if(mesh_type==3)ele_centroid(i,3)=centroidz/float(nele_type)
enddo
end subroutine