c************************** burn3d.f *******************************
c PROBLEM DEFINITION
c In a random multi-phase structure, a question that is important
c is whether a particular phase percolates through the microstructure
c or not. This program is designed to answer that question, for
c a general 3-D multi-phase random microstructure. The burning
c algorithm checks whether a percolation threshold exists in a periodic
c image. The program searches all three directions, using periodic
c boundary conditions in the two perpendicular directions for
c each burn.
c Variables
c There are a maximum of "ntot" phases possible, numbered 1,2,3,...
c The label of the phase being burned is "phase", and is an input variable.
c The value assigned to burned pixels is "burned" and is equal to ntot+1.
c To burn on more than one phase at a time, just use
c "phase2", "phase3", etc., and check for these values too, whenever
c the value "phase" is checked for. (See manual)
c A value of 0 is assigned to the variables percx, percy, and
c percz for non-continuity, and 1 for percolation (continuity) in the
c given direction.
c Dimensions
c (USER) The variable pix is dimensioned the size of the system.
c The variables old and new are dimensioned 1/10 the size of the system,
c but with three components. (A minimum of 1000 should be used
c to dimension these variables, so that small systems will have enough
c computation room.) These array dimensions should be changed
c simultaneously in all subroutines using a global replacement statement.
c (USER) Dimensions of main arrays
integer*2 pix(1000000),old(100000,3),new(100000,3)
integer*2 burned,phase
c (USER) Unit = 9 is the input file, unit 7 is the output file
open(unit=9,file='microstructure.dat')
open(unit=7,file='output.out')
c (USER) system size ns = nx x ny x nz
nx=100
ny=100
nz=100
ns=nx*ny*nz
c (USER) Identify the label of the phase to be burned
phase=1
c (USER) Total number of phases allowed in problem
ntot=100
c Label of burned pixel
burned=ntot+1
c Read in microstructure file
do 330 k=1,nz
do 330 j=1,ny
do 330 i=1,nx
m=nx*ny*(k-1)+nx*(j-1)+i
read(9,*) pix(m)
330 continue
c Call the subroutine that actually does the burning
call fire(pix,nx,ny,nz,percz,percy,percx,new,old,phase,burned)
c Output the values of perc*, showing the continuity of the three
c principal directions
write(7,*) ' percx = ',percx
write(7,*) ' percy = ',percy
write(7,*) ' percz = ',percz
end
c This subroutine does the actual burning. The burning starts with the
c k=1 or j=1 or i=1 plane, and then iteratively burnes through the system,
c until there are no more acessible pixels to be burned. To be burned,
c a pixel must be next to a previously burned pixel.
subroutine fire(pix,nx,ny,nz,percz,percy,percx,new,old,
& phase,burned)
integer*2 pix(1000000)
integer*2 old(100000,3),new(100000,3)
integer*2 in(6),jn(6),kn(6),phase,burned
c System size
ns=nx*ny*nz
c Direction labels to check for burning path (nearest neighbor information
c in 3-D digital system).
in(1)=-1
in(2)=1
in(3)=0
in(4)=0
in(5)=0
in(6)=0
jn(1)=0
jn(2)=0
jn(3)=-1
jn(4)=1
jn(5)=0
jn(6)=0
kn(1)=0
kn(2)=0
kn(3)=0
kn(4)=0
kn(5)=-1
kn(6)=1
c Initialize percolation flags
percx=0.0
percy=0.0
percz=0.0
do 3000 ijk=1,3
c Build up first burned pixels from i or j or k=1 layer, according to
c choice of ijk (ijk = 1, k = 1; ijk = 2, j = 1; ijk = 3, i = 1).
c Store(i,j,k) labels in array old().
iold=0
if(ijk.eq.1) then
n2=ny
n1=nx
end if
if(ijk.eq.2) then
n2=nz
n1=nx
end if
if(ijk.eq.3) then
n2=ny
n1=nz
end if
do 1000 jj=1,n2
do 1000 ii=1,n1
if(ijk.eq.1) then
i=ii
j=jj
k=1
end if
if(ijk.eq.2) then
i=ii
j=1
k=jj
end if
if(ijk.eq.3) then
i=1
j=jj
k=ii
end if
m=nx*ny*(k-1)+nx*(j-1)+i
if(pix(m).eq.phase) then
pix(m)=burned
iold=iold+1
old(iold,1)=i
old(iold,2)=j
old(iold,3)=k
end if
1000 continue
c If no pixels burned in first layer, then phase can't possibly percolate,
c so move to next direction
if(iold.eq.0) goto 3000
c Now start building up new burned pixels from old set of burned pixels,
c thus propagating the fire.
60 inew=0
do 100 i=1,iold
ii=old(i,1)
jj=old(i,2)
kk=old(i,3)
c check all six nearest neighbors of previously burned pixel
do 90 n=1,6
i1=ii+in(n)
j1=jj+jn(n)
k1=kk+kn(n)
c Periodic boundary conditions
c (USER) Can replace with hard boundary conditions
c if desired to remove periodicity. Keep the ijk if statements and the
c goto 90 statements, and change the other if statements to have
c goto 90 as well. That way the program does not allow "wrappping"
c around to find a neighbor. (See manual)
if(ijk.eq.1) then
if(k1.lt.1.or.k1.gt.nz) goto 90
if(i1.lt.1) i1=i1+nx
if(i1.gt.nx) i1=i1-nx
if(j1.lt.1) j1=j1+ny
if(j1.gt.ny) j1=j1-ny
end if
if(ijk.eq.2) then
if(j1.lt.1.or.j1.gt.ny) goto 90
if(i1.lt.1) i1=i1+nx
if(i1.gt.nx) i1=i1-nx
if(k1.lt.1) k1=k1+nz
if(k1.gt.nz) k1=k1-nz
end if
if(ijk.eq.3) then
if(i1.lt.1.or.i1.gt.nx) goto 90
if(k1.lt.1) k1=k1+nz
if(k1.gt.nz) k1=k1-nz
if(j1.lt.1) j1=j1+ny
if(j1.gt.ny) j1=j1-ny
end if
c Store (i,j,k) labels of newly burned pixels in array new().
m1=nx*ny*(k1-1)+nx*(j1-1)+i1
if(pix(m1).eq.phase) then
pix(m1)=burned
inew=inew+1
new(inew,1)=i1
new(inew,2)=j1
new(inew,3)=k1
end if
90 continue
100 continue
c If new pixels were burned, then transfer labels to old() array, start
c burning process over again.
if(inew.gt.0) then
iold=inew
do 150 i=1,inew
old(i,1)=new(i,1)
old(i,2)=new(i,2)
old(i,3)=new(i,3)
150 continue
goto 60
end if
c If no new burned pixels, then check to see if the last layer of the image
c has any burned pixels in it. If so, then there is continuity. If not,
c then there is no continuity.
if(ijk.eq.1) then
n2=ny
n1=nx
end if
if(ijk.eq.2) then
n2=nz
n1=nx
end if
if(ijk.eq.3) then
n2=ny
n1=nz
end if
do 30 jj=1,n2
do 30 ii=1,n1
if(ijk.eq.1) then
i=ii
j=jj
k=nz
end if
if(ijk.eq.2) then
i=ii
j=ny
k=jj
end if
if(ijk.eq.3) then
i=nx
j=jj
k=ii
end if
m=nx*ny*(k-1)+nx*(j-1)+i
if(pix(m).eq.burned) then
if(ijk.eq.1) percz=1.0
if(ijk.eq.2) percy=1.0
if(ijk.eq.3) percx=1.0
end if
30 continue
c Restore burned pixels back to their original label
call restore(pix,ns,phase,burned)
3000 continue
return
end
c This subroutine restores the burned pixels back to their original, unburned
c value (phase).
subroutine restore(pix,ns,phase,burned)
integer*2 pix(1000000),phase,burned
do 10 m=1,ns
if(pix(m).eq.burned) pix(m)=phase
10 continue
return
end