Fortran正交曲线网格生成程序.zip_PI9_正交曲线网格_网格fortran_网格生成_网格生成 正交

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)<max(xb0(ii+1),xb0(ii)))then !xx=xb0(ii) !yy=yb0(ii) !iffind=1 !cycle !endif !enddo !else !do ii=1,nb1-1 !if(x(i,j)>=xb1(ii).and.x(i,j)<xb1(ii+1))then !xx=xb1(ii) !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<ny)then !dj=1 !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