SUBROUTINE TFINBLIN(XB,YB,XL,YL,XT,YT,XR,YR,NB,NL,NT,NR, * NUMI,NUMJ,X,Y) C***BEGIN PROLOGUE TFINBLIN C***DATE WRITTEN 980828 (YYMMDD) C***REVISION DATE 981007 (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. C C***DESCRIPTION C C TFINBLIN generates a two dimensional boundary-fitted grid by C transfinite bilinear interpolation. In effect, TFINBLIN 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. Unlike routine C TFINXACT which maps the square to the exact boundary points input C by the user, TFINBLIN maps the square to the piece-wise linear C interpolant of the boundary obtained by connecting the boundary C points. The grid is obtained by evaluating the transfinite mapping C over an equally spaced (in each coordinate direction) mesh on the C unit square. C (Also see TFINLR, TFINTB, TFINXACT) C C Advantages: C 1. It's quick. Easy to create input file and output grid. C Number of boundary points on each side can be different. C C 2. User has complete flexibility in choosing grid C dimensions C C 3. Can interpolate boundary points as close as needed C by increasing mesh size C C Disadvantages: C 1. Large mesh may be required to obtain accurate interpolation C of boundary C C INPUT C C XB,YB REAL(NB), REAL(NB) 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(NT),REAL(NT) 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 NB INTEGER C Number of bottom boundary points. NB must be C 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 NT INTEGER C Number of top boundary points. NT 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 NUMI,NUMJ INTEGER, INTEGER C NUMI x NUMJ equal grid dimensions desired by user. C C C C OUTPUT C C X,Y REAL(NUMI,NUMJ),REAL(NUMI,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(NUMI,1),Y(NUMI,1)). C Last row of grid points would be C (X(1,NUMJ),Y(1,NUMJ)),....,(X(NUMI,NUMJ), C Y(NUMI,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 EXAMPLE PROGRAM C C PROGRAM TESTBLIN CC CC Sample driver for TFINBLIN. Computes NUMI x NUMJ grid for nonconvex figure CC whose boundary data is in file input.data . CC C PARAMETER(NUMI=5,NUMJ=5) C REAL X(NUMI,NUMJ),Y(NUMI,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 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) C CALL TFINBLIN(XB,YB,XL,YL,XT,YT,XR,YR,NB,NL,NT, C * NR,NUMI,NUMJ,X,Y) C DO 500 J=1,NUMJ C DO 400 I=1,NUMI 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.8798263 -1.319740 C 0.3875484 0.1937742 C 2.193774 1.096887 C 4.000000 2.000000 C C -2.731668 -0.8049936 C -1.263986 0.4974899 C 0.3140972 1.674914 C 2.296319 2.394537 C 4.278540 3.114161 C C -3.355364 1.421455 C -1.555123 2.333847 C 0.3187177 3.162865 C 2.461984 3.686684 C 4.605251 4.210501 C C -3.713016 3.704127 C -1.582738 4.182641 C 0.5843406 4.619469 C 2.886131 4.903695 C 5.187923 5.187923 C C -4.000000 6.000000 C -1.500000 6.000000 C 1.000000 6.000000 C 3.500000 6.000000 C 6.000000 6.000000 C 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***ROUTINES CALLED C C***END PROLOGUE TFINBLIN C PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) REAL XB(NB),YB(NB),XL(NL),YL(NL),XT(NT),YT(NT),XR(NR), * YR(NR),X(NUMI,NUMJ),Y(NUMI,NUMJ),A(NMAXBP),B(NMAXBP) INTEGER N(4) C C********TRANSFINITE MAPPING -- X COMPONENT*********** C TX(XI,ETA) = (1.0-XI)*XLEFT(ETA)+XI*XRIGHT(ETA) * + (1.0-ETA)*XBOT(XI)+ETA*XTOP(XI) * -(1.0-XI)*(1.0-ETA)*XBOT(0.0)-(1.0-XI)*ETA*XTOP(0.0) * -XI*(1.0-ETA)*XBOT(1.0)-XI*ETA*XTOP(1.0) C C********TRANSFINITE MAPPING -- Y COMPONENT*********** C TY(XI,ETA) = (1.0-XI)*YLEFT(ETA)+XI*YRIGHT(ETA) * +(1.0-ETA)*YBOT(XI)+ETA*YTOP(XI) * -(1.0-XI)*(1.0-ETA)*YBOT(0.0)-(1.0-XI)*ETA*YTOP(0.0) * -XI*(1.0-ETA)*YBOT(1.0)-XI*ETA*YTOP(1.0) C C C CURVES TO BE INTERPOLATED C DO 10 J=1,NB CX(1,J) = XB(J) CY(1,J) = YB(J) 10 CONTINUE DO 11 J=1,NL CX(2,J) = XL(J) CY(2,J) = YL(J) 11 CONTINUE DO 12 J=1,NT CX(3,J) = XT(J) CY(3,J) = YT(J) 12 CONTINUE DO 13 J=1,NR CX(4,J) = XR(J) CY(4,J) = YR(J) 13 CONTINUE C C CONSTRUCT ARCLENGTH PARAMETRIZATIONS OF BOUNDARY SIDES C N(1)=NB N(2)=NL N(3)=NT N(4)=NR DO 20 I=1,4 T(I,1)=0.0 20 CONTINUE DO 40 I=1,4 DO 30 J=2,N(I) T(I,J)= SQRT((CX(I,J)-CX(I,J-1))**2 + * (CY(I,J)-CY(I,J-1))**2) +T(I,J-1) 30 CONTINUE 40 CONTINUE DO 60 I=1,4 DO 50 J=2,N(I)-1 T(I,J) = T(I,J)/T(I,N(I)) 50 CONTINUE 60 CONTINUE DO 70 I=1,4 T(I,N(I))=1.0 70 CONTINUE C C C COMPUTE GRIDPOINTS C C C COMPUTE VECTORS A AND B, RESPECTIVELY. C A(1)=0.0 A(NUMI)=1.0 DO 160 I=2,NUMI-1 A(I)=FLOAT(I-1)/FLOAT(NUMI-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 DO 210 I=1,NUMI XI = A(I) DO 200 J=1,NUMJ ETA = B(J) X(I,J) = TX(XI,ETA) Y(I,J) = TY(XI,ETA) 200 CONTINUE 210 CONTINUE RETURN END C C REAL FUNCTION XBOT(U) PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) J=1 10 CONTINUE IF(U.GE.T(1,J).AND.U.LE.T(1,J+1)) THEN XBOT=(CX(1,J)*(T(1,J+1)-U)+CX(1,J+1)*(U-T(1,J)))/ * (T(1,J+1)-T(1,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION YBOT(U) PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) J=1 10 CONTINUE IF(U.GE.T(1,J).AND.U.LE.T(1,J+1)) THEN YBOT=(CY(1,J)*(T(1,J+1)-U)+CY(1,J+1)*(U-T(1,J)))/ * (T(1,J+1)-T(1,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION XLEFT(U) PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) J=1 10 CONTINUE IF(U.GE.T(2,J).AND.U.LE.T(2,J+1)) THEN XLEFT=(CX(2,J)*(T(2,J+1)-U)+CX(2,J+1)*(U-T(2,J)))/ * (T(2,J+1)-T(2,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION YLEFT(U) PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) J=1 10 CONTINUE IF(U.GE.T(2,J).AND.U.LE.T(2,J+1)) THEN YLEFT=(CY(2,J)*(T(2,J+1)-U)+CY(2,J+1)*(U-T(2,J)))/ * (T(2,J+1)-T(2,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION XTOP(U) PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) J=1 10 CONTINUE IF(U.GE.T(3,J).AND.U.LE.T(3,J+1)) THEN XTOP=(CX(3,J)*(T(3,J+1)-U)+CX(3,J+1)*(U-T(3,J)))/ * (T(3,J+1)-T(3,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION YTOP(U) PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) J=1 10 CONTINUE IF(U.GE.T(3,J).AND.U.LE.T(3,J+1)) THEN YTOP=(CY(3,J)*(T(3,J+1)-U)+CY(3,J+1)*(U-T(3,J)))/ * (T(3,J+1)-T(3,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION XRIGHT(U) PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) J=1 10 CONTINUE IF(U.GE.T(4,J).AND.U.LE.T(4,J+1)) THEN XRIGHT=(CX(4,J)*(T(4,J+1)-U)+CX(4,J+1)*(U-T(4,J)))/ * (T(4,J+1)-T(4,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END C C REAL FUNCTION YRIGHT(U) PARAMETER(NMAXBP=500) COMMON/POINTS/T(4,NMAXBP),CX(4,NMAXBP),CY(4,NMAXBP) J=1 10 CONTINUE IF(U.GE.T(4,J).AND.U.LE.T(4,J+1)) THEN YRIGHT=(CY(4,J)*(T(4,J+1)-U)+CY(4,J+1)*(U-T(4,J)))/ * (T(4,J+1)-T(4,J)) RETURN ELSE J=J+1 ENDIF GOTO 10 END