3d/operators/restrict3.f

c-----------------------------------------------------------------------
c     Three-dimensional restriction operator for AMROC.
c     A coarse cell value is overwritten by the mean value 
c     of all refined cells within this particular coarse cell.
c
c     Interface:
c        mfx,mfy,mfz := shape of fine grid
c        mcx,mcy,mcz := shape of coarse grid
c
c        uf(,,) := fine grid
c        uc(,,) := coarse grid
c
c        lbc(3) := lower bound for coarse grid
c        ubc(3) := upper bound for coarse grid
c        lbf(3) := lower bound for fine grid
c        ubf(3) := upper bound for fine grid
c        lbr(3) := lower bound for region restriction desired
c        ufr(3) := upper bound for region restriction desired
c        shaper(3) := shape of region restriction desired
c
c     Author: Ralf Deiterding
c-----------------------------------------------------------------------

      subroutine restrict3(uf,mfx,mfy,mfz,lbf,ubf, 
     &     uc,mcx,mcy,mcz,lbc,ubc,
     &     lbr,ubr,shaper,args,argc)
      implicit none
      common /ghosts/ mbc
      integer mbc,mfx,mfy,mfz,mcx,mcy,mcz
      
      integer shaper(3)
      integer args, argc            
      real*8  uf(args,mfx,mfy,mfz), uc(args,mcx,mcy,mcz)
      
      integer  lbf(3), ubf(3), 
     &     lbc(3), ubc(3), 
     &     lbr(3), ubr(3)
      
c      Local variables

      integer i, j, k, m, imin, imax, jmin, jmax, kmin, kmax, 
     &     ii, jj, kk, ifine, icoarse, jfine, jcoarse, kfine, kcoarse, 
     &     refine, stridec, stridef, meqn,
     &     getindx
      integer mbcf

      meqn = args
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c      See definition of member-function extents() in BBox.h 
c      for calculation of stride
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         
      stridec = (ubc(1) - lbc(1))/(mcx-1)
      stridef = (ubf(1) - lbf(1))/(mfx-1)
      refine = stridec/stridef
      mbcf = mbc * stridef

c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c     Find coarse domain over which to refine
c     Take three regions and select out intersection
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      imin = max(lbf(1)+mbcf, lbc(1), lbr(1))
      imax = min(ubf(1)-mbcf, ubc(1), ubr(1))
      jmin = max(lbf(2)+mbcf, lbc(2), lbr(2))
      jmax = min(ubf(2)-mbcf, ubc(2), ubr(2))
      kmin = max(lbf(3)+mbcf, lbc(3), lbr(3))
      kmax = min(ubf(3)-mbcf, ubc(3), ubr(3))
      
      if (mod(imin-lbc(1),stridec) .ne. 0) then
         imin = imin + stridec - mod(imin-lbc(1),stridec) 
      endif
      if (mod(jmin-lbc(2),stridec) .ne. 0) then
         jmin = jmin + stridec - mod(jmin-lbc(2),stridec) 
      endif
      if (mod(kmin-lbc(3),stridec) .ne. 0) then
         kmin = kmin + stridec - mod(kmin-lbc(3),stridec) 
      endif

c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c     Inject points to coarse grid from fine grid
c     Loop from lower bound to upper bound with stride of refine.
c     Convert the integer coordinates to fine and coarse grid absolute
c     coordinates...
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

      do 10 k=kmin, kmax, stridec
         kfine = getindx(k, lbf(3), stridef)
         kcoarse = getindx(k, lbc(3), stridec)
         
         do 10 j=jmin, jmax, stridec
            jfine = getindx(j, lbf(2), stridef)
            jcoarse = getindx(j, lbc(2), stridec)
            
            do 10 i=imin, imax, stridec            
               ifine = getindx(i, lbf(1), stridef)               
               icoarse = getindx(i, lbc(1), stridec)
               
!               if (icoarse .gt. mcx .or. 
!     &              jcoarse .gt. mcy .or. 
!     &              kcoarse .gt. mcz) then
!                  write(0,*)'ERROR in restriction: ',
!     &                 icoarse,jcoarse,kcoarse
!               end if

               do 10 m=1, meqn
                  uc(m,icoarse,jcoarse,kcoarse) = 0
                  do 20 kk=0, refine-1
                     do 20 jj=0, refine-1
                        do 20 ii=0, refine-1
                           uc(m,icoarse,jcoarse,kcoarse) = 
     &                          uc(m,icoarse,jcoarse,kcoarse) + 
     &                          uf(m,ifine+ii,jfine+jj,kfine+kk)
 20               continue
                  uc(m,icoarse,jcoarse,kcoarse) = 
     &                 uc(m,icoarse,jcoarse,kcoarse) / refine**3
                        
 10   continue

      return
      end