A algorithm problem meets in the FTC3D code
 

Hi all,
Recently,I download the FTC(Front-Tracking code)3D code from the internet. But the FTC3D had been modified already. Specially, I think the subroutine about updating the density field is not exact. Then I found the FTC2D code from Dr. Tryggvason's personal web:http://www3.nd.edu/~gtryggva/. And I rewrite the subroutine of the FTC3D record to the algorithm of the FTC2D, but I meet a problem!
The subroutine in the FTC2D is:
subroutine dens2(r,rtmp,tmp1,tmp2,nxp2,nyp2,iii,idim,pt,ptcon ,icp,elcon,
elprop,maxpt,maxel)

dimension tmp1(nxp2+2,nyp2+2), tmp2(nxp2+2,nyp2+2)
dimension rtmp(nxp2,nyp2)
dimension r(nxp2,nyp2)
dimension iii(idim,2)

dimension pt(maxpt,2), ptcon(maxpt)
dimension icp(maxel,2),elcon(maxel),elprop(maxel,3)

integer ptcon,elcon,np,ffp,lfp,fep, ne,ffe,lfe,fee

common/grid/hxi,hyi,nxp1,nyp1,nx,ny,xl,yl
common/bounds/ibdry,jbdry,nxlast,nylast
common/front/np,ffp,lfp,fep, ne,ffe,lfe,fee,amax,amin
common/fparam/xmv,mv,fxl,fyl
common/flprop/r1,r2,m1,m2,gx,gy,rro

xx=(xmv+0.5)*xl
yy=(xmv+0.5)*yl

c generate gradient vector

sc=1./2.
pi=4.*ATAN(1.)
pd2=0.5*pi
cnstx=(.25*hxi)**2
cnsty=(.25*hyi)**2

do 10 i=1,nxp2+2
do 10 j=1,nyp2+2
tmp1(i,j)=0.0
tmp2(i,j)=0.0
10 continue

k=ffe

do 30 ll=1,ne

p0x=pt(icp(k,1),1)
p0y=pt(icp(k,1),2)

p1x= p0x+amod((pt(icp(k,2),1)-p0x+xx),xl)-0.5*xl
p1y= p0y+amod((pt(icp(k,2),2)-p0y+yy),yl)-0.5*yl

xc=0.5*(p1x+p0x)
yc=0.5*(p1y+p0y)

x1=p1x-p0x
y1=p1y-p0y

c calculate components of n.dA

rx=-y1 *elprop(k,1)
ry= x1 *elprop(k,1)

c generate a grid value

xt=amod((xc+0.5/hxi)+xmv*xl,xl)
yt=amod((yc+0.5/hyi)+xmv*yl,yl)

ir=1+int(nx*xt/xl)
jr=1+int(ny*yt/yl)


do 40 i1=1,4
do 40 j1=1,4

ii=ir-2+i1
jj=jr-2+j1

drx = 1. + cos((xt*hxi - float(ii-1))*pd2)
dry = 1. + cos((yt*hyi - float(jj-1))*pd2)

tmp1(ii+1,jj+1)=tmp1(ii+1,jj+1)+rx*drx*dry*cnstx
tmp2(ii+1,jj+1)=tmp2(ii+1,jj+1)+ry*drx*dry*cnsty

40 continue

k=elcon(k)
30 continue

c correct boundaries:

call gbd1(tmp1,tmp2,ibdry,jbdry,nxp2,nyp2)

do 900 i=1,nxp2
do 900 j=1,nyp2
900 rtmp(i,j)=0.0

do 80 i=2,nxp1
do 80 j=2,nyp1
rtmp(i,j)= 0.5*(hxi*(tmp1(i+1+1,j+1)-tmp1(i-1+1,j+1))+
1 hyi*(tmp2(i+1,j+1+1)-tmp2(i+1,j-1+1)))
80 continue

c set connection
node=0
do 100 i=2,nxp1
do 100 j=2,nyp1
cc=0.0
do 105 io=1,3
do 105 jo=1,3
105 cc=cc+rtmp(i+io-2,j+jo-2)
if(abs(cc) .gt. 1.0e-03)then
node=node+1
iii(node,1)=i
iii(node,2)=j
end if
c open(13,file='cc',status='new')
c write(13,*) cc,i,j,node
100 continue
c pause

numnd=node

c then iterate

c-----------the following should be read in
itmax=1000
bet=1.2
errmax=1.0e-05
c------------------------------------------
cf=0.5/(hxi*hxi+hyi*hyi)
do 200 iter=1,itmax

err=0.0
do 210 node=1,numnd
i=iii(node,1)
j=iii(node,2)
c do 210 i=2,nxp1
c do 210 j=2,nyp1
rold=r(i,j)
r(i,j)=(1.0-bet)*r(i,j)+bet*cf*
1 (hxi*hxi*(r(i+1,j)+r(i-1,j))+
2 hyi*hyi*(r(i,j+1)+r(i,j-1))-rtmp(i,j))
err=err+abs(rold-r(i,j))
210 continue


c set boundary terms...
if(ibdry .eq. 0)then
do 110 j=1,nyp2
r(nxp2,j)=r(2,j)
r(1,j)=r(nxp1,j)
110 continue
else
do 115 j=1,nyp2
r(nxp2,j)=r(nxp1,j)
r(1,j)=r(2,j)
115 continue
end if

if(jbdry .eq. 0)then
do 120 i=1,nxp2
r(i,nyp2)=r(i,2)
r(i,1)=r(i,nyp1)
120 continue
else
do 125 i=1,nxp2
r(i,nyp2)=r(i,nyp1)
r(i,1)=r(i,2)
125 continue
end if

c check convergence......

if(err.lt.(errmax*float(numnd)))go to 300
200 continue
300 continue



write(*,*)' DENSITY: ',err,' after ',iter,' iterations',
1 ' using ',numnd,' nodes'

return
end

And I rewrite the 3D subroutine:
subroutine creater(r,rtmp,tmp1,tmp2,tmp3,nxp2,nyp2,nzp2,iii,i dim,pt,ptcon,
icp,elcon,maxpt,maxel)

dimension tmp1(nxp2,nyp2,nzp2),tmp2(nxp2,nyp2,nzp2)
dimension tmp3(nxp2,nyp2,nzp2),rtmp(nxp2,nyp2,nzp2)
dimension r(nxp2,nyp2,nzp2)
dimension iii(idim,3)

dimension pt(maxpt,3), ptcon(maxpt)
dimension icp(maxel,3),elcon(maxel)

dimension work(10000),bd(1,1)

integer ptcon,elcon,np,ffp,lfp,fep, ne,ffe,lfe,fee
integer p1,p2,p3

common/grid/hxi,hyi,hzi,nxp1,nyp1,nzp1,nx,ny,nz,xl,yl,zl
common/front/np,ffp,lfp,fep, ne,ffe,lfe,fee,amax,amin,aspmax
common/flprop/r1,r2,m1,m2,gx,gy,gz,st,rro

c generate gradient vector

sc=1./3.
pd2=0.5*pimach(1.0)

cnstx=(.25*hxi)**3
cnsty=(.25*hyi)**3
cnstz=(.25*hzi)**3

dxh=0.5/hxi
dyh=0.5/hyi
dzh=0.5/hzi


do 10 i=1,nxp2
do 10 j=1,nyp2
do 10 k=1,nzp2

tmp1(i,j,k)=0.0
tmp2(i,j,k)=0.0
tmp3(i,j,k)=0.0


10 continue

k=ffe


do 30 ll=1,ne

p1=icp(k,1)
p2=icp(k,2)
p3=icp(k,3)

xc=(pt(p1,1)+pt(p2,1)+pt(p3,1))*sc
yc=(pt(p1,2)+pt(p2,2)+pt(p3,2))*sc
zc=(pt(p1,3)+pt(p2,3)+pt(p3,3))*sc

x1=pt(p2,1)-pt(p1,1)
y1=pt(p2,2)-pt(p1,2)
z1=pt(p2,3)-pt(p1,3)

x2=pt(p3,1)-pt(p1,1)
y2=pt(p3,2)-pt(p1,2)
z2=pt(p3,3)-pt(p1,3)


c calculate components of n.dA

rx=(y1*z2-y2*z1)
ry=(z1*x2-z2*x1)
rz=(x1*y2-x2*y1)

c generate a grid value

xt=amod((xc+dxh)+5.0*xl,xl)
yt=amod((yc+dyh)+5.0*yl,yl)
zt=amod((zc+dzh)+5.0*zl,zl)

ir=1+int(nx*xt/xl)
jr=1+int(ny*yt/yl)
kr=1+int(nz*zt/zl)

do 40 i1=1,4
do 40 j1=1,4
do 40 k1=1,4

ii=ir-2+i1
jj=jr-2+j1
kk=kr-2+k1


drx = 1. + cos((xt*hxi - float(ii-1))*pd2)
dry = 1. + cos((yt*hyi - float(jj-1))*pd2)
drz = 1. + cos((zt*hzi - float(kk-1))*pd2)

ii=mod((ii-1+5*nx),nx)+1
jj=mod((jj-1+5*ny),ny)+1
kk=mod((kk-1+5*nz),nz)+1

tmp1(ii,jj,kk)=tmp1(ii,jj,kk)+rx*drx*dry*drz*cnstx
tmp2(ii,jj,kk)=tmp2(ii,jj,kk)+ry*drx*dry*drz*cnsty
tmp3(ii,jj,kk)=tmp3(ii,jj,kk)+rz*drx*dry*drz*cnstz


40 continue

k=elcon(k)

30 continue

c correct the ends of the box

call gbdry2(tmp1,nxp2,nyp2,nzp2)
call gbdry2(tmp2,nxp2,nyp2,nzp2)
call gbdry2(tmp3,nxp2,nyp2,nzp2)

do 900 i=1,nxp2
do 900 j=1,nyp2
do 900 k=1,nzp2
rtmp(i,j,k)=0.0
900 continue

do 80 i=2,nxp1
do 80 j=2,nyp1
do 80 k=2,nzp1

rtmp(i,j,k)=0.5*(hxi*hxi*(tmp1(i+1,j,k)-tmp1(i-1,j,k))+
1 hyi*hyi*(tmp2(i,j+1,k)-tmp2(i,j-1,k))+
2 hzi*hzi*(tmp3(i,j,k+1)-tmp3(i,j,k-1)))


80 continue


node=0

do 100 i=2,nxp1
do 100 j=2,nyp1
do 100 k=2,nzp1

cc=0.0

do 105 io=1,3
do 105 jo=1,3
do 105 ko=1,3

cc=cc+rtmp(i+io-2,j+jo-2,k+ko-2)


105 continue

if(abs(cc) .gt. 1.0e-3)then

node=node+1
iii(node,1)=i
iii(node,2)=j
iii(node,3)=k

end if

100 continue


numnd=node



itmax=1000
errmax=1.0e-3
beta=1.2

cf1=0.5/(hxi*hxi*hxi+hyi*hyi*hyi+hzi*hzi*hzi)


do 200 iter=1,itmax

err1=0.0



do 210 node=1,numnd


i=iii(node,1)
j=iii(node,2)
k=iii(node,3)

rold=r(i,j,k)

r(i,j,k)=(1.0-beta)*r(i,j,k)+beta*cf1*
1 (hxi*hxi*hxi*(r(i+1,j,k)+r(i-1,j,k))+
2 hyi*hyi*hyi*(r(i,j+1,k)+r(i,j-1,k))+
3 hzi*hzi*hzi*(r(i,j,k+1)+r(i,j,k-1))-rtmp(i,j,k))

err1=err1+abs(rold-r(i,j,k))

210 continue

call gbdry(r,nxp2,nyp2,nzp2)


if(err1.lt.errmax) go to 300

200 continue
300 continue


write(*,*)' DENSITY: ',err1,' after ',iter,' iterations',
1 ' using ',numnd,' nodes'


return
end

The results of the 2D and 3D code are in the attachment.
I readed the article of Dr.Tryggvason "A Front-Tracking Method for the Computations of Multiphase Flow'', the section 3.5 "Updating the Material Properties" explain the algorithm of Updating the density field. But I can't understand. So I need some help to modify the code. Thanks!
Yours'
Liaobin