SUBROUTINE LAPGRID(XB,YB,XL,YL,XT,YT,XR,YR,NB,NL,NT,NR, * NUMI,NUMJ,X,Y,ALPHA,BETA,GAMMA,CENTER,TOL,OMEGA) C***BEGIN PROLOGUE LAPGRID C***DATE WRITTEN 980828 (YYMMDD) C***REVISION DATE 000811 (YYMMDD) C***CATEGORY NO. C***KEYWORDS LAPLACE GRID GENERATOR, ELLIPTIC GRID GENERATION, C BOUNDARY-FITTED GRID GENERATION, NUMERICAL GRID GENERATION C C***AUTHOR SAUNDERS,B.V., (NIST) C***PURPOSE Uses Laplace's equation to generate a two dimensional C boundary-fitted grid. C C C***DESCRIPTION C C LAPGRID uses Laplace's equation to generate a two dimensional C boundary-fitted grid. Essentially, this is a simple application C of the homogeneous generator popularized by Thompson, Thames and C Mastin. LAPGRID creates a curvilinear coordinate system that maps C the unit square (in XI,ETA plane) onto a region (in X,Y plane) C whose boundary is input by the user. The boundary must be divided C into four sides: bottom, left, top, right. The coordinate system C mapping (XI,ETA) -> (X(XI,ETA),Y(XI,ETA)) is determined by C requiring that the inverse transformation C (X,Y) -> (XI(X,Y),ETA(X,Y)) satisfy the Laplace equations: C C DEL^2(XI)=0, DEL^2(ETA)=0 C C with DIRICHLET boundary conditions defined by C C ETA=0, 1 on bottom, top boundaries, respectively C XI=0, 1 on left, right boundaries, respectively. C C To solve this system, we invert the Laplace equations to C obtain equations in X and Y over the unit square. The Dirichlet C boundary conditions are defined by the requirement that the C boundary of the unit square be mapped to the boundary of the C region in the X,Y plane. The inverted system is solved using C a successive overrelaxation (SOR) algorithm. 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. C C NL INTEGER C Number of left boundary points. C C NT INTEGER C Number of top boundary points. C C NR INTEGER C Number of right boundary points. C C NUMI,NUMJ INTEGER, INTEGER C NUMI x NUMJ equal grid dimensions desired by user. C C TOL REAL C Error tolerance for SOR routine. C C OMEGA REAL C Relaxation parameter. 0.0 < OMEGA < 2.0 C C ALPHA,BETA, C GAMMA,CENTER REAL(NUMI,NUMJ) C Arrays containing coefficents for finite C difference representation of equations. C The data for these is computed by the C routine. The user only needs to create C the space by passing four work arrays of C size NUMI x NUMJ as parameters to the routine. 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 EXAMPLE PROGRAM C C PROGRAM TESTLAPGRID CC CC Sample driver for LAPGRID. Computes NUMI x NUMJ grid for nonconvex figure CC whose boundary data is in file input.data . C 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),WORK1(NUMI,NUMJ), C * WORK2(NUMI,NUMJ),WORK3(NUMI,NUMJ),WORK4(NUMI,NUMJ), C * TOL,OMEGA C INTEGER NB,NL,NT,NR C OPEN(10, FILE='input.data') C OPEN(17, FILE='output.data') CC CC INPUT NUMBER OF POINTS ON EACH BOUNDARY CURVE CC C READ(10,*) NB,NL,NT,NR CC CC INPUT CURVES TO BE INTERPOLATED CC 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 TOL=1.0E-5 C OMEGA=1.9 C CALL LAPGRID(XB,YB,XL,YL,XT,YT,XR,YR,NB,NL,NT, C * NR,NUMI,NUMJ,X,Y,WORK1,WORK2,WORK3, C * WORK4,TOL,OMEGA) C 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 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.243465 0.5907210 C 0.4470828 1.697235 C 2.373153 2.455189 C 4.278540 3.114161 C C -3.355364 1.421455 C -1.445167 2.439929 C 0.5432621 3.205339 C 2.620951 3.743342 C 4.605251 4.210501 C C -3.713016 3.704127 C -1.512565 4.252871 C 0.7202415 4.657767 C 2.984865 4.938241 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 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 C***ROUTINES CALLED C C TFINBLIN. TFINBLIN can be replaced by any of the other C transfinite routines with NB,NT,NL,NR,NUMI,NUMJ adjusted C as indicated below.: C C TFINXACT -- In calling program, assign NB=NT=NUMI and C NL=NR=NUMJ. C TFINLR ---- In calling program, assign NL=NR=NUMJ. C TFINTB ---- In calling program, assign NB=NT=NUMI. C C Call statement for each of these routines is given below. C User must uncomment only the one that is to be used. C C***END PROLOGUE LAPGRID REAL XB(NB),YB(NB),XL(NL),YL(NL),XT(NT),YT(NT),XR(NR), * YR(NR),X(NUMI,NUMJ),Y(NUMI,NUMJ),ALPHA(NUMI,NUMJ), * BETA(NUMI,NUMJ),GAMMA(NUMI,NUMJ),CENTER(NUMI,NUMJ), * H,K,TOL,OMEGA,X_ETA,Y_ETA,X_XI,Y_XI INTEGER MAXCOUNT,NUMI,NUMJ,KOUNT,NB,NL,NT,NR C H=1.0/FLOAT(NUMI-1) K=1.0/FLOAT(NUMJ-1) C C COMPUTE INITIAL GRID USING TRANSFINITE BILINEAR INTERPOLATION C C CALL TFINBLIN(XB,YB,XL,YL,XT,YT,XR,YR,NB,NL,NT,NR, * NUMI,NUMJ,X,Y) C C CALL TFINXACT(XB,YB,XL,YL,XT,YT,XR,YR,NUMI,NUMJ,X,Y) C C CALL TFINLR(XB,YB,XL,YL,XT,YT,XR,YR,NB,NT, C * NUMI,NUMJ,X,Y) C C CALL TFINTB(XB,YB,XL,YL,XT,YT,XR,YR,NL,NR,NUMI,NUMJ,X,Y) C DO 20 J=2,NUMJ-1 DO 10 I=2,NUMI-1 X_ETA = 1.0/(2.0*K)*(X(I,J+1)-X(I,J-1)) Y_ETA = 1.0/(2.0*K)*(Y(I,J+1)-Y(I,J-1)) X_XI = 1.0/(2.0*H)*(X(I+1,J)-X(I-1,J)) Y_XI = 1.0/(2.0*H)*(Y(I+1,J)-Y(I-1,J)) ALPHA(I,J) = (X_ETA*X_ETA+Y_ETA*Y_ETA)/(H*H) BETA(I,J) = (-2.0*(X_XI*X_ETA+Y_XI*Y_ETA))/ * (4.0*H*K) GAMMA(I,J) = (X_XI*X_XI+Y_XI*Y_XI)/(K*K) CENTER(I,J)= -2.0*ALPHA(I,J)-2.0*GAMMA(I,J) 10 CONTINUE 20 CONTINUE c TOL=1.0E-5 c OMEGA=1.9 MAXCOUNT=1000 CALL SOR(ALPHA,BETA,GAMMA,CENTER,TOL,OMEGA,MAXCOUNT, * NUMI,NUMJ,X) CALL SOR(ALPHA,BETA,GAMMA,CENTER,TOL,OMEGA,MAXCOUNT, * NUMI,NUMJ,Y) RETURN END SUBROUTINE SOR(ALPHA,BETA,GAMMA,CENTER,TOL,OMEGA, * MAXCOUNT,NUMI,NUMJ,U) REAL U(NUMI,NUMJ),ALPHA(NUMI,NUMJ),BETA(NUMI,NUMJ), * GAMMA(NUMI,NUMJ),CENTER(NUMI,NUMJ),TOL,OMEGA INTEGER MAXCOUNT,NUMI,NUMJ,KOUNT C KOUNT=0 10 CONTINUE KOUNT=KOUNT+1 ZERROR=0.0 DO 60 J=2,NUMJ-1 DO 50 I=2,NUMI-1 SUM=ALPHA(I,J)*(U(I+1,J)+U(I-1,J))+ * GAMMA(I,J)*(U(I,J+1)+U(I,J-1))+ * BETA(I,J)*(U(I+1,J+1)+U(I-1,J-1) * -U(I-1,J+1)-U(I+1,J-1)) TEMP=-SUM/CENTER(I,J) TEMP=U(I,J)+OMEGA*(TEMP-U(I,J)) ZERROR=MAX(ZERROR,ABS(TEMP-U(I,J))) U(I,J)=TEMP 50 CONTINUE 60 CONTINUE IF(ZERROR.LE.TOL) THEN PRINT*, 'TOLERANCE MET AFTER',KOUNT,' ITERATIONS' print*, 'error is ',zerror ELSEIF(KOUNT.LT.MAXCOUNT) THEN GOTO 10 ELSE PRINT*, 'TOLERANCE NOT MET AFTER',KOUNT,' ITERATIONS' print*, 'error is ',zerror ENDIF RETURN END