program main
implicit none
!----------------
integer,parameter::nx=101,ny=51
real*8::p=0.3,q=0.3
real*8::dks=1.0,dat=1.0
real*8,dimension(nx,ny)::x=0,y=0,xp=0,yp=0
real*8,allocatable::xb0(:),yb0(:),xb1(:),yb1(:)
!---------------
integer::nb0,nb1
integer::i,ii,j,times,di,dj
!----------------
real*8::errmax=0.5,err
real*8::xks,xat,yks,yat,xksat,yksat,arfa,beta,gama,jj
real*8::a,ap,b,bp,c,cp,d,dp,e,ep
real*8::l,lp,m,mp,n,np
real*8::xe,xw,xn,xs,ye,yw,yn,ys
real*8::xx,yy
!---------------
real*8::thta,thta1,thta2,k1,k2,dis,dis1,disa,disb,disc
real*8::res
integer::iffind
!------------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------------
!边界坐标
open(1,file='boundary.txt')
read(1,*) nb0 !左岸边界
allocate(xb0(nb0),yb0(nb0))
read(1,*) ((xb0(i),yb0(i)),i=1,nb0)
read(1,*) nb1 !右岸岸边界
allocate(xb1(nb1),yb1(nb1))
read(1,*) ((xb1(i),yb1(i)),i=1,nb1)
close(1)
!initial----------------------------------------------------------------------------------------
call cboundary(nx,nb0,xb0,yb0,x(:,1),y(:,1)) !(1,1) - (nx,1) 河床左岸边界坐标
call cboundary(nx,nb1,xb1,yb1,x(:,ny),y(:,ny)) !(1,ny) - (nx,ny) 河床右岸边界坐标
!------------------------------------
do i=1,nx,nx-1 !(1,1) - (1,ny),(nx,1) - (nx,ny) 河床上下边界坐标
do j=2,ny-1
x(i,j)=x(i,1)+(x(i,ny)-x(i,1))*(j-1)/(ny-1)
y(i,j)=y(i,1)+(y(i,ny)-y(i,1))*(j-1)/(ny-1)
enddo
enddo
do j=2,ny-1 !(x,y) 河床内部节点边界
do i=2,nx-1
x(i,j)=x(1,j)+(x(nx,j)-x(1,j))*(i-1)/(nx-1)
y(i,j)=y(1,j)+(y(nx,j)-y(1,j))*(i-1)/(nx-1)
enddo
enddo
do i=2,nx-1
do j=2,ny-1
x(i,j)=x(i,1)+(x(i,ny)-x(i,1))*(j-1)/(ny-1)
y(i,j)=y(i,1)+(y(i,ny)-y(i,1))*(j-1)/(ny-1)
enddo
enddo
!--------
open(1,file='initial.dat')
write(1,*) '"variable=x,y"'
write(1,*) 'zone i=',ny, 'j=',nx
do i=1,nx
do j=1,ny
write(1,*) x(i,j),y(i,j)
enddo
enddo
close(1)
!-------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
times=0
100 continue
err=0
do i=2,nx-1
do j=2,ny-1
xks=(x(i+1,j)-x(i-1,j))/(2*dks)
xat=(x(i,j+1)-x(i,j-1))/(2*dat)
yks=(y(i+1,j)-y(i-1,j))/(2*dks)
yat=(y(i,j+1)-y(i,j-1))/(2*dat)
xksat=(x(i+1,j+1)+x(i-1,j-1)-x(i+1,j-1)-x(i-1,j+1))/(4*dks*dat)
yksat=(y(i+1,j+1)+y(i-1,j-1)-y(i+1,j-1)-y(i-1,j+1))/(4*dks*dat)
!----------
arfa=xat*xat+yat*yat
gama=xks*xks+yks*yks
beta=xksat+yksat
jj=xks*yat-xat*yks
! 考虑P、Q-----------------------------------------------------------------------
! a=arfa*yks*yks/(xks*xks+yks*yks)
! b=gama*yat*yat/(xat*xat+yat*yat)
! c=-arfa*xks*yks/(xks*xks+yks*yks)
! d=-gama*xat*yat/(xat*xat+yat*yat)
! e=-2*beta*xksat
! ap=c
! bp=d
! cp=arfa*xks*xks/(xks*xks+yks*yks)
! dp=gama*xat*
xat/(xat*xat+yat*yat)
! ep=-2*beta*yksat
! !--------------------------------
! xe=x(i+1,j)
! xw=x(i-1,j)
! xn=x(i,j+1)
! xs=x(i,j-1)
! ye=y(i+1,j)
! yw=y(i-1,j)
! yn=y(i,j+1)
! ys=y(i,j-1)
! !--------------------------------
! l=a*(xe+xw)/dks/dks+b*(xn+xs)/dat/dat+c*(ye+yw)/dks/dks+d*(yn+ys)/dat/dat+e
! m=2*(a+b)
! n=2*(c+d)
! lp=ap*(xe+xw)/dks/dks+bp*(xn+xs)/dat/dat+cp*(ye+yw)/dks/dks+dp*(yn+ys)/dat/dat+ep
! mp=2*(ap+bp)
! np=2*(cp+dp)
! !----------------------------
! xp(i,j)=(l*np-lp*n)/(m*np-mp*n)
! yp(i,j)=(l*mp-lp*m)/(mp*n-m*np)
! if(abs(xp(i,j)-x(i,j))>err) err=abs(xp(i,j)-x(i,j))
! if(abs(yp(i,j)-y(i,j))>err) err=abs(yp(i,j)-y(i,j))
! x(i,j)=xp(i,j)
! y(i,j)=yp(i,j)
! 不考虑P，Q----------------------------------------------------------------------
p=0.
q=0.
!---------
xp(i,j)=arfa*(x(i-1,j)+x(i+1,j))/dks/dks+gama*(x(i,j+1)+x(i,j-1))/dat/dat-2*beta*xksat+jj*(p*xks+q*xat)
xp(i,j)=xp(i,j)/(2.0*(arfa/dks/dks+gama/dat/dat))
if(abs(xp(i,j)-x(i,j))>err) err=abs(xp(i,j)-x(i,j))
x(i,j)=xp(i,j)
!---------
yp(i,j)=arfa*(y(i-1,j)+y(i+1,j))/dks/dks+gama*(y(i,j+1)+y(i,j-1))/dat/dat-2*beta*yksat+jj*(p*yks+q*yat)
yp(i,j)=yp(i,j)/(2.0*(arfa/dks/dks+gama/dat/dat))
if(abs(yp(i,j)-y(i,j))>err) err=abs(yp(i,j)-y(i,j))
y(i,j)=yp(i,j)
enddo
enddo
!--------------------------------------------------------------------------------
goto 99
!边界滑移------------------------------------------------------------------------
!x boundary----
!do j=1,ny,ny-1
! if(j==1) dj=1 !bottom
! if(j==ny) dj=-1 !top
! do i=2,nx-1
! !print*,i,j
! !寻找边界控制点坐标
! xx=0
! yy=0
! iffind=0
! if(j==1)then
! do ii=1,nb0-1
! if(x(i,j)>=min(xb0(ii),xb0(ii+1)).and.x(i,j)
! yy=yb0(ii)
! iffind=1
! cycle
! endif
! enddo
! else
! do ii=1,nb1-1
! if(x(i,j)>=xb1(ii).and.x(i,j)
! yy=yb1(ii)
! iffind=1
! cycle
! endif
! enddo
! endif
! if(iffind/=1) then
! print*,'找不到边界控制点坐标'
! stop
! endif
! !-----------------------------------------------------------------
! ! x(i,j+dj)
! ! / *
! ! / / |
! ! c / / |
! ! / a / |
! ! / / |
! ! / b /thta |
! ! ---------*---------------*--------------*----------
! ! (xx,yy) x(i,j) (xx1,yy2)
! !------------------------------------------------------------------
! disa=sqrt((x(i,j+dj)-x(i,j))**2.0+(y(i,j+dj)-y(i,j))**2.0)
! disb=sqrt((xx-x(i,j))**2.0+(yy-y(i,j))**2.0)
! disc=sqrt((x(i,j+dj)-xx)**2.0+(y(i,j+dj)-yy)**2.0)
! if(disb/=0.and.disa/=0)then
! thta=(disa*disa+disb*disb-disc*disc)/(2*disa*disb)
! if(abs(thta)>1) thta=thta/abs(thta)
! thta=acos(thta)
! thta=3.1415926-thta
! x(i,j)=x(i,j)+(
x(i,j)-xx)*disa*cos(thta)/disb
! y(i,j)=y(i,j)+(y(i,j)-yy)*disa*cos(thta)/disb
! endif
!! disb=sqrt((x(i+di,j)-x(i,j))**2.0+(y(i+di,j)-y(i,j))**2.0)
!! disc=sqrt((x(i,j+dj)-x(i+di,j))**2.0+(y(i,j+dj)-y(i+di,j))**2.0)
!! thta=acos((disa*disa+disb*disb-disc*disc)/(2*disa*disb))
!! x(i,j)=x(i,j)+(x(i+di,j)-x(i,j))*disa*cos(thta)/disb
!! y(i,j)=y(i,j)+(y(i+di,j)-y(i,j))*disa*cos(thta)/disb
! enddo
!enddo
!y boundary
!do i=1,nx,nx-1
! if(i==1) di=1
! if(i==nx) di=-1
! do j=2,ny-1
! disa=sqrt((x(i+di,j)-x(i,j))**2.0+(y(i+di,j)-y(i,j))**2.0)
! if(j
! elseif(j==ny)then
! dj=-1
! endif
! disb=sqrt((x(i,j+dj)-x(i,j))**2.0+(y(i,j+dj)-y(i,j))**2.0)
! disc=sqrt((x(i+di,j)-x(i,j+dj))**2.0+(y(i+di,j)-y(i,j+dj))**2.0)
! thta=acos((disa*disa+disb*disb-disc*disc)/(2*disa*disb))
! x(i,j)=x(i,j)+(x(i,j+dj)-x(i,j))*disa*cos(thta)/disb
! y(i,j)=y(i,j)+(y(i,j+dj)-y(i,j))*disa*cos(thta)/disb
! enddo
!enddo
!----------------------------------------------------
99 continue
!----------------------------------------------------
if(err>errmax)then
times=times+1
if(mod(times,500)==0)then
print*,'times=',times,' err=',err
!if(times==10)then
! goto 101
!endif
endif
goto 100
endif
!---------------------------------------------------------------------------
101 continue
!---------------------------------------------------------------------------
open(1,file='res.dat')
write(1,*) '"variable=x,y"'
write(1,*) 'zone i=',ny, 'j=',nx
do i=1,nx
do j=1,ny
write(1,*) x(i,j),y(i,j)
enddo
enddo
close(1)
!----------------------------------------------------------------------------
print*,'-------finished-------'
print*,'err=',err,' times=',times
!----------------------------------------------------------------------------
end
subroutine cboundary(nx,nb,xb,yb,x,y)
implicit none
!--------------------------------------------------
integer::nx,nb
real*8,dimension(nb)::xb,yb
real*8,dimension(nx)::x,y
!--------------------------------------------------
real*8::dis,dis1,disa,disb,disc
integer::i,ii
!-----------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------
dis=0
do i=2,nb
dis=dis+sqrt((xb(i)-xb(i-1))**2.0+(yb(i)-yb(i-1))**2.0)
enddo
!---------------------------------------------------
x(1)=xb(1)
y(1)=yb(1)
x(nx)=xb(nb)
y(nx)=yb(nb)
dis1=dis/(nx-1)
disa=0
!-------------------
do i=2,nx-1
disa=(i-1)*dis1
disb=sqrt((xb(1)-xb(2))**2.0+(yb(1)-yb(2))**2.0)
if(disa<=disb)then
x(i)=xb(1)+(xb(2)-xb(1))*disa/disb
y(i)=yb(1)+(yb(2)-yb(1))*disa/disb
else
disb=0
do ii=2,nb-1
disb=disb+sqrt((xb(ii)-xb(ii-1))**2.0+(yb(ii)-yb(ii-1))**2.0)
disc=disb+sqrt((xb(ii+1)-xb(ii))**2.0+(yb(ii+1)-yb(ii))**2.0)
if(disa>disb.and.disa<=disc)then
x(i)=xb(ii)+(xb(ii+1)-xb(ii))*(disa-disb)/(disc-disb)
y(i)=yb
(ii)+(yb(ii+1)-yb(ii))*(disa-disb)/(disc-disb)
exit
endif
enddo
if(disa>disc)then
print*,'cannot find position'
print*,i,disa,disb,disc
stop
endif
endif
enddo
!----------------------------------------------------------------------------------------------
return
end subroutine cboundary