|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
[这个贴子最后由harry在 2004/07/15 11:20pm 第 1 次编辑]
c=======================================================================72
subroutine gridgen2d_TTM(nbx,nby,nblocx,nblocy,
& IM,JM,linex,liney,X,Y,ISW,iprint)
c=======================================================================72
c c
c UUU UUU HH HH c
c UU UU HH HH c
c UU UU HH HH c
c UU UU HH HH c
c UU UU HHHHHHHHHHH c
c UU UU HH HH c
c UU UU HH HH c
c UU UU HH HH c
c UUUUUU HH HH c
c c
c USTC Hydrodynamics # "GridGen2D-TTM" c
c c
c 2D Numerical Body-Fitted Grid Generation using c
c Elliptic Partial Differential Equations (P.D.E.) c
c by Thompson, Thames, Mastin (J.C.P. Vol.15, 1974) c
c c
c-----------------------------------------------------------------------c
c Acknowledge: Teachers, researchers and students at USTC :-) c
c c
c The user may use and modify the code or give it to third parties, c
c provided that an acknowledgement to the USTC source of the original c
c version is retained. c
c=======================================================================72
c
c 2-D grid generation subroutine using elliptic P.D.E.
c to find source functions on boundaries
c
c Valid for geometry (transformed plane):
c +--+ +-------+ +-----------------+
c | | | | +---+ +---+ | |
c +--+ +--+ | +---+ | | | | | +--------+ |
c | | | | | +----+ | | | | |
c +--+ +--+ | +---+ | | | +---+ | |
c | | | | +------------+ | | | |
c +--+ +-------+ +-------+ +----+
c
c Variables:
c nbx,nby : number of x lines and y lines
c nblocx,nblocy : location arrays of grid line (excluding boundaries)
c nblocx(num,1) start location of x line "num"
c nblocx(num,2) end location of x line "num"
c nblocx(num,3) y location of the special x line
c linex,liney : location arrays of grid points refering nblocx,nblocy
c reverse transformation of nblocx,nblocy on the grid
c im,jm : number of rectangle grid points which contains
c the physical domain
c x,y : arrays of grid points
c isw : improvement parameter for orthogonal property
c iprint : debug parameter to print the iteration information
c
c References:
c Thmpson JF, Thames FC, Mastin CW, Automatic Numerical Generation of
c Body-Fitted Curvilinear Coordinate System for Fields Containing
c Any Number of Arbitrary Two-Dimensional Bodies, J. Comput. Phys.
c Vol.15, pp.299-329, 1974.
c Thomas PD, Middlecoff JF, Direct Control of the Grid Point Distribution
c in Meshes Generated by Elliptic Equations, AIAA J. Vol.18, p.652
c 1980.
c Thompson JF, Warsi ZUA, Mastin CW, Numerical Grid Generation,
c Foundation and Application, North-Holland, New York, 1985.
c
c=======================================================================72
c subroutine gridgen2d_ttm(nbx,nby,nblocx,nblocy,
c & IM,JM,linex,liney,X,Y,ISW,iprint)
c=======================================================================72
parameter (NI=250,NK=200)
c Params
dimension nblocx(nbx,3),nblocy(nby,3)
& ,linex(im,jm), liney(im,jm)
dimension X(IM,JM),Y(IM,JM)
c Locals
dimension PHI(NI,NK),PSI(NI,NK)
dimension A(NI),B(NI),C(NI)
dimension D1(NI),D2(NI),D3(NI),BETA(NI)
dimension I1(NI),I2(NI),J1(NI),J2(NI),K1(NI),K2(NI)
dimension RTX(NI),RTY(NI),RTZ(NI),RTM1(NI),RTM2(NI),RTM3(NI)
c
c-----------------------------------------------------------------------72
c
c.....Set the constants
c
PI=4.*ATAN(1.)
OMG=1.615
EPS=0.001
ErrAngle = 10.
NumIterPDE = 200
NumIterOrth = 100
c
c.....Calculate the coordinate transformation and
c find source functions PHI, PSI on boundaries
c
do ib = 1, nby
I = nblocy(ib,3)
jst = nblocy(ib,1)
jed = nblocy(ib,2)
do IDN = 1, 2
J = jst + (IDN-1)*(jed-jst)
XXI=(X(I+1,J)-X(I-1,J))/2.
YXI=(Y(I+1,J)-Y(I-1,J))/2.
XXX=(X(I-1,J)-2.*X(I,J)+X(I+1,J))
YXX=(Y(I-1,J)-2.*Y(I,J)+Y(I+1,J))
PHI(I,J)=-(XXI*XXX+YXI*YXX)/(XXI**2+YXI**2)
end do
do J = jst+1, jed-1
ET=FLOAT(J-jst)/FLOAT(jed-jst)
PHI(I,J)=(1.-ET)*PHI(I,jst)+ET*PHI(I,jed)
c...........! Initial guess
c y(i,j) = (1.-et)*y(i,jst)+et*y(i,jed)
end do
end do
do ib = 1, nbx
J = nblocx(ib,3)
ist = nblocx(ib,1)
ied = nblocx(ib,2)
do IDN = 1, 2
I = ist + (IDN-1)*(ied-ist)
XET=(X(I,J+1)-X(I,J-1))/2.
YET=(Y(I,J+1)-Y(I,J-1))/2.
XEE=(X(I,J-1)-2.*X(I,J)+X(I,J+1))
YEE=(Y(I,J-1)-2.*Y(I,J)+Y(I,J+1))
PSI(I,J)=-(XET*XEE+YET*YEE)/(XET**2+YET**2)
end do
do I = ist+1, ied-1
XI=FLOAT(I-ist)/FLOAT(ied-ist)
PSI(I,J)=(1.-XI)*PSI(ist,J)+XI*PSI(ied,J)
c...........! Initial guess
c x(i,j) = (1.-xi)*x(ist,j)+xi*x(ied,j)
end do
end do
c
c.....Iteration for Elliptic P.D.E. Iteration for Elliptic P.D.E.
c-----------------------------------------------------------------------72
c
ITEOUT=1
40 CONTINUE
N=1
50 CONTINUE
EPXM=0.
EPYM=0.
DSMIN=10000.
ARMIN=10000.
ARMAX=0.
c
c.....SLOR method for Elliptic P.D.E. along Y directional line,
c require the X lines should be continuous such as the above sketch
c but have no requirement for Y lines
c
do 100 ib = 1, nby
I = nblocy(ib,3)
jst = nblocy(ib,1)
jed = nblocy(ib,2)
do 60 J = jst+1, jed-1
XXI=(X(I+1,J)-X(I-1,J))/2.
YXI=(Y(I+1,J)-Y(I-1,J))/2.
XET=(X(I,J+1)-X(I,J-1))/2.
YET=(Y(I,J+1)-Y(I,J-1))/2.
XXE=(X(I+1,J+1)-X(I+1,J-1)-X(I-1,J+1)+X(I-1,J-1))/4.
YXE=(Y(I+1,J+1)-Y(I+1,J-1)-Y(I-1,J+1)+Y(I-1,J-1))/4.
ALP=XET**2+YET**2
BT=XXI*XET+YXI*YET
ist = nblocx(linex(i,j),1)
ied = nblocx(linex(i,j),2)
IF(I.EQ.ist+1 .OR. I.EQ.ied-1 .OR.
& J.EQ.jst+1 .OR. J.EQ.jed-1) BT=0.
GM=XXI**2+YXI**2
A(J-1)=GM*(1.-.5*PSI(I,J))
B(J-1)=-2.*(ALP+GM)
C(J-1)=GM*(1.+.5*PSI(I,J))
DB=-ALP*(X(I-1,J)+X(I+1,J))+2.*BT*XXE-ALP*PHI(I,J)*XXI
D1(J-1)=DB
IF(J.EQ.jst+1) D1(J-1)=DB-A(J-1)*X(I,jst)
IF(J.EQ.jed-1) D1(J-1)=DB-C(J-1)*X(I,jed)
DB=-ALP*(Y(I-1,J)+Y(I+1,J))+2.*BT*YXE-ALP*PHI(I,J)*YXI
D2(J-1)=DB
IF(J.EQ.jst+1) D2(J-1)=DB-A(J-1)*Y(I,jst)
IF(J.EQ.jed-1) D2(J-1)=DB-C(J-1)*Y(I,jed)
60 CONTINUE
A(jst)=0.
C(jed-2)=0.
c
c........Tri-diagonal equation solver, call the subroutine
c
CALL TRI2D(jed-jst-1,A(jst),B(jst),C(jst),D1(jst),D2(jst),
& BETA(jst),RTM1(jst),RTM2(jst),RTX(jst),RTY(jst))
c
c........Compute the new grid coordinates by relax iteration,
c and post error estimating
c
DO 90 J = jst+1, jed-1
XPRV=X(I,J)
YPRV=Y(I,J)
X(I,J)=OMG*RTX(J-1)+(1.-OMG)*X(I,J)
Y(I,J)=OMG*RTY(J-1)+(1.-OMG)*Y(I,J)
DS1=SQRT((X(I,J)-X(I-1,J))**2+(Y(I,J)-Y(I-1,J))**2)
DS2=SQRT((X(I,J)-X(I,J-1))**2+(Y(I,J)-Y(I,J-1))**2)
AR=DS1/DS2
DSLOC=AMIN1(DS1,DS2)
DSMIN=AMIN1(DSMIN,DSLOC)
ARMIN=AMIN1(ARMIN,AR)
ARMAX=AMAX1(ARMAX,AR)
IF(J.NE.jed-1) GO TO 70
DS22=SQRT((X(I,J)-X(I,J+1))**2+(Y(I,J)-Y(I,J+1))**2)
AR=DS1/DS22
DSLOC=AMIN1(DSLOC,DS22)
DSMIN=AMIN1(DSMIN,DSLOC)
ARMIN=AMIN1(ARMIN,AR)
ARMAX=AMAX1(ARMAX,AR)
70 IF(I.NE.nblocx(linex(i,j),2)-1) GO TO 80
DS12=SQRT((X(I,J)-X(I+1,J))**2+(Y(I,J)-Y(I+1,J))**2)
AR=DS12/DS2
DSLOC=AMIN1(DSLOC,DS12)
DSMIN=AMIN1(DSMIN,DSLOC)
ARMIN=AMIN1(ARMIN,AR)
ARMAX=AMAX1(ARMAX,AR)
IF(J.EQ.jed-1) ARMIN=AMIN1(ARMIN,DS12/DS22)
IF(J.EQ.jed-1) ARMAX=AMAX1(ARMAX,DS12/DS22)
80 EPX=ABS(X(I,J)-XPRV)/DSLOC
EPY=ABS(Y(I,J)-YPRV)/DSLOC
IF(EPX.GT.EPXM) EPXM=EPX
IF(EPY.GT.EPYM) EPYM=EPY
90 CONTINUE
100 CONTINUE
c
c.....Internal output options
c
c if (iprint.eq.1) WRITE (*,*) N,EPXM,EPYM,ADFM
IF(N.GT.NumIterPDE) GO TO 110
IF(EPXM.LT.EPS.AND.EPYM.LT.EPS) GO TO 110
N=N+1
GO TO 50
110 CONTINUE
c
c.....Y lines iteration end
c-----------------------------------------------------------------------72
c
IF(ISW.EQ.1) GO TO 150
c
c.....Orthogonal property requirement and improvement,
c need more outer iteration
c
ADFM=0.
ADFV=0.
do ib = 1, nby
I = nblocy(ib,3)
jst = nblocy(ib,1)
jed = nblocy(ib,2)
do IDN = 1, 2
J = jst + (IDN-1)*(jed-jst)
JV= jst + (IDN-1)*(jed-jst-1)
XXI=(X(I+1,J)-X(I-1,J))/2.
YXI=(Y(I+1,J)-Y(I-1,J))/2.
XET=(X(I,JV+1)-X(I,JV))
YET=(Y(I,JV+1)-Y(I,JV))
RX=SQRT(XXI**2+YXI**2)
RE=SQRT(XET**2+YET**2)
RXRE=XXI*XET+YXI*YET
COSA=RXRE/(RX*RE)
ANG=ACOS(COSA)
DIFA=0.5*PI-ANG
ADFM=AMAX1(ADFM,ABS(DIFA))
ADFV=ADFV+ABS(DIFA)
SGN=FLOAT(3-2*IDN)
PHI(I,J)=PHI(I,J)-SGN*0.1*TANH(DIFA)
end do
do J = jst+1, jed-1
ET=FLOAT(J-jst)/FLOAT(jed-jst)
PHI(I,J)=(1.-ET)*PHI(I,jst)+ET*PHI(I,jed)
end do
end do
do ib = 1, nbx
J = nblocx(ib,3)
ist = nblocx(ib,1)
ied = nblocx(ib,2)
do IDN = 1, 2
I = ist + (IDN-1)*(ied-ist)
IV= ist + (IDN-1)*(ied-ist-1)
XXI=(X(IV+1,J)-X(IV,J))
YXI=(Y(IV+1,J)-Y(IV,J))
XET=(X(I,J+1)-X(I,J-1))/2.
YET=(Y(I,J+1)-Y(I,J-1))/2.
RX=SQRT(XXI**2+YXI**2)
RE=SQRT(XET**2+YET**2)
RXRE=XXI*XET+YXI*YET
COSA=RXRE/(RX*RE)
ANG=ACOS(COSA)
DIFA=0.5*PI-ANG
ADFM=AMAX1(ADFM,ABS(DIFA))
SGN=FLOAT(3-2*IDN)
PSI(I,J)=PSI(I,J)-SGN*0.1*TANH(DIFA)
end do
do I = ist+1, ied-1
XI=FLOAT(I-ist)/FLOAT(ied-ist)
PSI(I,J)=(1.-XI)*PSI(ist,J)+XI*PSI(ied,J)
end do
end do
ADFM=ADFM*180./PI
ADFV=ADFV/(2.*IM)
if (iprint.eq.1)
& WRITE (*,*) 'IOUT=',ITEOUT,' DIF=',ADFM,' DIV=',ADFV
IF(ADFM.LT.ErrAngle.or.iteout.gt.NumIterOrth) GO TO 150
ITEOUT=ITEOUT+1
GO TO 40
150 continue
return
end
c=======================================================================72
SUBROUTINE TRI2D(N,A,B,C,F1,F2,BETA,Y1,Y2,X1,X2)
c=======================================================================72
DIMENSION A(N),B(N),C(N),X1(N),X2(N),F1(N),F2(N)
DIMENSION BETA(N),Y1(N),Y2(N)
c
BETA(1)=C(1)/B(1)
NM1=N-1
DO 10 I=2,NM1
BETA(I)=C(I)/(B(I)-A(I)*BETA(I-1))
10 CONTINUE
Y1(1)=F1(1)/B(1)
Y2(1)=F2(1)/B(1)
DO 20 I=2,N
Y1(I)=(F1(I)-A(I)*Y1(I-1))/(B(I)-A(I)*BETA(I-1))
Y2(I)=(F2(I)-A(I)*Y2(I-1))/(B(I)-A(I)*BETA(I-1))
20 CONTINUE
X1(N)=Y1(N)
X2(N)=Y2(N)
DO 30 K=1,NM1
X1(N-K)=Y1(N-K)-BETA(N-K)*X1(N-K+1)
X2(N-K)=Y2(N-K)-BETA(N-K)*X2(N-K+1)
30 CONTINUE
RETURN
END
上边是程序，请问程序中求解对三角用的是什么方法，也就是TRI2D(N,A,B,C,F1,F2,BETA,Y1,Y2,X1,X2)对应的子程序 |
|
IP卡
狗仔卡
发表于 2004-7-15 23:18:28
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
发表于 2004-7-15 23:23:11
楼主
