SUBROUTINE TFINTB(XB,YB,XL,YL,XT,YT,XR,YR,NL,NR,NBT, * NUMJ,X,Y) C***BEGIN PROLOGUE TFINTB C***DATE WRITTEN 000614 (YYMMDD) C***REVISION DATE 000614 (YYMMDD) C***CATEGORY NO. C***KEYWORDS TRANSFINITE BILINEAR INTERPOLATION, BLENDING FUNCTION C INTERPOLATION, ALGEBRAIC GRID GENERATION, BOUNDARY-FITTED C GRID GENERATION, NUMERICAL GRID GENERATION, CURVILINEAR C COORDINATE SYSTEM C C***AUTHOR SAUNDERS,B.V., (NIST) C***PURPOSE Generates a two dimensional boundary-fitted grid by C using transfinite bilinear interpolation. Grid points C on top and bottom boundaries are fixed to be the points C input by the user. C C***DESCRIPTION C C TFINTB generates a two dimensional boundary-fitted grid by C transfinite bilinear interpolation. In effect, TFINTB creates C a curvilinear coordinate system that maps the unit square onto C a region whose boundary is input by the user. The boundary must be C divided into four sides: bottom, left, top, right. Grid points C on top and bottom boundaries are fixed to be the points C input by the user so that the top and bottom boundaries are C interpolated accurately. The number of points on the top and bottom C boundaries must be the same. C (Also see TFINLR, TFINBLIN,TFINXACT) C C Advantage: C Top and bottom boundary grid points interpolate user inputted C top and bottom boundary points exactly. C C Disadvantage: C Number of grid points along top and bottom boundaries restricted. C C C INPUT C C XB,YB REAL(NBT), REAL(NBT) C Real vectors containing abcissas, ordinates C for bottom boundary points. C C XL,YL REAL(NL),REAL(NL) C Real vectors containing abcissas, ordinates C for left boundary points. C C XT,YT REAL(NBT),REAL(NBT) C Real vectors containing abcissas, ordinates C for top boundary points. C C XR,YR REAL(NR),REAL(NR) C Real vectors containing abcissas, ordinates C for right boundary points. C C NBT INTEGER C Number of bottom(top) boundary points. NBT C must be less than or equal to NMAXBP. C C NL INTEGER C Number of left boundary points. NL must be C less than or equal to NMAXBP. C C NR INTEGER C Number of right boundary points. NR must be C less than or equal to NMAXBP. C C NUMJ INTEGER C NBT x NUMJ equal the dimensions of the outputted C grid. C C C C OUTPUT C C X,Y REAL(NBT,NUMJ),REAL(NBT,NUMJ) C Matrices containing abcissas, ordinates of C grid points. First row of grid points would C be (X(1,1),Y(1,1)),....,(X(NBT,1),Y(NBT,1)). C Last row of grid points would be C (X(1,NUMJ),Y(1,NUMJ)),....,(X(NBT,NUMJ), C Y(NBT,NUMJ)) C C C PARAMETER C CONSTANT C C NMAXBP INTEGER C Maximum number of boundary points allowed C as input on any side of the domain. If size C insufficient for the needs of the user, user C can globally change parameter statement and C recompile. C C C C EXAMPLE PROGRAM C C PROGRAM TESTTB CC CC Sample driver for TFINTB. Computes NBT x NUMJ grid for nonconvex figure CC whose boundary data is in file input.data . CC C PARAMETER(NBT=3,NUMJ=5) C REAL X(NBT,NUMJ),Y(NBT,NUMJ),XB(100),YB(100),XL(100), C * YL(100),XT(100),YT(100),XR(100),YR(100) C OPEN(10, FILE='input.data') C OPEN(17, FILE='output.data') CC CC Number of points on top and bottom curves must be CC the same value, NBT. Consequently, only NL CC and NR need to be in the input file, but placing CC values for all 4 sides allows the same input CC format to be used for other transfinite routines CC TFINBLIN and TFINLR. Routine TFINXACT can also CC use the same format, but the first read ignores CC NT and NR. CC CC CC INPUT BOUNDARY POINTS CC C READ(10,*) NB,NL,NT,NR C READ(10,*) (XB(J),J=1,NB) C READ(10,*) (YB(J),J=1,NB) C READ(10,*) (XL(J),J=1,NL) C READ(10,*) (YL(J),J=1,NL) C READ(10,*) (XT(J),J=1,NT) C READ(10,*) (YT(J),J=1,NT) C READ(10,*) (XR(J),J=1,NR) C READ(10,*) (YR(J),J=1,NR) CC CHECK TO MAKE SURE NB = NT = NBT. IF NOT PRINT CC ERROR MESSAGE AND STOP. CC C IF(NB.NE.NBT.OR.NT.NE.NBT) THEN C PRINT*,'NUMBER OF POINTS ON TOP,BOTTOM BOUNDARIES MUST BE NBT' C STOP C ENDIF C CALL TFINTB(XB,YB,XL,YL,XT,YT,XR,YR,NL, C * NR,NBT,NUMJ,X,Y) C DO 500 J=1,NUMJ C DO 400 I=1,NBT C WRITE(17,*) X(I,J),Y(I,J) C 400 CONTINUE C WRITE(17,*) C 500 CONTINUE C STOP C END C C C input.data C C 3 4 3 4 C -2 0 4 C -3 0 2 C -2 -3 -3.5 -4 C -3 0 2 6 C -4 0 6 C 6 6 6 C 4 4.5 5 6 C 2 4 5 6 C C C output.data C C -2.000000 -3.000000 C 0.0000000E+00 0.0000000E+00 C 4.000000 2.000000 C C -2.731668 -0.8049936 C -0.2265644 1.529583 C 4.278540 3.114161 C C -3.355364 1.421455 C -0.3750565 3.065978 C 4.605251 4.210501 C C -3.713016 3.704127 C -0.2625465 4.571025 C 5.187923 5.187923 C C -4.000000 6.000000 C 0.0000000E+00 6.000000 C 6.000000 6.000000 C***REFERENCES Knupp,P., Steinberg, S., Fundamentals of Grid C Generation, CRC Press, 1993. C Thompson, J.F., Warsi, Z.U.A., Mastin, C.W., C Numerical Grid Generation: Foundations and C Applications, Elsevier, 1985. C Gordon, W.J. and Hall, C.A., "Constructions of C Curvilinear Coordinate Systems and Applications C to Mesh Generation," Int. Journal for Numer. C Methods in Eng., Vol. 7, 461-477, 1963. C C C***END PROLOGUE TFINTB C PARAMETER(NMAXBP=500) COMMON/LEFTRIGHT/S(2,NMAXBP),XLR(2,NMAXBP),YLR(2,NMAXBP) REAL XB(NBT),YB(NBT),XL(NL),YL(NL),XT(NBT),YT(NBT),XR(NR), * YR(NR),X(NBT,NUMJ),Y(NBT,NUMJ),A(NMAXBP),B(NMAXBP) C C FLEXIBLE LEFT, RIGHT BOUNDARIES C DO 5 J=1,NL XLR(1,J) = XL(J) YLR(1,J) = YL(J) 5 CONTINUE DO 10 J=1,NR XLR(2,J) = XR(J) YLR(2,J) = YR(J) 10 CONTINUE C C CONSTRUCT ARCLENGTH PARAMETRIZATIONS OF LEFT,RIGHT BOUNDARIES C DO 20 I=1,2 S(I,1)=0.0 20 CONTINUE DO 30 J=2,NL S(1,J)= SQRT((XL(J)-XL(J-1))**2 + * (YL(J)-YL(J-1))**2) +S(1,J-1) 30 CONTINUE DO 40 J=2,NR S(2,J)= SQRT((XR(J)-XR(J-1))**2 + * (YR(J)-YR(J-1))**2) +S(2,J-1) 40 CONTINUE DO 50 J=2,NL-1 S(1,J) = S(1,J)/S(1,NL) 50 CONTINUE DO 60 J=2,NR-1 S(2,J) = S(2,J)/S(2,NL) 60 CONTINUE S(1,NL)=1.0 S(2,NR)=1.0 C C COMPUTE VECTORS A AND B C A(1)=0.0 A(NBT)=1.0 DO 160 I=2,NBT-1 A(I)=FLOAT(I-1)/FLOAT(NBT-1) 160 CONTINUE B(1)=0.0 B(NUMJ)=1.0 DO 170 J=2,NUMJ-1 B(J)=FLOAT(J-1)/FLOAT(NUMJ-1) 170 CONTINUE C C C******************TRANSFINITE MAPPING ******************** C DO 500 J=1,NUMJ ETA = B(J) DO 400 I=1,NBT XI = A(I) X(I,J) = (1.0-XI)*XLEFT(ETA)+XI*XRIGHT(ETA) * + (1.0-ETA)*XB(I)+ETA*XT(I) * -(1.0-XI)*(1.0-ETA)*XB(1)-(1.0-XI)*ETA*XT(1) * -XI*(1.0-ETA)*XB(NBT)-XI*ETA*XT(NBT) Y(I,J) = (1.0-XI)*YLEFT(ETA)+XI*YRIGHT(ETA) * +(1.0-ETA)*YB(I)+ETA*YT(I) * -(1.0-XI)*(1.0-ETA)*YB(1)-(1.0-XI)*ETA*YT(1) * -XI*(1.0-ETA)*YB(NBT)-XI*ETA*YT(NBT) 400 CONTINUE 500 CONTINUE RETURN END REAL FUNCTION XLEFT(U) PARAMETER(NMAXBP=500) COMMON/LEFTRIGHT/S(2,NMAXBP),XLR(2,NMAXBP),YLR(2,NMAXBP) J=1 10 CONTINUE IF(U.GE.S(1,J).AND.U.LE.S(1,J+1)) THEN XLEFT=(XLR(1,J)*(S(1,J+1)-U)+XLR(1,J+1)*(U-S(1,J)))/ * (S(1,J+1)-S(1,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C REAL FUNCTION YLEFT(U) PARAMETER(NMAXBP=500) COMMON/LEFTRIGHT/S(2,NMAXBP),XLR(2,NMAXBP),YLR(2,NMAXBP) J=1 10 CONTINUE IF(U.GE.S(1,J).AND.U.LE.S(1,J+1)) THEN YLEFT=(YLR(1,J)*(S(1,J+1)-U)+YLR(1,J+1)*(U-S(1,J)))/ * (S(1,J+1)-S(1,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION XRIGHT(U) PARAMETER(NMAXBP=500) COMMON/LEFTRIGHT/S(2,NMAXBP),XLR(2,NMAXBP),YLR(2,NMAXBP) J=1 10 CONTINUE IF(U.GE.S(2,J).AND.U.LE.S(2,J+1)) THEN XRIGHT=(XLR(2,J)*(S(2,J+1)-U)+XLR(2,J+1)*(U-S(2,J)))/ * (S(2,J+1)-S(2,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION YRIGHT(U) PARAMETER(NMAXBP=500) COMMON/LEFTRIGHT/S(2,NMAXBP),XLR(2,NMAXBP),YLR(2,NMAXBP) J=1 10 CONTINUE IF(U.GE.S(2,J).AND.U.LE.S(2,J+1)) THEN YRIGHT=(YLR(2,J)*(S(2,J+1)-U)+YLR(2,J+1)*(U-S(2,J)))/ * (S(2,J+1)-S(2,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END