c-----------------------------------------------------------------------
c grid.f
c
c Grid-to-Grid operations for AMR
c - Prolongation operators
c - Restriction operators
c-----------------------------------------------------------------------
c Routine to convert from global array coordinates to local array coordinates
integer function getindx(loc, lb, stride )
implicit none
integer loc, lb, stride
getindx = (loc - lb)/stride + 1
return
end
c-----------------------------------------------------------------------
c Prolongation operator: prolong2d2 (2D 2nd order)
c USES INTEGER Lower and Upper bounds for use with DAGH
c
c - Uses linear interpolation where necessary
c Interface:
c shapef(2) := shape of fine grid function
c shapec(2) := shape of coarse grid function
c
c uf(,) := fine grid function
c uc(,) := coarse grid function
c
c lbc(2) := lower bound for coarse grid
c ubc(2) := upper bound for coarse grid
c lbf(2) := lower bound for fine grid
c ubf(2) := upper bound for fine grid
c lbr(2) := lower bound for region prolongation desired
c ufr(2) := upper bound for region prolongation desired
subroutine prolong2d2(
. uc, lbc, ubc, shapec,
. uf, lbf, ubf, shapef,
. lbr, ubr, args, argc)
implicit none
integer shapef(2), shapec(2)
real*8 uf(shapef(1), shapef(2)),
. uc(shapec(1), shapec(2))
integer lbf(2), ubf(2),
. lbc(2), ubc(2),
. lbr(2), ubr(2),
. getindx
integer argc
real*8 args(argc)
c Local variables
integer i, j, ic,jc,
. ilbf(2), iubf(2),
. ilbc(2), iubc(2),
. refine,
. ifine, icoarse,
. jfine, jcoarse, stridec, stridef
real*8 hf, hc, rem(2), x, y,
. a11, a12, a21, a22
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Get strides coarse grid
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
stridec = (ubc(1) - lbc(1))/(shapec(1)-1)
stridef = (ubf(1) - lbf(1))/(shapef(1)-1)
refine = stridec/stridef
c write(45,*)'In prolong: stridec=',stridec,' stridef=',stridef
c write(46,*)'In prolong: stridec=',
c . (ubc(1) - lbc(1))/(shapec(1)-1),
c . ' stridef=',
c . (ubf(1) - lbf(1))/(shapef(1)-1)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Prolongation region is defined on the domain of the fine grid.
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Loop through points in prolongation region
c write(55,*)shapec(1), shapec(2)
c write(55,*)shapef(1), shapef(2)
c write(55,*)'stride coarse ',stridec
c write(55,*)'stride fine ',stridef
c write(55,*)'Region l/u x l/u y-----------------------'
c write(55,*)lbr(1), ubr(1),lbr(2), ubr(2)
c write(55,*)'Coarse l/u x l/u y-----------------------'
c write(55,*)lbc(1), ubc(1),lbc(2), ubc(2)
c write(55,*)'Fine l/u x l/u y-----------------------'
c write(55,*)lbf(1), ubf(1),lbf(2), ubf(2)
c write(55,*)'********************'
do i = lbr(1), ubr(1), stridef
do j = lbr(2), ubr(2), stridef
c calculate coarse grid integer coordinates difference
ic = i - lbc(1)
jc = j - lbc(2)
ifine = getindx(i, lbf(1), stridef)
jfine = getindx(j, lbf(2), stridef)
icoarse = getindx(i, lbc(1), stridec)
jcoarse = getindx(j, lbc(2), stridec)
c if(trace)then
c write(6,*)icoarse, jcoarse, ifine, jfine
c write(6,*)ic, jc, mod(ic,refine), mod(jc,refine)
c write(6,*)'**********'
c end if
c if(trace)then
c write(6,*)i,j,ic,jc
c write(6,*)'m ',mod(i-1,refine),mod(j-1,refine)
c write(6,*)'********************************** '
c end if
if(mod(ic,stridec) .eq. 0 .and. mod(jc,stridec) .eq.0)then
uf(ifine,jfine) = uc(icoarse,jcoarse)
end if
if(mod(ic,stridec) .eq. 0 .and. mod(jc,stridec) .ne.0)then
uf(ifine,jfine) = 0.5 * uc(icoarse,jcoarse) +
. 0.5 * uc(icoarse,jcoarse+1)
end if
if(mod(ic,stridec) .ne. 0 .and. mod(jc,stridec) .eq.0)then
uf(ifine,jfine) = 0.5 * uc(icoarse,jcoarse) +
. 0.5 * uc(icoarse+1,jcoarse)
end if
if(mod(ic,stridec) .ne. 0 .and. mod(jc,stridec).ne. 0)then
uf(ifine,jfine) = 0.25 * uc(icoarse,jcoarse) +
. 0.25 * uc(icoarse+1,jcoarse) +
. 0.25 * uc(icoarse,jcoarse+1) +
. 0.25 * uc(icoarse+1,jcoarse+1)
end if
end do
end do
return
end
c-----------------------------------------------------------------------
c Restriction operator: restrict2d (pure injections)
c Interface:
c shapef(2) := shape of fine grid function
c shapec(2) := shape of coarse grid function
c
c uf(,) := fine grid function
c uc(,) := coarse grid function
c
c lbc(2) := lower bound for coarse grid
c ubc(2) := upper bound for coarse grid
c lbf(2) := lower bound for fine grid
c ubf(2) := upper bound for fine grid
c lbr(2) := lower bound for region restriction desired
c ufr(2) := upper bound for region restriction desired
subroutine restrict2d2(uf, lbf, ubf, shapef,
. uc, lbc, ubc, shapec,
. lbr, ubr, args, argc)
implicit none
integer shapef(2), shapec(2)
real*8 uf(shapef(1), shapef(2)),
. uc(shapec(1), shapec(2))
integer lbf(2), ubf(2),
. lbc(2), ubc(2),
. lbr(2), ubr(2)
integer argc
real*8 args(argc)
c Local variables
integer i, j, imin, imax, jmin, jmax,
. ifine, icoarse,
. jfine, jcoarse, refine, stridec, stridef,
. ilbc(2), iubc(2), getindx
real*8 hf, hc, x, y
stridec = (ubc(1) - lbc(1))/(shapec(1)-1)
stridef = (ubf(1) - lbf(1))/(shapef(1)-1)
refine = stridec/stridef
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Find coarse domain over which to refine
c Take three regions and select out intersection
imin = max(lbf(1), lbc(1), lbr(1))
imax = min(ubf(1), ubc(1), ubr(1))
jmin = max(lbf(2), lbc(2), lbr(2))
jmax = min(ubf(2), ubc(2), ubr(2))
if (mod(imin-lbc(1),stridec) .ne. 0) then
imin = imin + mod(imin-lbc(1),stridec)
endif
if (mod(jmin-lbc(2),stridec) .ne. 0) then
jmin = jmin + mod(jmin-lbc(2),stridec)
endif
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...
do i = imin, imax,stridec
do j = jmin, jmax,stridec
ifine = getindx(i, lbf(1), stridef)
jfine = getindx(j, lbf(2), stridef)
icoarse = getindx(i, lbc(1), stridec)
jcoarse = getindx(j, lbc(2), stridec)
if(icoarse .gt. shapec(1) .or. jcoarse .gt. shapec(2))then
write(0,*)'ERROR in restriction: ',icoarse,jcoarse
end if
uc(icoarse,jcoarse) = uf(ifine,jfine)
end do
end do
return
end
c-----------------------------------------------------------------------