c  *************************  dc3d.f  *********************************
c  BACKGROUND

c  This program accepts as input a 3-d digital image, converting it
c  into a real conductor network. The conjugate gradient method
c  is used to solve this finite difference representation of Laplace's 
c  equation for real conductivity problems.
c  Periodic boundary conditions are maintained.
c  In the comments below, (USER) means that this is a section of code
c  that the user might have to change for his particular problem.
c  Therefore the user is encouraged to search for this string.

c  PROBLEM AND VARIABLE DEFINITION

c  The mathematical problem that the conjugate gradient algorithm solves
c  is the minimization of the quadratic form 1/2 uAu, where
c  u is the vector of voltages, and A is generated from the bond
c  conductances between pixels. Nodes are thought of as being in the 
c  center of pixels. The minimization is constrained by maintaining an
c  general direction applied electric field across the sample.
c  The vectors gx,gy,gz are bond conductances, u is the voltage array, 
c  and gb,h, and Ah are auxiliary variables, used in subroutine dembx.
c  The vector pix contains the phase labels for each pixel.
c  The small vector a(i) is the volume fraction
c  of the i'th phase, and currx, curry, currz are the total volume-averaged
c  currents in the x,y, and z directions.

c  DIMENSIONS

c  The vectors gx,gy,gz,u,gb,h,Ah,list,pix are all dimensioned 
c  ns2 = (nx+2)*(ny+2)*(nz+2).  This number is used, rather than the
c  system size nx x ny x nz, because an extra layer of pixels is 
c  put around the system to be able to maintain periodic boundary
c  conditions (see manual, Sec. 3.3). The arrays pix and list are also
c  dimensioned this way.
c  At present the program is set up for up to 100 
c  phases, but that can easily be changed by the user, by changing the
c  dimension of sigma, a, and be. Note that be has both dimensions
c  equal to each other. The parameter nphase gives the number of
c  phases being considered. The parameter ntot is the total number
c  of phases possible in the program, and should be equal to the 
c  dimension of sigma, a, and be.
c  All arrays are passed to subroutines in the call statements.

c  STRONGLY RECOMMENDED:  READ MANUAL BEFORE USING THE PROGRAM!!

c  (USER)  Change these dimensions for different system sizes.  All
c  dimensions in the subroutines are passed, so do not need to be changed.
c  The dimensions of sigma, a, and be should be equal to the value of ntot.
      real gx(10648),gy(10648),u(10648),gz(10648)
      real gb(10648),h(10648),Ah(10648)
      real currx,curry,currz,sigma(100,3)
      real a(100),be(100,100,3)
      integer*2 pix(10648)
      integer*4 list(10648)

c  (USER) Unit 9 is the microstructure input file, unit 7 is the 
c  results output file.
      open(unit=9,file='microstructure.dat')
      open(unit=7,file='outputfile.out')

c  (USER) real image size is nx x ny x nz
      nx=20
      ny=20
      nz=20
      write(7,1111) nx,ny,nz,nx*ny*nz
1111  format(' Image is ',3i6,' No. of real sites = ',i8)
c  auxiliary variables involving the system size
      nx1=nx+1
      ny1=ny+1
      nz1=nz+1
      nx2=nx+2
      ny2=ny+2
      nz2=nz+2
      L22=nx2*ny2
c  computational image size ns2 is nx2 x ny2 x nz2
      ns2=nx2*ny2*nz2

c  (USER) set cutoff for norm squared of gradient, gtest.  gtest is
c  the stopping criterion, compared to gb*gb. When gb*gb is less
c  than gtest=abc*ns2, then the rms value of the gradient at a pixel
c  is less than sqrt(abc).
      gtest=1.0e-16*ns2

c  (USER) nphase is the number of phases being considered in the problem.
c  The values of pix(m) will run from 1 to nphase. ntot is the total 
c  number of phases possible in the program, and is the dimension of
c  sigma, a, and be.
      nphase=2
      ntot=100

c  Make list of real (interior) sites, used in subroutine dembx.  The 1-d
c  labelling scheme goes over all ns2 sites, so a list of the real sites
c  is needed.
      nlist=0
      do 103 i=2,nx1
      do 102 j=2,ny1
      do 101 k=2,nz1
      m=i+(j-1)*nx2+(k-1)*L22
      nlist=nlist+1
      list(nlist)=m
101   continue
102   continue
103   continue

c  Compute average current in each pixel.
c (USER) npoints is the number of microstructures to use.

      npoints=1
      do 8000 micro=1,npoints
c Read in a microstructure in subroutine ppixel, and set up pix(m)
c  with the appropriate phase assignments.
      call ppixel(pix,nx2,ny2,nz2,a,ns2,nphase,ntot)
c  output phase volume fractions
      do 99 i=1,nphase
      write(7,299) i,a(i)
299   format(' Phase fraction of ',i3,' = ',f12.6)
99    continue

c  (USER) Set components of applied field, E = (ex,ey,ez)
      ex=1.0
      ey=1.0
      ez=1.0
      write(7,*) 'Applied field components:'
      write(7,*) 'ex = ',ex,' ey = ',ey,' ez = ',ez
c  Initialize the voltage distribution by putting on uniform field.
      do 30 k=1,nz2
      do 30 j=1,ny2
      do 30 i=1,nx2
      m=(k-1)*nx2*ny2+nx2*(j-1)+i
      u(m)=-ex*i-ey*j-ez*k
30    continue
c  (USER) input value of real conductivity tensor for each phase
c  (diagonal only). 1,2,3 = x,y,z, respectively.
      sigma(1,1)=1.0
      sigma(1,2)=1.0
      sigma(1,3)=1.0
      sigma(2,1)=0.5
      sigma(2,2)=0.5
      sigma(2,3)=0.5

c  Subroutine bond sets up conductor network in gx,gy,gz 1-d arrays
       call bond(pix,gx,gy,gz,nx2,ny2,nz2,ns2,sigma,be,nphase,ntot)

c  Subroutine dembx accepts gx,gy,gz and solves for the voltage field
c  that minimizes the dissipated energy.  
      call dembx(nx2,ny2,nz2,ns2,gx,gy,gz,u,ic,gb,h,Ah,list,nlist,gtest)

c  find final current after voltage solution is done
      call current(nx2,ny2,nz2,ns2,currx,curry,currz,u,gx,gy,gz)
      write(7,*) 'Average current in x direction= ',currx
      write(7,*) 'Average current in y direction= ',curry
      write(7,*) 'Average current in z direction= ',currz
      write(7,*) ic,' number of conjugate gradient cycles needed'
      call flush(7)
8000  continue
      end

c  Subroutine that performs the conjugate gradient solution routine to 
c  find the correct set of nodal voltages

      subroutine dembx(nx2,ny2,nz2,ns2,gx,gy,gz,u,ic,gb,h,Ah,
     &  list,nlist,gtest)
      real gx(ns2),gy(ns2),u(ns2),gb(ns2)
      real Ah(ns2),h(ns2),gz(ns2)
      real gg,hAh,lambda,gglast,gamma,ravg,currx,curry,currz
      integer*4 list(ns2),ncheck

c  Note:  voltage gradients are maintained because in the conjugate gradient
c  relaxation algorithm, the voltage vector is only modified by adding a
c  periodic vector to it.
      L22=nx2*ny2
c  First stage, compute initial value of gradient (gb), initialize h, the 
c  conjugate gradient direction, and compute norm squared of gradient vector.

      call prod(nx2,ny2,nz2,ns2,gx,gy,gz,u,gb)
      do 20 i=1,ns2
      h(i)=gb(i)
20    continue
c  Variable gg is the norm squared of the gradient vector
      gg=0.0
      do 105 k=1,nlist
      m=list(k)
      gg=gb(m)*gb(m)+gg
105   continue

c  Second stage, initialize Ah variable, compute parameter lamdba, 
c  make first change in voltage array, update gradient (gb) vector

      if(gg.lt.gtest) goto 44
      call prod(nx2,ny2,nz2,ns2,gx,gy,gz,h,Ah)
      hAh=0.0
      do 205 k=1,nlist
      m=list(k)
      hAh=hAh+h(m)*Ah(m)
205   continue
      lambda=gg/hAh
      do 50 i=1,ns2
      u(i)=u(i)-lambda*h(i)
      gb(i)=gb(i)-lambda*Ah(i)
50    continue

c  third stage:  iterate conjugate gradient solution process until
c  gg < gtest criterion is satisfied.  
c  (USER) The parameter ncgsteps is the total number of conjugate gradient steps
c  to go through.  Only in very unusual problems, like when the conductivity
c  of one phase is much higher than all the rest, will this many steps be 
c  used.
      ncgsteps=30000

      do 33 icc=1,ncgsteps
      gglast=gg
      gg=0.0
      do 305 k=1,nlist
      m=list(k)
      gg=gb(m)*gb(m)+gg
305   continue
      call flush(7)
      if(gg.lt.gtest) goto 44
      gamma=gg/gglast
c  update conjugate gradient direction
      do 70 i=1,ns2
      h(i)=gb(i)+gamma*h(i)
70    continue
      call prod(nx2,ny2,nz2,ns2,gx,gy,gz,h,Ah)
      hAh=0.0
      do 401 k=1,nlist
      m=list(k)
      hAh=hAh+h(m)*Ah(m)
401   continue
      lambda=gg/hAh
c  update voltage, gradient vectors
      do 90 i=1,ns2
      u(i)=u(i)-lambda*h(i)
      gb(i)=gb(i)-lambda*Ah(i)
90    continue

c  (USER) This piece of code forces dembx to write out the total current and
c  the norm of the gradient squared, every ncheck conjugate gradient steps,
c  in order to see how the relaxation is proceeding. If the currents become
c  unchanging before the relaxation is done, then gtest was picked to be
c  smaller than was necessary.  
      ncheck=30 

      if(ncheck*(icc/ncheck).eq.icc) then
      write(7,*) icc
      write(7,*) ' gg = ',gg
c  call current subroutine
      call current(nx2,ny2,nz2,ns2,currx,curry,currz,u,gx,gy,gz)
      write(7,*) ' currx = ',currx
      write(7,*) ' curry = ',curry
      write(7,*) ' currz = ',currz
      end if
      call flush(7)
33    continue
      write(7,*) ' Iteration failed to converge after',ncgsteps,' steps'
44    continue
      ic=icc

      return
      end

c The matrix product subroutine

      subroutine prod(nx2,ny2,nz2,ns2,gx,gy,gz,xw,yw)
      real gx(ns2),gy(ns2),gz(ns2),xw(ns2),yw(ns2)

c  xw is the input vector, yw = (A)(xw) is the output vector

c  auxiliary variables involving the system size
      nx1=nx2-1
      ny1=ny2-1
      nz1=nz2-1
      nx=nx2-2
      ny=ny2-2
      nz=nz2-2
      L22=nx2*ny2

c  Perform basic matrix multiplication, results in incorrect information at 
c  periodic boundaries.
      do 10 i=1,ns2
      yw(i)=0.0
10    continue
      do 20 i=L22+1,ns2-L22
      yw(i)=-xw(i)*(gx(i-1)+gx(i)+gz(i-L22)+gz(i)+gy(i)+gy(i-nx2))
      yw(i)=yw(i)+gx(i-1)*xw(i-1)+gx(i)*xw(i+1)
     + +gz(i-L22)*xw(i-L22)+gz(i)*xw(i+L22)+gy(i)*xw(i+nx2)
     + +gy(i-nx2)*xw(i-nx2)
20    continue

c  Correct terms at periodic boundaries (Section 3.3 in manual)

c  x faces
      do 30 k=1,nz2
      do 30 j=1,ny2
      yw((k-1)*L22+nx2*(j-1)+nx2)=yw((k-1)*L22+nx2*(j-1)+2)
      yw((k-1)*L22+nx2*(j-1)+1)=yw((k-1)*L22+nx2*(j-1)+nx1)
30    continue

c   y faces
      do 40 k=1,nz2
      do 40 i=1,nx2
      yw((k-1)*L22+i)=yw((k-1)*L22+ny*nx2+i)
      yw((k-1)*L22+ny1*nx2+i)=yw((k-1)*L22+nx2+i)
40    continue

c  z faces 
      do 50 m=1,L22
      yw(m)=yw(m+nz*L22)
      yw(m+nz1*L22)=yw(m+L22)
50    continue

      return
      end

c  Subroutine that determines the correct bond conductances that are used
c  to compute multiplication by the matrix A

      subroutine bond(pix,gx,gy,gz,nx2,ny2,nz2,ns2,sigma,be,nphase,ntot)
      real gx(ns2),gy(ns2),gz(ns2),sigma(ntot,3),be(ntot,ntot,3)
      integer*2 pix(ns2)

c  auxiliary variables involving the system size
      nx=nx2-2
      ny=ny2-2
      nz=nz2-2
      nx1=nx2-1
      ny1=ny2-1
      nz1=nz2-1
      L22=nx2*ny2

c  Set values of conductor for phase(i,m)--phase(j,m) interface,
c  store in array be(i,j,m), m=1,2,3. If either phase i or j 
c  has zero conductivity in the m'th direction, then be(i,j,m)= 0.0.

      do 10 m=1,3
      do 10 i=1,nphase
      do 10 j=1,nphase
      if(sigma(i,m).eq.0.0) then
      be(i,j,m)=0.0
      goto 10
      end if
      if(sigma(j,m).eq.0.0) then
      be(i,j,m)=0.0
      goto 10
      end if
      be(i,j,m)=1./(0.5/sigma(i,m)+0.5/sigma(j,m))
10    continue 

c  Trim off x and y faces so that no current can flow past periodic
c  boundaries.  This step is not really necessary, as the voltages on the 
c  periodic boundaries will be matched to the corresponding real voltages
c  in each conjugate gradient step. 
      do 20 k=1,nz2
      do 15 j=1,ny2
      gx((k-1)*L22+nx2*(j-1)+nx2)=0.0
15    continue
      do 16 i=1,nx2
      gy((k-1)*L22+ny1*nx2+i)=0.0
16    continue
20    continue

c  Set up conductor network
c   bulk--gz
      do 30 i=1,nx2
      do 30 j=1,ny2
      do 30 k=1,nz1
      m=(k-1)*L22+(j-1)*nx2+i
      i1=i
      j1=j
      k1=k+1
      m1=(k1-1)*L22+(j1-1)*nx2+i1
      gz(m)=be(pix(m),pix(m1),3)
30    continue

c  bulk---gy
      do 40 i=1,nx2
      do 40 j=1,ny1
      do 40 k=2,nz1
      m=(k-1)*L22+(j-1)*nx2+i
      j1=j+1
      i1=i
      k1=k
      m1=(k1-1)*L22+(j1-1)*nx2+i1
      gy(m)=be(pix(m),pix(m1),2)
40    continue

c  bulk--gx
      do 50 i=1,nx1
      do 50 j=1,ny2
      do 50 k=2,nz1
      m=(k-1)*L22+(j-1)*nx2+i
      i1=i+1
      j1=j
      k1=k
      m1=(k1-1)*L22+(j1-1)*nx2+i1
      gx(m)=be(pix(m),pix(m1),1)
50    continue

      return
      end

c  Subroutine that sets up the image, either by reading it from file,
c  or generating it internally

      subroutine ppixel(pix,nx2,ny2,nz2,a,ns2,nphase,ntot)
      real a(ntot)
      integer*2 pix(ns2)

c  (USER) If you want to set up a test image inside the program, instead
c  of reading it in from a file, this should be done inside this subroutine.

c  auxiliary variables involving the system size
      nx=nx2-2
      ny=ny2-2
      nz=nz2-2
      L22=nx2*ny2

c  Initialize phase fraction array.
      do 120 i=1,nphase
      a(i)=0.0
120   continue

c Use 1-d labelling scheme as shown in manual
      do 100 k=2,nz2-1
      do 100 j=2,ny2-1
      do 100 i=2,nx2-1 
      m=(k-1)*L22+(j-1)*nx2+i
      read(9,*) pix(m)
      a(pix(m))=a(pix(m))+1.0
100   continue
      do 220 i=1,nphase
      a(i)=a(i)/float(nx*ny*nz)
220   continue

c  now map periodic boundaries of pix (see Section 3.3, Figure 3 in manual)
      do 110 k=1,nz2
      do 110 j=1,ny2
      do 110 i=1,nx2
      if(i.gt.1.and.i.lt.nx2) then
         if(j.gt.1.and.j.lt.ny2) then
           if(k.gt.1.and.k.lt.nz2) then
           goto 110
           end if
         end if
      end if

      k1=k
      if(k.eq.1) k1=k+nz
      if(k.eq.nz2) k1=k-nz
      j1=j
      if(j.eq.1) j1=j+ny
      if(j.eq.ny2) j1=j-ny
      i1=i
      if(i.eq.1) i1=i+nx
      if(i.eq.nx2) i1=i-nx

      m=(k-1)*L22+(j-1)*nx2+i
      m1=(k1-1)*L22+(j1-1)*nx2+i1

      pix(m)=pix(m1)
110    continue

c  Check for wrong phase labels--less than 1 or greater than nphase
       do 500 m=1,ns2
       if(pix(m).lt.1) then
        write(7,*) 'Phase label in pix < 1--error at ',m
       end if
       if(pix(m).gt.nphase) then
        write(7,*) 'Phase label in pix > nphase--error at ',m
       end if
500    continue

      return
      end

c  Subroutine to compute the total current in the x, y, and z directions

      subroutine current(nx2,ny2,nz2,ns2,currx,curry,currz,u,gx,gy,gz)
      real u(ns2),gx(ns2),gy(ns2),gz(ns2),currx,curry,currz

c  auxiliary variables involving the system size
      nx=nx2-2
      ny=ny2-2
      nz=nz2-2
      L22=nx2*ny2
c  initialize the volume averaged currents
      currx=0.0
      curry=0.0
      currz=0.0

c  Only loop over real sites and bonds in order to get true total current
      do 10 k=2,nz2-1 
      do 10 j=2,ny2-1 
      do 10 i=2,nx2-1 
      m=L22*(k-1)+nx2*(j-1)+i
c  cur1, cur2, cur3 are the currents in one pixel
      cur1=0.5*( ( u(m-1) - u(m) )*gx(m-1)+( u(m)-u(m+1) )*gx(m) )
      cur2=0.5*( ( u(m-nx2) - u(m) )*gy(m-nx2)+( u(m)-u(m+nx2) )*gy(m) )
      cur3=0.5*( ( u(m-L22) - u(m) )*gz(m-L22)+( u(m)-u(m+L22) )*gz(m) )
c  sum pixel currents into volume averages
      currx=currx+cur1
      curry=curry+cur2
      currz=currz+cur3
10    continue

      currx=currx/float(nx*ny*nz)
      curry=curry/float(nx*ny*nz)
      currz=currz/float(nx*ny*nz)

      return
      end
