C/////////////////////////////////////////////////////////////////////// C************* ALGORITHMS INTCOL & HERMCOL **************************** C//////////////////////////////////////////////////////////////////////// C//////////////// LOGICAL FILE 1 : DRIVER2 ///////////////////////////// C/////////////////////////////////////////////////////////////////////// C C PROGRAM DRIVE2 C C C ---- T E S T D R I V E R F O R R C T C O L ----- C C C MODIFIABLE PARAMETERS FOR DIMENSION SETTING. C THE FOLLOWING DIMENSIONS LIMIT THE SIZE OF THE PROBLEM THAT C CAN BE SOLVED. C C NGXMAX - MAXIMUM NUMBER OF X GRID LINES C NGYMAX - MAXIMUM NUMBER OF Y GRID LINES C NTBXMX - MAXIMUM NUMBER OF X GRID LINES IN TABLED OUTPUT ON A C GRID OTHER THAN THE DISCRETIZATION GRID C NTBYMX - MAXIMUM NUMBER OF Y GRID LINES IN TABLED OUTPUT ON A C GRID OTHER THAN THE DISCRETIZATION GRID C NOUTMX - MAXIMUM NUMBER OF OUTPUT REQUESTS C NEQMAX - MAXIMUM NUMBER OF EQUATIONS (= 4*NGXMAX*NGYMAX) ) C NCOEMX - MAXIMUM NUMBER OF COEFFICIENTS C NWKMAX - MAXIMUM AMOUNT OF WORKSPACE (= 2+(12*NGYMAX + 23)*NEQMAX ) ) C C C MAJOR ARRAYS C IDENTIFIER DIMENSION C GRIDX NGXMAX C GRIDY NGYMAX C COEF NEQMAX,NCOEMX C UNKVCT NEQMAX C WRKSP NWKMAX C TABX NTBXMX C TABY NTBYMX C IDCOEF NEQMAX,NCOEMX C NUMUNK NEQMAX C OUTFNC NOUTMX C OUTTYP NOUTMX C C NOTES C - THE DRIVER IS SET TO RUN HERMCOL ON EXAMPLE1 C - TO CHANGE TO INTCOL, SET BCTYPE TO UNCU IN THE DATA STMNT C - TO CHANGE THE GRID, CHANGE IGRID C - FOR EXAMPLE 2, CHANGE NOUT TO 4 AND USE FILE 7 INSTEAD OF FILE 6 C C C DIMENSIONS C REAL COEF(64,17),UNKVCT(64),WRKSP(4546),TABX(10),TABY(10) INTEGER IDCOEF(64,17),NUMUNK(64),OUTFNC(4),OUTTYP(4) LOGICAL HOMOBC CHARACTER *4 BCTYPE C C C********************* NON-STANDARD COMMON BLOCK USE **************** C MANY LABELED COMMON BLOCKS ARE DIMENSIONED WITH LENGTH 1 * C IN THIS ALGORITHM. THE CORRECT LENGTHS MUST BE SET IN THE * C DRIVER(MAIN) PROGRAM. THIS IS NOT STANDARD FORTRAN 77 BUT * C IT WORKS ON ALL SYSTEMS FORTRAN KNOWN TO US. IT ALLOWS ONE * C TO USE THE PROGRAMS REPEATEDLY WITHOUT RECOMPILING THEM. * C EXAMPLE. * C COMMON / GRIDXZ / GRIDX(1) * C IS USED INSTEAD OF * C COMMON / GRIDXZ / GRIDX(8) * C IF YOUR COMPILER ENFORCES THE STANDARD THEN YOU MUST CHANGE * C THESE COMMON BLOCK LENGTHS AND RECOMPILE EACH TIME * C********************* NON-STANDARD COMMON BLOCK USE **************** C C COMMON /PROBR/AXZ,BXZ,AYZ,BYZ COMMON /GRIDXZ/GRIDX(8) COMMON /GRIDYZ/GRIDY(8) COMMON /PROBI/NGRDXZ,NGRDYZ COMMON /INTEGS/NUMBEQ,NUMCOE,LEVELZ,MOUTPT,NROW C C DATA STATEMENTS TO SELECT PROBLEM TO RUN C C IGRID - INTEGER NOT 0 OR 1 C IF POSITIVE, THEN UNIFORM IGRID X IGRID C IF NEGATIVE, THEN VALUES LISTED IN SUBROUTINE SETGRD C AX,BX,AY,BY - LOWER AND UPPER LIMITS OF RECTANGULAR DOMAIN C DATA IGRID,AX,BX,AY,BY/4,0.,1.,0.,1./ C C TYPEBC - SELECT BOUNDARY CONDITION TYPE C FOR HERMCOL, MUST BE MIXD C FOR INTCOL, MUST BE NEUM OR UNCU FOR NEUMAN OR UNCOUPLED C HOMOBC - .TRUE. IF BOUNDARY CONDITIONS ARE HOMOGENEOUS C .FALSE. IF BOUNDARY CONTITIONS ARE NOT HOMOGEOUS C P1,P2 - PARAMETERS FOR BOUNDARY COLLOCATION POINTS(HERMCOL ONLY) C MUST HAVE (0 .LE. P1 .LT. P2 .LE 1.) C .AND. C .NOT. (P1 .EQ. 0. AND. P2 .EQ. 1.) C P1 = P2 = 0. SELECTS THE GAUSS POINTS C LEVEL - LEVEL OF OUTPUT DURING EXECUTION C DATA BCTYPE,HOMOBC,P1,P2,LEVEL/'MIXD',.FALSE.,0.,0.,1/ C MXNEQ = 64 C AXZ = AX BXZ = BX AYZ = AY BYZ = BY LEVELZ = LEVEL MOUTPT = 6 C CALL SETGRD(GRIDX,GRIDY,NGRIDX,NGRIDY,AX,BX,AY,BY,IGRID) NGRDXZ = NGRIDX NGRDYZ = NGRIDY C C SET OUTPUT TYPES C NOUT = 2 FOR EXAMPLE 1, = 4 FOR EXAMPLE 2 C NOUT = 4 C TABLE APPROXIMATE SOLUTION OUTFNC(1) = 1 OUTTYP(1) = 3 C TABLE RESIDUAL ON 10 X 10 GRID OUTFNC(2) = 3 OUTTYP(2) = 4 NTABX = 10 NTABY = 10 DO 10 I = 1,10 AIM1 = I - 1 TABX(I) = AX + (BX-AX)*AIM1/9. TABY(I) = AY + (BY-AY)*AIM1/9. 10 CONTINUE C MAXIMUM ERROR (MEANINGLESS FOR EXAMPLE 1 WHERE TRUE IS UNKNOWN) OUTFNC(3) = 2 OUTTYP(3) = 1 C MAXIMUM ERROR ON 10 X 10 GRID (MEANINGLESS FOR EXAMPLE 1) OUTFNC(4) = 2 OUTTYP(4) = 2 C CALL RCTCOL(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ,UNKVCT, . NUMUNK,WRKSP,HOMOBC,BCTYPE,P1,P2,OUTFNC,OUTTYP,NOUT, . TABX,TABY,NTABX,NTABY) STOP END SUBROUTINE SETUP(COEF,IDCOEF,MXNEQ,MXNCOE,NUMBEQ,NUMCOE,ABD,LDA, . UNKNWN,NBANDU,NBANDL) C C INTERFACE TO SOLUTION MODULE C REAL COEF(MXNEQ,*),UNKNWN(*),ABD(LDA,*) INTEGER IDCOEF(MXNEQ,*) C C FIND BANDWIDTHS AND SET MATRIX SIZE C IF (NBANDU*NBANDL.EQ.0) CALL BNDWTH(IDCOEF,MXNEQ,NUMBEQ,NUMCOE, . NBANDU,NBANDL) NROW = 2*NBANDL + NBANDU + 1 NCOL = NUMBEQ KWORK = NROW*NCOL + NCOL WRITE (6,9001) NUMBEQ,NBANDL,NBANDU,KWORK C C ZERO OUT ABD C DO 20 J = 1,NCOL DO 10 I = 1,NROW ABD(I,J) = 0.0 10 CONTINUE 20 CONTINUE C C LOAD ABD AND RIGHT SIDE C M = NBANDL + NBANDU + 1 DO 40 I = 1,NUMBEQ UNKNWN(I) = COEF(I,MXNCOE) DO 30 JJ = 1,NUMCOE J = IDCOEF(I,JJ) IF (J.EQ.0) GO TO 30 K = I - J + M ABD(K,J) = COEF(I,JJ) 30 CONTINUE 40 CONTINUE C C RETURN C 9001 FORMAT (/1H ,38H INTERFACE TO LINEAR EQUATION SOLVER/1H , . 22H NUMBER OF EQUATIONS,I9/1H ,18H LOWER BANDWIDTH,I13/1H , . 18H UPPER BANDWIDTH,I13/1H ,21H REQUIRED WORKSPACE,I10) END SUBROUTINE BNDWTH(IDCOEF,MXNEQ,NUMBEQ,NUMCOE,NBANDU,NBANDL) C C BNDWTH COMPUTES THE BANDWIDTH OF THE LINEAR SYSTEM WHOSE IDS C ARE STORED IN THE ARRAY IDCOEF C C NBANDU = NUMBER OF BANDS ABOVE DIAGONAL C NBANDL = NUMBER OF BANDS BELOW DIAGONAL C INTEGER IDCOEF(MXNEQ,*) C NBANDU = 0 NBANDL = 0 NCOEF = NUMCOE - 1 C DO 20 I = 1,NUMBEQ DO 10 JJ = 1,NCOEF J = IDCOEF(I,JJ) IF (J.EQ.0) GO TO 10 JIDIFF = J - I NBANDU = MAX0(NBANDU,JIDIFF) NBANDL = MIN0(NBANDL,JIDIFF) 10 CONTINUE 20 CONTINUE NBANDL = IABS(NBANDL) NBANDU = IABS(NBANDU) C RETURN END SUBROUTINE SOLVE(ABD,LDA,NUMBEQ,UNKNWN,PIVOTS,NBANDU,NBANDL) C C WAYNE R. DYKSEN, JANUARY 1982 C REAL ABD(LDA,*),UNKNWN(*),PIVOTS(*) C C FACTOR THE MATRIX C WRITE (6,9001) CALL FACTR(ABD,LDA,NUMBEQ,NBANDL,NBANDU,PIVOTS,INFO) IF (INFO.NE.0) GO TO 10 C C SOLVE THE SYSTEM C CALL BACKSU(ABD,LDA,NUMBEQ,NBANDL,NBANDU,PIVOTS,UNKNWN,0) C WRITE (6,9011) RETURN C C ERROR EXIT -- DIVISION BY ZERO IN BACKSU C 10 CONTINUE WRITE (MOUTPT,9021) INFO,INFO STOP C C 9001 FORMAT (/1H ,25H BAND GAUSS ELIMINATION) 9011 FORMAT (1H ,23H EXECUTION SUCCESSFUL) 9021 FORMAT (1H ,36H THE BAND FACTOR ROUTINE EXECUTED/1H , . 35H SUCCESSFULLY, BUT THE BAND BACK/1H , . 40H SOLVE ROUTINE WILL DIVIDE BY ZERO IF/1H , . 36H CALLED. THE DIAGONAL ELEMENT IN/1H ,8H ROW ,I10, . 20H OF THE UPPER FACTOR/1H ,12H IS ZERO.) END SUBROUTINE BACKSU(ABD,LDA,N,ML,MU,PIVOTS,B,JOB) C C P U R P O S E C C BACKSU SOLVES THE REAL BAND SYSTEM A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY FACTOR. C C D E S C R I P T I O N C C BACKSU IS A MODIFED VERSION OF THE LINPACK GENERAL BAND SOLVE C ROUTINE SGBSL WHICH USES AN INTEGER RATHER THAN A REAL PIVOT C VECTOR. C C R E F E R E N C E S C C J. J. DONGARRA, C. B. MOLER, J. R. BUNCH, AND G. W. STEWART, C LINPACK USERS' GUIDE. SIAM, PHILADELPHIA, 1979. C C A U T H O R C C WAYNE R. DYKSEN C DEPARTMENT OF COMPUTER SCIENCE C PURDUE UNIVERSITY C WEST LAFAYETTE, INDIANA 47097 C 317-494-6001 C C V E R S I O N C C JANUARY 1982 C C--------------------------------------------------------------------- C INTEGER LDA,N,ML,MU,JOB REAL ABD(LDA,*),B(*),PIVOTS(*) C C BACKSU SOLVES THE REAL BAND SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY FACTOR. C C ON ENTRY C C ABD REAL(LDA, N) C THE OUTPUT FROM FACTOR C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C C PIVOTS REAL(N) C THE PIVOT VECTOR FROM FACTOR. C C B REAL(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B , WHERE C TRANS(A) IS THE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF FACTOR HAS SET INFO .EQ. 0 . C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL FACTOR(ABD,LDA,N,ML,MU,PIVOTS,INFO) C IF (INFO .NE. 0) GO TO ... C DO 10 J = 1, P C CALL BACKSU(ABD,LDA,N,ML,MU,PIVOTS,C(1,J),0) C 10 CONTINUE C C SUBROUTINES AND FUNCTIONS C C BLAS SAXPY,SDOT C FORTRAN MIN0 C C INTERNAL VARIABLES C REAL SDOT,T INTEGER K,KB,L,LA,LB,LM,M,NM1 C M = MU + ML + 1 NM1 = N - 1 IF (JOB.NE.0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (ML.EQ.0) GO TO 30 IF (NM1.LT.1) GO TO 30 DO 20 K = 1,NM1 LM = MIN0(ML,N-K) L = PIVOTS(K) T = B(L) IF (L.EQ.K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL SAXPY(LM,T,ABD(M+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1,N K = N + 1 - KB B(K) = B(K)/ABD(M,K) LM = MIN0(K,M) - 1 LA = M - LM LB = K - LM T = -B(K) CALL SAXPY(LM,T,ABD(LA,K),1,B(LB),1) 40 CONTINUE GO TO 100 C 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1,N LM = MIN0(K,M) - 1 LA = M - LM LB = K - LM T = SDOT(LM,ABD(LA,K),1,B(LB),1) B(K) = (B(K)-T)/ABD(M,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (ML.EQ.0) GO TO 90 IF (NM1.LT.1) GO TO 90 DO 80 KB = 1,NM1 K = N - KB LM = MIN0(ML,N-K) B(K) = B(K) + SDOT(LM,ABD(M+1,K),1,B(K+1),1) L = PIVOTS(K) IF (L.EQ.K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE FACTR(ABD,LDA,N,ML,MU,PIVOTS,INFO) C C P U R P O S E C C FACTOR FACTORS A REAL BAND MATRIX BY GAUSS ELIMINATION. C C D E S C R I P T I O N C C FACTOR FACTORS A REAL BAND MATRIX BY GAUSS ELIMINATION WITH C PARTIAL PIVOTING AND ROW EQUILIBRATION. FACTOR IS A MODIFIED C VERSION OF THE LINPACK GENERAL BAND FACTOR ROUTINE SGBFA WHICH C DOES NOT DO ROW EQUILIBRATION. C C R E F E R E N C E S C C J. J. DONGARRA, C. B. MOLER, J. R. BUNCH, AND G. W. STEWART, C LINPACK USERS' GUIDE. SIAM, PHILADELPHIA, 1979. C C A U T H O R C C WAYNE R. DYKSEN C DEPARTMENT OF COMPUTER SCIENCE C PURDUE UNIVERSITY C WEST LAFAYETTE, INDIANA 47097 C 317-494-6001 C C V E R S I O N C C JANUARY 1982 C C----------------------------------------------------------------------- C C INTEGER LDA,N,ML,MU,INFO REAL ABD(LDA,*),PIVOTS(*) C C ON ENTRY C C ABD REAL(LDA, N) C CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS C OF THE MATRIX ARE STORED IN THE COLUMNS OF ABD AND C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS C ML+1 THROUGH 2*ML+MU+1 OF ABD . C SEE THE COMMENTS BELOW FOR DETAILS. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY ABD . C LDA MUST BE .GE. 2*ML + MU + 1 . C C N INTEGER C THE ORDER OF THE ORIGINAL MATRIX. C C ML INTEGER C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL. C 0 .LE. ML .LT. N . C C MU INTEGER C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL. C 0 .LE. MU .LT. N . C MORE EFFICIENT IF ML .LE. MU . C ON RETURN C C ABD AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C PIVOTS REAL(N) C A REAL VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT EGBSL WILL DIVIDE BY ZERO IF C CALLED. C C BAND STORAGE C C IF A IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT C WILL SET UP THE INPUT. C C ML = (BAND WIDTH BELOW THE DIAGONAL) C MU = (BAND WIDTH ABOVE THE DIAGONAL) C M = ML + MU + 1 C DO 20 J = 1, N C I1 = MAX0(1, J-MU) C I2 = MIN0(N, J+ML) C DO 10 I = I1, I2 C K = I - J + M C ABD(K,J) = A(I,J) C 10 CONTINUE C 20 CONTINUE C C THIS USES ROWS ML+1 THROUGH 2*ML+MU+1 OF ABD . C IN ADDITION, THE FIRST ML ROWS IN ABD ARE USED FOR C ELEMENTS GENERATED DURING THE TRIANGULARIZATION. C THE TOTAL NUMBER OF ROWS NEEDED IN ABD IS 2*ML+MU+1 . C THE ML+MU BY ML+MU UPPER LEFT TRIANGLE AND THE C ML BY ML LOWER RIGHT TRIANGLE ARE NOT REFERENCED. C C C SUBROUTINES AND FUNCTIONS C C BLAS SAXPY,SSCAL C MODIFIED BLAS ISWMAX C FORTRAN MAX0,MIN0 C C INTERNAL VARIABLES C REAL T INTEGER I,ISAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1 INTEGER KBEG,KEND,ISWMAX,IPVTK C C M = ML + MU + 1 INFO = 0 C C FIND RECIPROCAL OF LARGEST ELEMENT IN EACH ROW C FOR ROW EQUILIBRATION C DO 10 I = 1,N PIVOTS(I) = 0. 10 CONTINUE C DO 30 J = 1,N KBEG = M - MIN0(J-1,MU) KEND = M + MIN0(N-J,ML) DO 20 K = KBEG,KEND I = K + J - M PIVOTS(I) = AMAX1(PIVOTS(I),ABS(ABD(K,J))) 20 CONTINUE 30 CONTINUE C DO 40 I = 1,N IF (PIVOTS(I).NE.0.) PIVOTS(I) = 1./PIVOTS(I) 40 CONTINUE C C ZERO INITIAL FILL-IN COLUMNS C J0 = MU + 2 J1 = MIN0(N,M) - 1 IF (J1.LT.J0) GO TO 70 DO 60 JZ = J0,J1 I0 = M + 1 - JZ DO 50 I = I0,ML ABD(I,JZ) = 0.0E0 50 CONTINUE 60 CONTINUE 70 CONTINUE JZ = J1 JU = 0 C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C NM1 = N - 1 IF (NM1.LT.1) GO TO 170 DO 160 K = 1,NM1 KP1 = K + 1 C C ZERO NEXT FILL-IN COLUMN C JZ = JZ + 1 IF (JZ.GT.N) GO TO 90 IF (ML.LT.1) GO TO 90 DO 80 I = 1,ML ABD(I,JZ) = 0.0E0 80 CONTINUE 90 CONTINUE C C FIND L = PIVOT INDEX C LM = MIN0(ML,N-K) ISAMAX = ISWMAX(LM+1,ABD(M,K),PIVOTS(K),1) L = ISAMAX + M - 1 K1 = K + ISAMAX - 1 PIVOTS(K1) = PIVOTS(K) PIVOTS(K) = L + K - M C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (ABD(L,K).EQ.0.0E0) GO TO 140 C C INTERCHANGE IF NECESSARY C IF (L.EQ.M) GO TO 100 T = ABD(L,K) ABD(L,K) = ABD(M,K) ABD(M,K) = T 100 CONTINUE C C COMPUTE MULTIPLIERS C T = -1.0E0/ABD(M,K) CALL SSCAL(LM,T,ABD(M+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C IPVTK = PIVOTS(K) JU = MIN0(MAX0(JU,MU+IPVTK),N) MM = M IF (JU.LT.KP1) GO TO 130 DO 120 J = KP1,JU L = L - 1 MM = MM - 1 T = ABD(L,J) IF (L.EQ.MM) GO TO 110 ABD(L,J) = ABD(MM,J) ABD(MM,J) = T 110 CONTINUE CALL SAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1) 120 CONTINUE 130 CONTINUE GO TO 150 C 140 CONTINUE INFO = K 150 CONTINUE 160 CONTINUE 170 CONTINUE PIVOTS(N) = N IF (ABD(M,N).EQ.0.0E0) INFO = N RETURN END SUBROUTINE FNDINT(X1A,Y1A,X1B,Y1B,X2A,Y2A,X2B,Y2B,X0,Y0,NOINT) C C ............................................................... C C E L L P A C K O U T P U T M O D U L E C C SUBROUTINE FNDINT C C PURPOSE C TO FIND THE INTERSECTION POINT (X0,Y0) OF LINE L1 C DEFINED BY (X1A,Y1A) AND (X1B,Y1B) AND LINE L2 C DEFINED BY (X2A,Y2A) AND (X2B,Y2B). C C METHOD C THE SLOPE INTERCEPT FORM ( Y = M*X + B ) IS USED. C C ............................................................... C REAL X1A,Y1A,X1B,Y1B,X2A,Y2A,X2B,Y2B,X0,Y0,M1,B1,M2,B2,DX,DM, . EPSGRD COMMON /NUMCON/EPSGRD C NOINT = 0 NCASE = 1 C C IS L1 A VERTICAL LINE C DX = X1B - X1A IF (ABS(DX).GT.EPSGRD) GO TO 10 NCASE = 2 GO TO 20 C 10 CONTINUE M1 = (Y1B-Y1A)/DX B1 = Y1A - M1*X1A 20 CONTINUE C C IS L2 A VERTICAL LINE C DX = X2B - X2A IF (ABS(DX).GT.EPSGRD) GO TO 30 NCASE = NCASE + 2 GO TO 40 C 30 CONTINUE M2 = (Y2B-Y2A)/DX B2 = Y2A - M2*X2A C C BRANCH ON NCASE C C NCASE = 1 - BOTH LINES ARE NOT VERTICAL C NCASE = 2 - L1 IS VERTICAL, L2 IS NOT C NCASE = 3 - L2 IS VERTICAL, L1 IS NOT C NCASE = 4 - BOTH LINES ARE VERTICAL C 40 CONTINUE GO TO (50,60,70,80),NCASE C C BOTH LINES ARE NOT VERTICAL C 50 CONTINUE DM = M2 - M1 IF (ABS(DM).LE.EPSGRD) GO TO 80 X0 = (B1-B2)/DM Y0 = M1*X0 + B1 GO TO 90 C C L1 IS VERTICAL, L2 IS NOT C 60 CONTINUE X0 = X1A Y0 = M2*X0 + B2 GO TO 90 C C L2 IS VERTICAL, L1 IS NOT C 70 CONTINUE X0 = X2A Y0 = M1*X0 + B1 GO TO 90 C C L1 AND L2 ARE PARALLEL C 80 CONTINUE NOINT = 1 C 90 CONTINUE RETURN END SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY) C C OVERWRITE SINGLE PRECISION SY WITH SINGLE PRECISION SA*SX +SY. C REAL SX(*),SY(*),SA C IF (N.LE.0 .OR. SA.EQ.0.E0) RETURN IF (INCX.EQ.INCY) IF (INCX-1) 10,30,70 10 CONTINUE C C CODE FOR NONEQUAL OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 20 I = 1,N SY(IY) = SY(IY) + SA*SX(IX) IX = IX + INCX IY = IY + INCY 20 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 4. C 30 M = N - (N/4)*4 IF (M.EQ.0) GO TO 50 DO 40 I = 1,M SY(I) = SY(I) + SA*SX(I) 40 CONTINUE IF (N.LT.4) RETURN 50 MP1 = M + 1 DO 60 I = MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I+1) = SY(I+1) + SA*SX(I+1) SY(I+2) = SY(I+2) + SA*SX(I+2) SY(I+3) = SY(I+3) + SA*SX(I+3) 60 CONTINUE RETURN C C CODE FOR EQUAL, POSITIVE, NONUNIT INCREMENTS. C 70 CONTINUE NS = N*INCX DO 80 I = 1,NS,INCX SY(I) = SA*SX(I) + SY(I) 80 CONTINUE RETURN END REAL FUNCTION SDOT(N,SX,INCX,SY,INCY) C C RETURNS THE DOT PRODUCT OF SINGLE PRECISION SX AND SY. C REAL SX(1),SY(1) SDOT = 0.0E0 IF (N.LE.0) RETURN IF (INCX.EQ.INCY) IF (INCX-1) 10,30,70 10 CONTINUE C C CODE FOR UNEQUAL INCREMENTS OR NONPOSITIVE INCREMENTS. C IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO 20 I = 1,N SDOT = SDOT + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 20 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP SO REMAINING VECTOR LENGTH IS A MULTIPLE OF 5. C 30 M = N - (N/5)*5 IF (M.EQ.0) GO TO 50 DO 40 I = 1,M SDOT = SDOT + SX(I)*SY(I) 40 CONTINUE IF (N.LT.5) RETURN 50 MP1 = M + 1 DO 60 I = MP1,N,5 SDOT = SDOT + SX(I)*SY(I) + SX(I+1)*SY(I+1) + SX(I+2)*SY(I+2) + . SX(I+3)*SY(I+3) + SX(I+4)*SY(I+4) 60 CONTINUE RETURN C C CODE FOR POSITIVE EQUAL INCREMENTS .NE.1. C 70 CONTINUE NS = N*INCX DO 80 I = 1,NS,INCX SDOT = SDOT + SX(I)*SY(I) 80 CONTINUE RETURN END SUBROUTINE SSCAL(N,SA,SX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO 1. C JACK DONGARRA, LINPACK, 3/11/78. C REAL SA,SX(*) INTEGER I,INCX,M,MP1,N,NINCX C IF (N.LE.0) RETURN IF (INCX.EQ.1) GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX SX(I) = SA*SX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF (M.EQ.0) GO TO 40 DO 30 I = 1,M SX(I) = SA*SX(I) 30 CONTINUE IF (N.LT.5) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SX(I) = SA*SX(I) SX(I+1) = SA*SX(I+1) SX(I+2) = SA*SX(I+2) SX(I+3) = SA*SX(I+3) SX(I+4) = SA*SX(I+4) 50 CONTINUE RETURN END INTEGER FUNCTION ISWMAX(N,SX,WEIGHT,INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. WEIGHTED ABSOLUTE VALUE. C WAYNE DYKSEN, JANUARY 1982 C REAL SX(1),WEIGHT(1),SMAX INTEGER I,INCX,IX,N C ISWMAX = 0 IF (N.LT.1) RETURN ISWMAX = 1 IF (N.EQ.1) RETURN IF (INCX.EQ.1) GO TO 30 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 SMAX = ABS(SX(1)*WEIGHT(1)) IX = IX + INCX DO 20 I = 2,N IF (ABS(SX(IX)*WEIGHT(I)).LE.SMAX) GO TO 10 ISWMAX = I SMAX = ABS(SX(IX)*WEIGHT(I)) 10 IX = IX + INCX 20 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 30 SMAX = ABS(SX(1)*WEIGHT(1)) DO 40 I = 2,N IF (ABS(SX(I)*WEIGHT(I)).LE.SMAX) GO TO 40 ISWMAX = I SMAX = ABS(SX(I)*WEIGHT(I)) 40 CONTINUE RETURN END SUBROUTINE SETGRD(GRIDX,GRIDY,NGRIDX,NGRIDY,AX,BX,AY,BY,IGRID) C REAL GRIDX(*),GRIDY(*) C IF (IGRID.LT.0) GO TO 20 C C CODE TO SET A UNIFORM IGRID X IGRID GRID C DX = BX - AX DY = BY - AY NGRIDX = IGRID NGRIDY = IGRID C AIGRM1 = IGRID - 1 DO 10 I = 1,IGRID AIM1 = I - 1 GRIDX(I) = AX + (AIM1)*DX/ (AIGRM1) GRIDY(I) = AY + (AIM1)*DY/ (AIGRM1) 10 CONTINUE C RETURN C C C IF A NONUNIFORM GRID IS DESIRED, INSERT CODE TO SET C NGRIDX, NGRIDY, GRIDX(1..NGRIDX), GRIDY(1..NGRIDY) C 20 CONTINUE RETURN END SUBROUTINE TABLER(OUTFNC,TABX,NTABX,TABY,NTABY,TPRINT,UNKNWN) C C C PURPOSE C C PRINTS THE VALUES OF THE OUTPUT FUNCTION ON THE GRID C DEFINED BY THE ARRAYS TABX AND TABY IF A PAR- C TICULAR GRID POINT IS INSIDE THE REGION. C C PARAMETERS C C OUTFNC - INDICATES WHAT FUNCTION IS TO BE OUTPUT C =1 APPROXIMATE SOLUTION C =2 ERROR (TRUE MUST BE SUPPLIED) C =3 RESIDUAL C C TABX, TABY - X AND Y COORDINATES FOR OUTTYP 2 AND 4 C C NTABX, NTABY - NUMBER OF VALUES IN TABX, TABY C C TPRINT - WORKSPACE OF DIMENSION >= NTABX C C ALL OTHER PARAMETERS ARE PASSED TO SUBROUTINES C C REAL TABX(NTABX),TABY(NTABY),DERVSL(6),UNKNWN(*),TPRINT(NTABX) C INTEGER OUTFNC C COMMON /GRIDXZ/GRIDX(1) COMMON /GRIDYZ/GRIDY(1) COMMON /INTEGS/NUMBEQ,NUMCOE,LEVEL,MOUTPT,NROW C C C PRINT HEADING C IF (OUTFNC.EQ.1) WRITE (MOUTPT,9001) NTABX,NTABY IF (OUTFNC.EQ.2) WRITE (MOUTPT,9011) NTABX,NTABY IF (OUTFNC.EQ.3) WRITE (MOUTPT,9021) NTABX,NTABY WRITE (MOUTPT,9031) WRITE (MOUTPT,9051) (TABX(I),I=1,NTABX) C C DO 20 JUP = 1,NTABY J = NTABY - JUP + 1 Y = TABY(J) WRITE (MOUTPT,9041) Y C DO 10 I = 1,NTABX TPRINT(I) = 0.0 X = TABX(I) IPNTYP = 0 IF (IPNTYP.EQ.3) GO TO 10 IF (OUTFNC.NE.3) TPRINT(I) = COLAPR(X,Y,UNKNWN,DERVSL) IF (OUTFNC.EQ.2) TPRINT(I) = TPRINT(I) - TRUE(X,Y) IF (OUTFNC.EQ.3) TPRINT(I) = RESID(X,Y,IPNTYP,UNKNWN) 10 CONTINUE C WRITE (MOUTPT,9051) (TPRINT(I),I=1,NTABX) 20 CONTINUE C RETURN C 9001 FORMAT (///1H ,10X,46 (1H+)/1H ,10X,1H+,44X,1H+/1H ,10X,1H+,4X, . 9HTABLE OF ,8HSOLUTION,3H ON,I4,2H X,I4, . 6H GRID,4X,1H+) 9011 FORMAT (///1H ,10X,46 (1H+)/1H ,10X,1H+,44X,1H+/1H ,10X,1H+,4X, . 9HTABLE OF ,8HERROR ,3H ON,I4,2H X,I4, . 6H GRID,4X,1H+) 9021 FORMAT (///1H ,10X,46 (1H+)/1H ,10X,1H+,44X,1H+/1H ,10X,1H+,4X, . 9HTABLE OF ,9HRESIDUALS,3H ON,I4,2H X,I4, . 6H GRID,4X,1H+) 9031 FORMAT (1H ,10X,1H+,44X,1H+/1H ,10X,46 (1H+)///1H ,4X, . 15HX-ABSCISSAE ARE/1H ,4X,15 (1H-)) 9041 FORMAT (/1H ,7X,3HY =,E13.6/1H ,7X,16 (1H-)) 9051 FORMAT ((1H ,3X,E13.6,3X,E13.6,3X,E13.6,3X,E13.6)) END SUBROUTINE FNCMAX(OUTFNC,TABX,NTABX,TABY,NTABY,UNKNWN) C C C PURPOSE C C COMPUTE AND PRINT MAX, L1 AND L2 NORMS OF THE OUTPUT C FUNCTION BASED ON THE GRID DEFINED BY TABX AND TABY C C PARAMETERS C C OUTFNC - INDICATES WHAT FUNCTION IS TO BE OUTPUT C =1 APPROXIMATE SOLUTION C =2 ERROR (TRUE MUST BE SUPPLIED) C =3 RESIDUAL C C TABX, TABY - GRID ON WHICH TO MEASURE NORMS C C NTABX, NTABY - NUMBER OF VALUES IN TABX, TABY C REAL TABX(NTABX),TABY(NTABY),DERVSL(6),UNKNWN(*) C INTEGER OUTFNC C COMMON /GRIDXZ/GRIDX(1) COMMON /GRIDYZ/GRIDY(1) COMMON /INTEGS/NUMBEQ,NUMCOE,LEVEL,MOUTPT,NROW C C R1NRMI = 0.0 R1NRM1 = 0.0 R1NRM2 = 0.0 GPTS = 0 C C DO 20 I = 1,NTABX X = TABX(I) DO 10 J = 1,NTABY Y = TABY(J) IPTYPE = 0 IF (IPTYPE.EQ.3) GO TO 10 IF (OUTFNC.NE.3) FVALUE = COLAPR(X,Y,UNKNWN,DERVSL) IF (OUTFNC.EQ.2) FVALUE = ABS(FVALUE-TRUE(X,Y)) IF (OUTFNC.EQ.3) FVALUE = RESID(X,Y,IPTYPE,UNKNWN) GPTS = GPTS + 1. R1NRMI = AMAX1(R1NRMI,FVALUE) R1NRM1 = R1NRM1 + FVALUE R1NRM2 = R1NRM2 + FVALUE**2 10 CONTINUE 20 CONTINUE R1NRM1 = R1NRM1/GPTS R1NRM2 = SQRT(R1NRM2/GPTS) IF (OUTFNC.EQ.1) WRITE (MOUTPT,9001) NTABX,NTABY,R1NRMI,NTABX, . NTABY,R1NRM1,NTABX,NTABY,R1NRM2 IF (OUTFNC.EQ.2) WRITE (MOUTPT,9011) NTABX,NTABY,R1NRMI,NTABX, . NTABY,R1NRM1,NTABX,NTABY,R1NRM2 IF (OUTFNC.EQ.3) WRITE (MOUTPT,9021) NTABX,NTABY,R1NRMI,NTABX, . NTABY,R1NRM1,NTABX,NTABY,R1NRM2 RETURN C C 9001 FORMAT (//5X,60 (1H+)/5X,1H+,58X,1H+/5X,12H+ MAX( ABS(, . 8HSOLUTION,7H) ) ON ,I3,3H X ,I3,7H GRID =, . 1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 12H+ L1 NORM( ,8HSOLUTION,7H) ON ,I3,3H X , . I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 12H+ L2 NORM( ,8HSOLUTION,7H) ON ,I3,3H X , . I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 60 (1H+)) 9011 FORMAT (//5X,60 (1H+)/5X,1H+,58X,1H+/5X,12H+ MAX( ABS(, . 8HERROR ,7H) ) ON ,I3,3H X ,I3,7H GRID =, . 1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 12H+ L1 NORM( ,8HERROR ,7H) ON ,I3,3H X , . I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 12H+ L2 NORM( ,8HERROR ,7H) ON ,I3,3H X , . I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 60 (1H+)) 9021 FORMAT (//5X,60 (1H+)/5X,1H+,58X,1H+/5X,12H+ MAX( ABS(, . 8HRESIDUAL,7H) ) ON ,I3,3H X ,I3,7H GRID =, . 1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 12H+ L1 NORM( ,8HRESIDUAL,7H) ON ,I3,3H X , . I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 12H+ L2 NORM( ,8HRESIDUAL,7H) ON ,I3,3H X , . I3,7H GRID =,1PE14.7,2X,1H+/5X,1H+,58X,1H+/5X, . 60 (1H+)) END SUBROUTINE ONLINE(X,Y,GT,GY,XBOUND,YBOUND,NBNDPT,IPNTYP,IERR) C C ............................................................... C C E L L P A C K O U T P U T M O D U L E C C SUBROUTINE ONLINE C C PURPOSE C TO DETERMINE IF THE POINT (X,Y) IS INSIDE THE REGION, C ON THE BOUNDARY, OR OUTSIDE THE REGION WHEN (X,Y) IS C ON A VERTICAL GRID LINE. C C PARAMETERS C X - THE X COORDINATE OF THE POINT C Y - THE Y COORDINATE OF THE POINT C GT - GT(1) IS THE GTYPE OF THE LOWER GRID POINT C (X,GY(1)) WHILE GT(2) IS THE GTYPE OF THE C THE UPPER GRID POINT (X,GY(2)) C GY - GY(1) IS THE Y COORDINATE OF THE LOWER GRID C POINT WHILE GY(2) IS THE Y COORDINATE OF THE C UPPER GRID POINT. C IPNTYP - 1, IF (X,Y) IS INSIDE THE REGION C 2, IF (X,Y) IS ON THE BOUNDARY C 3, IF (X,Y) IS OUTSIDE THE REGION C IERR - ERROR RETURN CODE. C C ............................................................... C REAL X,Y,GY(2),XBOUND(NBNDPT),YBOUND(NBNDPT),PLEN,BLEN INTEGER GT(2),IGT(2) COMMON /NUMCON/EPSGRD C C NINSID - THE NUMBER OF NEIGHBORING GRID POINTS C INSIDE THE DOMAIN C NOUTSD - THE NUMBER OF NEIGHBORING GRID POINTS C OUTSIDE THE DOMAIN C NINSID = 0 NOUTSD = 0 C DO 20 K = 1,2 IF (GT(K).LT.999) GO TO 10 NINSID = NINSID + 1 GO TO 20 10 CONTINUE IF (GT(K).LE.0) NOUTSD = NOUTSD + 1 20 CONTINUE C NCASE = 3*NOUTSD + NINSID + 1 C C BRANCH ON NCASE C ----------------------------------- C NOUTSD NINSID NCASE (X,Y) C ----------------------------------- C C 0 0 1 BOUNDARY C 0 1 2 INSIDE C 0 2 3 INSIDE C 1 0 4 OUTSIDE C 1 1 5 UNKNOWN C 1 2 6 IMPOSS. C 2 0 7 OUTSIDE C GO TO (100,90,90,110,30,120,110),NCASE C C THIS IS THE CASE WHERE ONE OF THE NEIGHBORING GRID C POINTS IS INSIDE WHILE THE OTHER IS OUTSIDE. C C FIRST, TRY TO USE GTYPE INFO SAVED IN GT TO FIND THE C BOUNDARY - GRID LINE INTERSECTION POINT. C 30 CONTINUE C DO 40 K = 1,2 KB = MOD(IABS(GT(K)),1000) IF (ABS(XBOUND(KB)-X).GT.EPSGRD) GO TO 40 IF (YBOUND(KB).LT.GY(1)) GO TO 40 IF (YBOUND(KB).GT.GY(2)) GO TO 40 GO TO 60 C 40 CONTINUE C C USING GTYPE HAS FAILED. C NOW TRY A BRUTE FORCE SEARCH. C DO 50 KB = 1,NBNDPT IF (ABS(XBOUND(KB)-X).GT.EPSGRD) GO TO 50 IF (YBOUND(KB).LT.GY(1)) GO TO 50 IF (YBOUND(KB).GT.GY(2)) GO TO 50 GO TO 60 50 CONTINUE C IERR = 4 GO TO 120 C C PLEN - THE DISTANCE FROM THE INTERIOR POINT TO (X,Y) C C BLEN - THE DISTANCE FROM THE INTERIOR POINT TO THE C BOUNDARY - GRID LINE INTERSECTION POINT C 60 CONTINUE IF (GT(2).GE.999) GO TO 70 PLEN = ABS(GY(1)-Y) BLEN = ABS(GY(1)-YBOUND(KB)) GO TO 80 C 70 CONTINUE PLEN = ABS(GY(2)-Y) BLEN = ABS(GY(2)-YBOUND(KB)) C C NOW DETERMINE THE RELATIONSHIP C OF (X,Y) TO THE BOUNDARY. C 80 CONTINUE IF (ABS(PLEN-BLEN).LE.EPSGRD) GO TO 100 IF (PLEN.LT.BLEN) GO TO 90 GO TO 110 C C (X,Y) IS INSIDE C 90 CONTINUE IPNTYP = 1 GO TO 120 C C (X,Y) IS ON THE BOUNDARY C 100 CONTINUE IPNTYP = 2 GO TO 120 C C (X,Y) IS OUTSIDE C 110 CONTINUE IPNTYP = 3 C 120 CONTINUE RETURN END REAL FUNCTION RESID(X,Y,IPTCOM,UNKNWN) C C C PURPOSE C TO DETERMINE THE VALUE OF THE RESIDUAL AT C THE POINT (X,Y). C C PARAMETERS C C X,Y - THE POINT AT WHICH TO COMPUTE THE RESIDUAL C IPTCOM - RESULT FROM SUBROUTINE OUTSID C =1 POINT IS INTERIOR TO DOMAIN C =2 POINT IS ON BOUNDARY OF DOMAIN C C REAL UNKNWN(*),BVALUS(4),DERVSL(6),CVALUS(7) C C COMMON /BRANZZ/BRANGE(2,1) COMMON /GRIDXZ/GRIDX(1) COMMON /GRIDYZ/GRIDY(1) COMMON /NUMCON/EPSGRD COMMON /PROBR/AX,BX,AY,BY C RESID = 0. C C CHECK TO SEE IF (X,Y) IS ON THE BOUNDARY. C IN THIS CASE THE RESIDUAL OF THE BOUNDARY C OPERATOR MUST BE COMPUTED. C IPTCOM = 1 IF (X.EQ.AX .OR. X.EQ.BX .OR. Y.EQ.AY .OR. Y.EQ.BY) IPTCOM = 2 C 10 IF (IPTCOM.EQ.2) GO TO 20 C C EVALUATE THE COEFFICIENTS OF THE PDE OPERATOR. C CALL PDE(X,Y,CVALUS) GO TO 30 C C RECTANGULAR CASE C 20 IF (X.EQ.BX) IP = 1 IF (Y.EQ.AY) IP = 2 IF (X.EQ.AX) IP = 3 IF (Y.EQ.BY) IP = 4 TEMP = BCOND(IP,X,Y,BVALUS) C C EVALUATE COMPUTED SOLUTION AND ITS DERIVATIVES AT (X,Y) C 30 CONTINUE TEMP = COLAPR(X,Y,UNKNWN,DERVSL) C C DERIVATIVES OF SOLUTION ARE NOW IN THE ARRAY DERVSL C C COMPUTE THE RESIDUAL FOR ALL CASES C IF (IPTCOM.EQ.2) GO TO 40 C C RESIDUAL OF PDE OPERATOR IS C RESID = CVALUS(1)*DERVSL(1) + CVALUS(2)*DERVSL(2) + . CVALUS(3)*DERVSL(3) + CVALUS(4)*DERVSL(4) + . CVALUS(5)*DERVSL(5) + CVALUS(6)*DERVSL(6) - PDERHS(X,Y) GO TO 50 C C THE RESIDUAL OF THE BOUNDARY OPERATOR IS C 40 CONTINUE RESID = BVALUS(1)*DERVSL(6) + BVALUS(2)*DERVSL(4) + . BVALUS(3)*DERVSL(5) - BVALUS(4) C 50 CONTINUE RETURN END SUBROUTINE OUTPUT(OUTFNC,OUTTYP,TABX,NTABX,TABY,NTABY,WKSPAC, . UNKNWN) C C PURPOSE C C CONTROLS THE TYPE OF OUTPUT OF THE COLLOCATION APPROX. C C PARAMETERS C C OUTFNC - INDICATES WHAT FUNCTION IS TO BE OUTPUT C =1 APPROXIMATE SOLUTION C =2 ERROR (TRUE MUST BE SUPPLIED) C =3 RESIDUAL C C OUTTYP - INDICATES WHAT INFORMATION ABOUT THE C OUTPUT FUNCTION IS TO BE PRINTED C =1 PRINT MAX, L1, L2 NORMS OF FUNCTION C BASED ON DISCRETIZATION GRID C =2 PRINT MAX, L1, L2 NORMS OF FUNCTION C BASED ON GIVEN GRID (TABX,TABY) C =3 PRINT TABLE OF FUNCTION ON DISC. GRID C =4 PRINT TABLE OF FUNCTION ON GIVEN GRID C C TABX, TABY - X AND Y COORDINATES FOR OUTTYP 2 AND 4 C C NTABX, NTABY - NUMBER OF VALUES IN TABX, TABY C C ALL OTHER PARAMETERS ARE PASSED TO SUBROUTINES C C REAL TABX(NTABX),TABY(NTABY),WKSPAC(NTABX),UNKNWN(*) INTEGER OUTFNC,OUTTYP C COMMON /GRIDXZ/GRIDX(1) COMMON /GRIDYZ/GRIDY(1) COMMON /PROBI/NGRIDX,NGRIDY C GO TO (10,20,30,40),OUTTYP C C PRINT NORMS OF THE FUNCTION ON DISCRETIZATION GRID C 10 CALL FNCMAX(OUTFNC,GRIDX,NGRIDX,GRIDY,NGRIDY,UNKNWN) RETURN C C PRINT NORMS OF FUNCTION ON GIVEN GRID C 20 CALL FNCMAX(OUTFNC,TABX,NTABX,TABY,NTABY,UNKNWN) RETURN C C PRINT TABLE OF FUNCTION ON DISCRETIZATION GRID C 30 CALL TABLER(OUTFNC,GRIDX,NGRIDX,GRIDY,NGRIDY,WKSPAC,UNKNWN) RETURN C C PRINT TABLE OF FUNCTION ON GIVEN GRID C 40 CALL TABLER(OUTFNC,TABX,NTABX,TABY,NTABY,WKSPAC,UNKNWN) RETURN END C////////////////////////////////////////////////////////////////// C*********** ALGORITHMS INTCOL & HERMCOL ************************* C////////////////////////////////////////////////////////////////// C>>>>>>>>>>> LOGICAL FILE 2: ALGORITHM INTCOL <<<<<<<<<<<<<<<<<<<< C////////////////////////////////////////////////////////////////// C SUBROUTINE BCINTR(ISIDE,VISIDE,Q,NQD,R,NRD,T,KNOTS) C C **** DETERMINES THE 1-D HERMITE INTERPOLANT OF DIRICHLET OR C NEUMANN TYPE BOUNDARY CONDITIONS ***** C C ---- INPUT VARIABLES:: ISIDE : BOUNDARY SIDE C RÕ1:NRDå : ARRAY OF GAUSS POINTS ,BASIS C FUNCTIONS AND THEIR DERIVATIVES C DEFINED IN S/R INTPTS C TÕ1:KNOTSå : NODES OF A SIDE BOUNDARY C C ---- OUTPUT VARIABLES:: QÕ1:NQDå : ARRAY WITH THE DEGREES OF FREEDOM C OF THE HERMITE INTERPOLANT OF B.C C C **NOTE** WE ASSUME THAT 'U' OR 'DU/DN' ARE DEFINED IN C SUBINTERVALS OF A BOUNDARY SIDE WHOSE END POINTS C ARE NODAL POINTS C C C ---- DECLARATIONS ------ C REAL Q(NQD),R(NRD),T(KNOTS),BVALUS(4),BVLUS1(4),RS(4) INTEGER TYPEBC COMMON /INTEGS/NUMBEQ,NUMCOE,LEVEL,MOUTPT,NROW C C DETERMINES HERMITE INTERPOLANT IN THE INTERVAL (T(1),T(1+1/2)) C USING BOUNDARY VALUES AT T(1), TAU1, TAU2, TM = (T(1)+T(I+1))/2 C WHERE TAU1,TAU2 ARE THE GAUSS POINTS OF THE INTERVAL C C CASE I : NODE = 1 C RS(1) = BRHS(ISIDE,VISIDE,T(1),BVALUS) C TAU1L = (T(1)+R(1))*.5 TAU2L = (T(1)+R(2))*.5 RS(2) = BRHS(ISIDE,VISIDE,TAU1L,BVALUS) RS(3) = BRHS(ISIDE,VISIDE,TAU2L,BVALUS) C TMDLL = (T(1)+T(2))*.5 VTMDLL = BRHS(ISIDE,VISIDE,TMDLL,BVLUS1) RS(4) = VTMDLL C DETINV = 1./ (R(5)*R(18)-R(17)*R(6)) C C DETERMINES THE TWO DEGREES OF FREEDOM ASSOCIATED WITH C NODE = 1 C Q(1) = RS(1) B1 = RS(2) - RS(1)*R(3) - RS(4)*R(4) B2 = RS(3) - RS(1)*R(15) - RS(4)*R(16) Q(2) = 2.* (B1*R(18)-B2*R(6))*DETINV C C DETERMINES THE TWO DEGREES OF FREEDOM ASSOCIATED WITH C NODE I C C CASE NODE = 2,KNOTS -1 C KM1 = KNOTS - 1 DO 30 I = 2,KM1 C VBCND = BRHS(ISIDE,VISIDE,T(I),BVALUS) TMDLR = (T(I)+T(I+1))*.5 VTMDLR = BRHS(ISIDE,VISIDE,TMDLR,BVLUS1) C C CHECK THE TYPE OF BOUNDARY CONDITION AT THE NODE I C TYPEBC = 0 IF (BVALUS(1).NE.0.) TYPEBC = 1 IF (BVALUS(2).NE.0.) TYPEBC = 2 IF (BVALUS(3).NE.0.) TYPEBC = 3 IF (TYPEBC.EQ.0) WRITE (MOUTPT,9001) C C IF THE TYPE OF B.C AT THE NODE I IS THE SAME WITH C THE TYPE OF B.C AT THE RIGHT SUBINTERVAL THEN IT C USES THE INTERVAL (T(I),T(I+1/2)) TO DETERMINE THE C THE HERMITE INTERPOLANT OTHERWISE THE INTERVAL C (T(I-1/2),T(I)) C IF (BVALUS(TYPEBC).NE.BVLUS1(TYPEBC)) GO TO 10 C USE RIGHT HALF INTERVAL IBLK0 = 26*I - 26 IBLK1 = IBLK0 + 2 IBLK2 = IBLK0 + 14 TAU1R = (T(I)+R(IBLK0+1))*.5 TAU2R = (T(I)+R(IBLK0+2))*.5 RS(1) = VBCND RS(2) = BRHS(ISIDE,VISIDE,TAU1R,BVALUS) RS(3) = BRHS(ISIDE,VISIDE,TAU2R,BVALUS) RS(4) = VTMDLR C DETINV = 1./ (R(IBLK1+3)*R(IBLK2+4)-R(IBLK2+3)*R(IBLK1+4)) IT2 = 2*I Q(IT2-1) = RS(1) B1 = RS(2) - RS(1)*R(IBLK1+1) - RS(4)*R(IBLK1+2) B2 = RS(3) - RS(1)*R(IBLK2+1) - RS(4)*R(IBLK2+2) Q(IT2) = 2.* (B1*R(IBLK2+4)-B2*R(IBLK1+4))*DETINV GO TO 20 C 10 CONTINUE C USE LEFT INTERVAL IBLK0 = 26*I - 52 IBLK1 = IBLK0 + 2 IBLK2 = IBLK0 + 14 TMDLL = (T(I)+T(I-1))*.5 TAU1L = (T(I)+R(IBLK0+1))*.5 TAU2L = (T(I)+R(IBLK0+2))*.5 VBCMDL = BRHS(ISIDE,VISIDE,TMDLL,BVALUS) RS(1) = VBCMDL RS(2) = BRHS(ISIDE,VISIDE,TAU1L,BVALUS) RS(3) = BRHS(ISIDE,VISIDE,TAU2L,BVALUS) RS(4) = VBCND C DETINV = 1./ (R(IBLK1+3)*R(IBLK2+4)-R(IBLK2+3)*R(IBLK1+4)) IT2 = 2*I Q(IT2-1) = RS(4) B1 = RS(2) - RS(1)*R(IBLK1+1) - RS(4)*R(IBLK1+2) B2 = RS(3) - RS(1)*R(IBLK2+1) - RS(4)*R(IBLK2+2) Q(IT2) = 2.* (B2*R(IBLK1+3)-B1*R(IBLK2+3))*DETINV C 20 CONTINUE C 30 CONTINUE C C DETERMINE THE DEGREES OF FREEDOM AT NODE = KNOTS USING C THE INTERVAL (T(KNOTS-1/2),T(KNOTS)) C C CASE NODE = KNOTS C C I = KNOTS TMDLL = (T(I)+T(I-1))*.5 IBLK0 = 26*I - 52 IBLK1 = IBLK0 + 2 IBLK2 = IBLK0 + 14 TAU1L = (T(I)+R(IBLK0+1))*.5 TAU2L = (T(I)+R(IBLK0+2))*.5 VBCMDL = BRHS(ISIDE,VISIDE,TMDLL,BVALUS) RS(1) = VBCMDL RS(2) = BRHS(ISIDE,VISIDE,TAU1L,BVALUS) RS(3) = BRHS(ISIDE,VISIDE,TAU2L,BVALUS) RS(4) = BRHS(ISIDE,VISIDE,T(KNOTS),BVALUS) C DETINV = 1./ (R(IBLK1+3)*R(IBLK2+4)-R(IBLK2+3)*R(IBLK1+4)) C IT2 = 2*KNOTS Q(IT2-1) = RS(4) B1 = RS(2) - RS(1)*R(IBLK1+1) - RS(4)*R(IBLK1+2) B2 = RS(3) - RS(1)*R(IBLK2+1) - RS(4)*R(IBLK2+2) Q(IT2) = 2.* (B2*R(IBLK1+3)-B1*R(IBLK2+3))*DETINV C IF (LEVEL.LT.4) GO TO 50 C DO 40 I = 1,IT2,2 WRITE (MOUTPT,9011) Q(I),Q(I+1) 40 CONTINUE 50 CONTINUE RETURN C 9001 FORMAT (1H ,32HINCOMPATIBLE BOUNDARY CONDITIONS) 9011 FORMAT (1H ,E13.6,3X,E13.6) END FUNCTION BRHS(ISIDE,VISIDE,T,BVALUS) C C DEFINES THE BOUNDARY CONDITIONS OVER THE ENTIRE C REGION C REAL BVALUS(4) C GO TO (10,20,30,40),ISIDE C 10 XBP = VISIDE YBP = T BRHS = BCOND(ISIDE,XBP,YBP,BVALUS) GO TO 50 C 20 XBP = T YBP = VISIDE BRHS = BCOND(ISIDE,XBP,YBP,BVALUS) GO TO 50 C 30 XBP = VISIDE YBP = T BRHS = BCOND(ISIDE,XBP,YBP,BVALUS) GO TO 50 C 40 XBP = T YBP = VISIDE BRHS = BCOND(ISIDE,XBP,YBP,BVALUS) 50 CONTINUE RETURN END SUBROUTINE EQGNHM(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,MXNEQ,RX,NRDX,RY, . NRDY,LEVEL) C **** DISCRETIZATION MODULE COLLOCATION-HOMOGENEOUS **** C C ***** THIS MODULE FORMS THE COLLOCATION EQUATIONS FOR C THE GENERAL SECOND ORDER ELLIPTIC PDE DEFINED C ON THE RECTANGLE (AX,BX)X(AY,BY). THE BOUNDARY C CONDITIONS ARE ASSUMED HOMOGENOUS. THE TYPE OF C BOUNDARY CONDITIONS CAN BE DIRICHLET OR NEUMANN C OR MIXED ON EACH BOUNDARY SIDE. IN CASE OF MIXED C B.C IT IS ASSUMED THAT ON A PART OF EACH BOUNDARY C SIDE - WHOSE END POINTS ARE NODES OF THE DISRCETIZED C REGION - THE BOUNDARY CONDITIONS CAN BE DIRICHLET C OR NEUMANN TYPE.***** C C ***** DECLARATIONS ***** C REAL GRIDX(NGRIDX),GRIDY(NGRIDY),COEF(MXNEQ,*),RX(NRDX),RY(NRDY) C C C ***** DEFINITION OF CONSTANTS AND LOCAL VARIABLES ***** C NX = NGRIDX - 1 NY = NGRIDY - 1 NROW = 1 C C ***** GENERATE COLLOCATION EQUATIONS CORRESPONDING TO C INTERIOR COLLOCATION POINTS ***** C DO 20 IX = 1,NX C IXBLK = 26*IX - 25 C DO 10 JY = 1,NY C JYBLK = 26*JY - 25 + NRDX C C ***** GENERATE A BLOCK OF FOUR EQUATIONS CORRESPONDING TO C COLLOCATION POINTS IN THE ELEMENT (IX,JY) ***** C C C CALL EQBLCK(NROW,RX(IXBLK),RX(IXBLK+2),RX(IXBLK+14), . RY(JYBLK),RY(JYBLK+2),RY(JYBLK+14),COEF,MXNEQ, . MXNCOE) C 10 CONTINUE 20 CONTINUE RETURN END C C SUBROUTINE HRMCOL(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ, C . UNKVCT,NUMUNK,WRKSP,HOMOBC,BCTYPE,P1,P2,LEVEL) C C DUMMY TO USE IN DIRECTORY OF INTERIOR COLLOCATION C C RETURN C END C SUBROUTINE INTCOL(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ, . UNKVCT,NUMUNK,WRKSP,HOMOBC,BCTYPE,LEVEL) C C -- PURPOSE -- C C APPLY COLLOCATION METHOD BASED ON HERMITE BICUBICS C TO DISCRETIZE THE PDE PROBLEM (1),(2) IN CASE B.C C ARE UNCOUPLED OR NEUMANN TYPE. C C -- DESCRIPTION -- C C GENERATE BOUNDARY AND INTERIOR COLLOCATION EQUATIONS, C SOLVES THE UNCOUPLED BOUNDARY EQUATIONS WITH RESPECT C TO SOME DOFS, ELIMINATES THEM FROM INTERIOR EQUATIONS, C NUMBERS THE ACTIVE DOFS, AND FOR EACH INTERIOR EQUA- C TION STORES THE UNKNOWN INDECIES WHOSE COEFFICIENTS C ARE NONZERO. C -- REFERENCES -- C C E.N.HOUSTIS,W.F.MITCHELL AND J.R.RICE,COLLOCATION C SOFTWARE FOR SECOND ORDER ELLIPTIC PARTIAL DIFFERENTIAL C TRANSACTIONS OF MATH. SOFT. TO APPEAR. C C -- DECLARATIONS -- C REAL GRIDX(*),GRIDY(*),COEF(MXNEQ,*),UNKVCT(4,*),WRKSP(*) INTEGER IDCOEF(MXNEQ,*),NUMUNK(4,*),LEVEL CHARACTER *4 BCTYPE LOGICAL HOMOBC C C -- INPUT ARGUMENTS -- C C GRIDX,NGRIDX,GRIDY,NGRIDY,BCTYPE,HOMOBC C (DEFINED IN S/R RCTCOL) C C -- OUTPUT ARGUMENTS -- C C COEF,IDCOEF,UNKVCT,NUMUNK,WRKSP C (DEFINED IN S/R RCTCOL) C C C /* ALGORITHM INTERIOR COLLOCATION */ C C PROCESS P1: GENERATE DETAIL GRID INFORMATION C C *NOTE* IN THIS CASE THE USER SUPPLIES THIS INFORMATION C C PROCESS P2: FIND INTERIOR X - COLLOCATION POINTS, EVALUATE C BASIS FUNCTIONS AND DERIVATIVES AT THESE POINTS C NRDX = 26* (NGRIDX-1) C CALL INTPTS(GRIDX,NGRIDX,WRKSP(1),NRDX,LEVEL) C C PROCESS P3: FIND INTERIOR Y - COLLOCATION POINTS, EVALUATE C BASIS FUNCTIONS AND DERIVATIVES AT THESE POINTS C NRDY = 26* (NGRIDY-1) LRS = NRDX + 1 C CALL INTPTS(GRIDY,NGRIDY,WRKSP(LRS),NRDY,LEVEL) C C PROCESS P4: LOOP OVER ELEMENTS , K=1 TO NX*NY , OF DOMAIN C AND COMPUTE THE COEFFICIENTS FROM COLLOCATION C OF PDE AT 4 INTERIOR POINTS OF EACH ELEMENT C NUMBEQ = 4* (NGRIDX-1)* (NGRIDY-1) NUMCOE = 17 C CALL EQGNHM(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,MXNEQ,WRKSP(1),NRDX, . WRKSP(LRS),NRDY,LEVEL) C C PROCESS P5: C STEP1: DETERMINES TWO DEGREES OF FREEDOM AT EACH C BOUNDARY NODE FROM THE HERMITE INTERPOLANT C OF 'G' ON EACH BOUNDARY SIDE C STEP2: ENFORCES UNIQUENESS ,IN CASE OF NEUMANN C BOUNDARY CONDITIONS, USING THE INFORMATION C FROM USER SUPPLIED FUNCTION 'UNQU' C STEP3: INDEX THE UNDETERMINED UNKNOWNS BOTTOM TO C TOP , THEN LEFT TO RIGHT C STEP4: ELIMINATES THE COMPUTED BOUNDARY DEGREES C OF FREEDOM FROM THE INTERIOR COLLOCATION C EQUATIONS BY UPDATING THE RIGHT SIDE C LRS1 = NRDX + NRDY + 1 NQD = 2*MAX0(NGRIDX,NGRIDY) C CALL ORDRHM(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ,UNKVCT, . NUMUNK,WRKSP(1),WRKSP(LRS),WRKSP(LRS1),HOMOBC,BCTYPE, . LEVEL) RETURN END SUBROUTINE ORDRHM(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ, . UNKVCT,NUMUNK,RX,RY,Q,HOMOBC,BCTYPE,LEVEL) C C -- PURPOSE -- C C GENERATE BOUNDARY COLLOCATION EQUATIONS, SOLVES C THEM WITH RESPECT TO CERTAIN DOFS, ELIMINATES THEM C FROM INTERIOR COLLOCATION EQUATIONS, NUMBERS THE C REST OF DOFS AND FOR EACH INTERIOR EQUATION C STORES THE INDECIES OF ACTIVE DOFS WHOSE C COEFFICIENTS ARE NONZERO. C C -- DESCRIPTION -- C C ORDRHM APPLIES THE FOLLOWING 4 STEPS : C STEP1: DETERMINES TWO DEGREES OF FREEDOM AT EACH C BOUNDARY NODE FROM THE HERMITE INTERPOLANT C OF 'G' ON EACH BOUNDARY SIDE C STEP2: ENFORCES UNIQUENESS ,IN CASE OF NEUMANN C BOUNDARY CONDITIONS, USING THE INFORMATION C FROM USER SUPPLIED FUNCTION 'UNQU' C STEP3: INDEX THE UNDETERMINED UNKNOWNS BOTTOM TO C TOP , THEN LEFT TO RIGHT C STEP4: ELIMINATES THE COMPUTED BOUNDARY DEGREES C OF FREEDOM FROM THE INTERIOR COLLOCATION C EQUATIONS BY UPDATING THE RIGHT SIDE C C C C -- INPUT ARGUMENTS -- C C COEF COEFÕ1:MXNEQ;1:17å:COEFFICIENTS C OF COLLOCATION EQUATIONS C GRIDX,GRIDY GRIDXÕ1:NGRIDXå,GRIDYÕ1:NGRIDYå: C GRID POINT COORDINATES C RX,RY RXÕ1:NRDXå,RYÕ1:NRDYå: INTERIOR C COLLOCATION POINTS AND C VALUES OF THE BASIS FUNCTIONS C AND THEIR DERIVATIVES C (DEFINED IN S/R INTPTS) C LEVEL LEVEL: OUTPUT LEVEL C HOMOBC HOMOBC : LOGICAL VARIABLE C WHICH DEFINES CHARACTERISTICS C OF THE PDE PROBLEM C BCTYPE BCTYPE : CHARACTER VARIABLE C DEFINES THE TYPE OF B.C C C -- OUTPUT ARGUMENTS-- C COEF COEFÕNROW,17å: THE RIGHT HAND SIDE C HAS BEEN MODIFIED AFTER THE C ELIMINATION OF BOUNDARY DOF C THAT MAKE 'UC' TO SATISFY C EXACTLY THE HERMITE INTERPOLANT C OF 'U' AT THE BOUNDARY C IDCOEF IDCOEFÕNROW,1:16å:THE INDECIES OF C THE UNKNOWNS FOR EQUATION 'NROW' C IF IDCOEFÕNROW,.å IS 0 THAT IMPLIES C THAT THE CORRESPONDING UNKNOWN HAS C BEEN ELIMINATED. C NUMUNK NUMUNKÕ1:4,NODå: IF ZERO IMPLIES THAT C THE CORRESPONDING DOF IS DETERMINED C BY INTERPOLATION. C UNKVCT UNKVCTÕ1:4,NODå: THE VALUES OF DEGREES C OF FREEDOM AT NODE 'NOD'.AT THIS C STAGE ARE ZERO BUT THOSE THAT HAVE C PRECOMPUTED BY INTERPOLATING THE B.C C ------ DECLARATIONS ------- C C REAL GRIDX(*),GRIDY(*),BVALUS(4),COEF(MXNEQ,*),RX(*),RY(*),Q(*), . UNKVCT(4,*) INTEGER IDCOEF(MXNEQ,*),NUMUNK(4,*),UL,UR CHARACTER *4 BCTYPE,ISYMBL LOGICAL HOMOBC DATA ISYMBL/'NEUM'/ C C C -------- TEMPORARY VARIABLES -------- C NX = NGRIDX - 1 NY = NGRIDY - 1 NQ3D = 2*NGRIDY NQ1D = NQ3D NQ4D = 2*NGRIDX NQ2D = NQ4D NRYD = 26*NY NRXD = 26*NX NQD = MAX0(NQ3D,NQ4D) NTEMP = NGRIDY* (NGRIDX-1) NNODES = NTEMP + NGRIDY C C INITIALIZATIONS C DO 20 NOD = 1,NNODES DO 10 NDF = 1,4 NUMUNK(NDF,NOD) = 0 UNKVCT(NDF,NOD) = 0. 10 CONTINUE 20 CONTINUE C C DO 30 K = 1,NQD Q(K) = 0. 30 CONTINUE C C APPLIES STEP 1 ON THE FOUR SIDES OF THE RECTANGLE C C CASE ISIDE = 3 C IF (HOMOBC) GO TO 40 CALL BCINTR(3,GRIDX(1),Q,NQ3D,RY,NRYD,GRIDY,NGRIDY) 40 CONTINUE C DO 70 J = 1,NGRIDY NOD = J JT2 = 2*J TEMP = BCOND(3,GRIDX(1),GRIDY(J),BVALUS) IF (BVALUS(1).EQ.0.) GO TO 50 C THE B.C IS DIRICHLET TYPE : DETERMINES UNKNOWNS U,UY NUMUNK(1,NOD) = 1 NUMUNK(2,NOD) = 1 UNKVCT(1,NOD) = Q(JT2-1) UNKVCT(2,NOD) = Q(JT2) C 50 IF (BVALUS(2).EQ.0.) GO TO 60 C THE B.C IS NEUMANN TYPE : DETERMINES UNKNOWNS UX,UXY NUMUNK(3,NOD) = 1 NUMUNK(4,NOD) = 1 UNKVCT(3,NOD) = Q(JT2-1) UNKVCT(4,NOD) = Q(JT2) 60 CONTINUE 70 CONTINUE C C CASE ISIDE = 1 C IF (HOMOBC) GO TO 80 CALL BCINTR(1,GRIDX(NGRIDX),Q,NQ1D,RY,NRYD,GRIDY,NGRIDY) 80 CONTINUE DO 110 J = 1,NGRIDY NOD = J + NTEMP JT2 = 2*J TEMP = BCOND(1,GRIDX(NGRIDX),GRIDY(J),BVALUS) IF (BVALUS(1).EQ.0.) GO TO 90 C THE B.C IS DIRICHLET TYPE : DETERMINES UNKNOWNS U,UY NUMUNK(1,NOD) = 1 NUMUNK(2,NOD) = 1 UNKVCT(1,NOD) = Q(JT2-1) UNKVCT(2,NOD) = Q(JT2) 90 IF (BVALUS(2).EQ.0.) GO TO 100 C THE B.C IS NEUMANN TYPE : DETERMINES UNKNOWNS UX,UXY NUMUNK(3,NOD) = 1 NUMUNK(4,NOD) = 1 UNKVCT(3,NOD) = Q(JT2-1) UNKVCT(4,NOD) = Q(JT2) 100 CONTINUE 110 CONTINUE C C CASE ISIDE = 4 C IF (HOMOBC) GO TO 120 CALL BCINTR(4,GRIDY(NGRIDY),Q,NQ4D,RX,NRXD,GRIDX,NGRIDX) 120 CONTINUE C DO 150 J = 1,NGRIDX NOD = J*NGRIDY JT2 = 2*J TEMP = BCOND(4,GRIDX(J),GRIDY(NGRIDY),BVALUS) IF (BVALUS(1).EQ.0.) GO TO 130 C THE B.C IS DIRICHLET TYPE : DETERMINES UNKNOWNS U,UX NUMUNK(1,NOD) = 1 NUMUNK(3,NOD) = 1 UNKVCT(1,NOD) = Q(JT2-1) UNKVCT(3,NOD) = Q(JT2) 130 IF (BVALUS(3).EQ.0.) GO TO 140 C THE B.C IS NEUMANN TYPE : DETERMINES UNKNOWNS UY,UXY NUMUNK(2,NOD) = 1 NUMUNK(4,NOD) = 1 UNKVCT(2,NOD) = Q(JT2-1) UNKVCT(4,NOD) = Q(JT2) 140 CONTINUE 150 CONTINUE C C 11 CASE ISIDE = 2 C IF (HOMOBC) GO TO 160 CALL BCINTR(2,GRIDY(1),Q,NQ2D,RX,NRXD,GRIDX,NGRIDX) 160 CONTINUE DO 190 J = 1,NGRIDX TEMP = BCOND(2,GRIDX(J),GRIDY(1),BVALUS) JT2 = 2*J NOD = J*NGRIDY - NGRIDY + 1 IF (BVALUS(1).EQ.0.) GO TO 170 NUMUNK(1,NOD) = 1 NUMUNK(3,NOD) = 1 UNKVCT(1,NOD) = Q(JT2-1) UNKVCT(3,NOD) = Q(JT2) 170 IF (BVALUS(3).EQ.0.) GO TO 180 NUMUNK(2,NOD) = 1 NUMUNK(4,NOD) = 1 UNKVCT(2,NOD) = Q(JT2-1) UNKVCT(4,NOD) = Q(JT2) 180 CONTINUE 190 CONTINUE C C APPLIES STEP 2 C C IN CASE OF NEUMANN CONDITIONS OVER THE ENTIRE BOUNDARY C SETS U AT (AX,BX) WITH THE SUPPLIED VALUE AT THIS POINT C C IF (BCTYPE.NE.ISYMBL) GO TO 200 NUMUNK(1,1) = 1 UNKVCT(1,1) = UNQU(GRIDX(1),GRIDY(1)) NUMUNK(3,1) = 0 UNKVCT(3,1) = 0. 200 CONTINUE C C NUMBER THE UNDETERMINED DEGREES OF FREEDOM C NUN = 0 DO 230 NOD = 1,NNODES DO 220 NDF = 1,4 IF (NUMUNK(NDF,NOD).EQ.1) GO TO 210 NUN = NUN + 1 NUMUNK(NDF,NOD) = NUN GO TO 220 C 210 NUMUNK(NDF,NOD) = 0 220 CONTINUE 230 CONTINUE C C APPLIES STEP 4 , 5 C C FORM ARRAY " IDCOEF" , ELIMINATES THE PREDETERMINED BOUNDARY C CONDITIONS BY UPDATING THE RIGHT SIDE C NROW = 0 C DO 270 IX = 1,NX DO 260 JY = 1,NY C FIND ELEMENT INCIDENCIES LL = JY + (IX-1)*NGRIDY LR = LL + NGRIDY UR = LR + 1 UL = LL + 1 C DO 250 K = 1,4 C INCREASE EQUATION COUNTER C NROW = NROW + 1 DO 240 NDF = 1,4 C DETERMINES IDCOEF IDCOEF(NROW,NDF) = NUMUNK(NDF,LL) C UPDATES RIGHT SIDE OF INTERIOR EQUATIONS IF (NUMUNK(NDF,LL).EQ.0) COEF(NROW,17) = COEF(NROW, . 17) - COEF(NROW,NDF)*UNKVCT(NDF,LL) C NDF1 = NDF + 4 C DETERMINES IDCOEF IDCOEF(NROW,NDF1) = NUMUNK(NDF,LR) C UPDATES RIGHT SIDE OF INTERIOR EQUATIONS IF (NUMUNK(NDF,LR).EQ.0) COEF(NROW,17) = COEF(NROW, . 17) - COEF(NROW,NDF1)*UNKVCT(NDF,LR) C NDF2 = NDF1 + 4 C DETERMINES IDCOEF IDCOEF(NROW,NDF2) = NUMUNK(NDF,UR) C UPDATES RIGHT SIDE OF INTERIOR EQUATIONS IF (NUMUNK(NDF,UR).EQ.0) COEF(NROW,17) = COEF(NROW, . 17) - COEF(NROW,NDF2)*UNKVCT(NDF,UR) C NDF3 = NDF2 + 4 C DETERMINES IDCOEF IDCOEF(NROW,NDF3) = NUMUNK(NDF,UL) C UPDATES RIGHT SIDE OF INTERIOR EQUATIONS IF (NUMUNK(NDF,UL).EQ.0) COEF(NROW,17) = COEF(NROW, . 17) - COEF(NROW,NDF3)*UNKVCT(NDF,UL) C 240 CONTINUE 250 CONTINUE 260 CONTINUE 270 CONTINUE RETURN END C C//////////////////////////////////////////////////////////////// C**************** ALGORITHMS INTCOL & HERMCOL ******************* C/////////////////////////////////////////////////////////////// C>>>>>>>>>>>>> LOGICAL FILE 3 : RECTANGULAR UTILITIES <<<<<<<<<< C//////////////////////////////////////////////////////////////// C SUBROUTINE EQBLCK(NROW,X,BXC1,BXC2,Y,BYC1,BYC2,COEF,MXNEQ,MXNCOE) C C ***** GENERATES FOUR ALGEBRAIC EQUATIONS BY FORCING C THE COLLOCATION APPROXIMATION ( A HERMITE BICUBIC C PIECEWISE POLYNOMIAL ) 'UC(X,Y)' TO SATISFY THE C THE ELLIPTIC EQUATION AT THE FOUR GAUSS POINTS OF A C RECTANGULAR ELEMENT OF THE PARTITION OF THE REGION ***** C C ***** COLLOCATION APPROXIMATION : C C UC(X,Y) = Q( 1)BX1*BY1 + Q( 2)BX1*BY3 + Q( 3)BX3BY1 + Q( 4)BX3BY3 C C + Q( 5)BX2*BY1 + Q( 6)BX2*BY3 + Q( 7)BX4BY1 + Q( 8)BX4BY3 C C + Q( 9)BX2*BY2 + Q(10)BX2*BY4 + Q(11)BX4BY2 + Q(12)BX4BY4 C C + Q(13)BX1*BY2 + Q(14)BX1*BY4 + Q(15)BX3BY2 + Q(16)BX3BY4 C C--LOCAL DEFINITION OF HERMITE BASIS FUNCTIONS IN X-DIRECTION C ( SIMILARLY ONE CAN DEFINE THE BASIS FUNCTIONS IN Y-DIRECTION ) C C C BX1 I++ ++ ** ** I BX2 I BX3 C I ++ ** I I + + C I ** ++ I I + + C I** ++ I I+ + C ------X------------------X----- ----X----------------X------ C X(I) X(I+1) X(I) + + X(I+1) C + + C BX4 + + C--NUMBERING OF TWO-DIMENSIONAL BASIS FUNCTIONS OR DEGREES OF C FREEDOM (DOF) -- C C HERMITE ELEMENT C NODE: UL ------------------- NODE :UR , C DOF :Q(I),I=13,16 I I DOF : Q(I),I=9,12 C I I C I I C I I C NODE: LL ------------------- NODE :LR , C DOF :Q(I),I=1,4 DOF : Q(I),I=5,8 C C--NATURAL RELATION OF DOF'S AND UC-- C C Q(1) : UC(X(LL),Y(LL)), Q(2) : DYUC(X(LL),Y(LL)) C C Q(3) : DXUC(X(LL),Y(LL)), Q(4) : DXYUC(X(LL),Y(LL)) C C NOTE : THE SAME RELATION HOLDS FOR OTHER DEGREES OF FREEDOM C C--COLLOCATION EQUATIONS -- C C Õ CUXX*DXXUC + CUXY*DXYUC + CUYY*DYYUC + CUX*DXUC + CUY*DYU + CU*U C = F å (.,.) C ************************************************* C C ------ INPUT ARGUMENTS :: NROW : INDEX OF THE NEXT COLLOCATION C EQUATION, C X,YÕ1:2å : ARRAYS WITH THE X AND Y C COORDINATES OF THE FOUR GAUSS C POINTS, C BXC1,BXC2Õ1:12å : THE VALUES OF THE BASIS C FUNCTIONS AND THEIR DERIVATIVES C AT POINTS X(1),X(2), C BYC1,BYC2Õ1:12å : THE VALUES OF THE BASIS C FUNCTIONS AND THEIR DERIVATIVES C AT POINTS Y(1),Y(2), C C ------ OUTPUT ARGUMENTS :: COEFÕNROW:NROW+3,1:17å : COEFFICIENTS C AND RIGHT SIDE C C ***** DECLARATIONS ***** C REAL X(2),BX1(12),BXC1(12),BXC2(12),Y(2),BY1(12),BYC1(12), . BYC2(12),COEF(MXNEQ,*),CVALUS(7) C C LOOP OVER 4 COLLOCATION POINTS (X(I),Y(J) TO C COMPUTE COEFFICIENTS FROM COLLOCATION EQUATIONS C DO 140 K = 1,4 C C SET COLLOCATION POINT INDECES I,J AND PUT VALUES C OF BASIS FUNCTIONS AND THEIR DERIVATIVES INTO C APPROPRIATE TEMPORARY ARRAYS BX1,BY1 C GO TO (10,30,70,50),K C 10 I = 1 J = 1 DO 20 L = 1,12 BY1(L) = BYC1(L) BX1(L) = BXC1(L) 20 CONTINUE GO TO 90 C 30 J = 1 I = 2 DO 40 L = 1,12 BX1(L) = BXC2(L) 40 CONTINUE GO TO 90 C 50 J = 2 I = 1 DO 60 L = 1,12 BX1(L) = BXC1(L) 60 CONTINUE GO TO 90 C 70 J = 2 I = 2 DO 80 L = 1,12 BY1(L) = BYC2(L) 80 CONTINUE 90 CONTINUE C C C EVALUATE PDE COEFFICIENTS AT (X(I),Y(J)) C CALL PDE(X(I),Y(J),CVALUS) C TEMPORARY VARIABLES T11 = CVALUS(1)*BX1(9) T12 = CVALUS(1)*BX1(10) T13 = CVALUS(1)*BX1(11) T14 = CVALUS(1)*BX1(12) C IF THE COEFFICIENT OF UX = 0 SKIP RELATED COMPUTATIONS IF (CVALUS(4).EQ.0.0) GO TO 100 C T11 = T11 + CVALUS(4)*BX1(5) T12 = T12 + CVALUS(4)*BX1(6) T13 = T13 + CVALUS(4)*BX1(7) T14 = T14 + CVALUS(4)*BX1(8) C 100 CONTINUE C IF THE COEFFICIENT OF U = 0 SKIP RELATED COMPUTATIONS IF (CVALUS(6).EQ.0.0) GO TO 110 C T11 = T11 + CVALUS(6)*BX1(1) T12 = T12 + CVALUS(6)*BX1(2) T13 = T13 + CVALUS(6)*BX1(3) T14 = T14 + CVALUS(6)*BX1(4) C 110 CONTINUE T21 = CVALUS(3)*BY1(9) T22 = CVALUS(3)*BY1(10) T23 = CVALUS(3)*BY1(11) T24 = CVALUS(3)*BY1(12) C IF THE COEFFICIENT OF UY = 0 SKIP RELATED COMPUTATIONS IF (CVALUS(5).EQ.0.0) GO TO 120 C T21 = T21 + CVALUS(5)*BY1(5) T22 = T22 + CVALUS(5)*BY1(6) T23 = T23 + CVALUS(5)*BY1(7) T24 = T24 + CVALUS(5)*BY1(8) 120 CONTINUE C C COMPUTE COEFFICIENTS AND RIGHT SIDE OF THE COLLOCATION C EQUATIONS FOR THE OPERATOR C C L(UC) = CUXX*DXXUC+CUYY*DYYUC+CUX*DXUC+CUY*DYUC+CU*UC C COEF(NROW,1) = BY1(1)*T11 + BX1(1)*T21 COEF(NROW,2) = BY1(3)*T11 + BX1(1)*T23 COEF(NROW,3) = BY1(1)*T13 + BX1(3)*T21 COEF(NROW,4) = BY1(3)*T13 + BX1(3)*T23 C COEF(NROW,5) = BY1(1)*T12 + BX1(2)*T21 COEF(NROW,6) = BY1(3)*T12 + BX1(2)*T23 COEF(NROW,7) = BY1(1)*T14 + BX1(4)*T21 COEF(NROW,8) = BY1(3)*T14 + BX1(4)*T23 C COEF(NROW,9) = BY1(2)*T12 + BX1(2)*T22 COEF(NROW,10) = BY1(4)*T12 + BX1(2)*T24 COEF(NROW,11) = BY1(2)*T14 + BX1(4)*T22 COEF(NROW,12) = BY1(4)*T14 + BX1(4)*T24 C COEF(NROW,13) = BY1(2)*T11 + BX1(1)*T22 COEF(NROW,14) = BY1(4)*T11 + BX1(1)*T24 COEF(NROW,15) = BY1(2)*T13 + BX1(3)*T22 COEF(NROW,16) = BY1(4)*T13 + BX1(3)*T24 C IF THE COEFFICIENT OF UXY IS 0 SKIP RELATED COMPUTATIONS IF (CVALUS(2).EQ.0.0) GO TO 130 C T31 = CVALUS(2)*BX1(5) T32 = CVALUS(2)*BX1(6) T33 = CVALUS(2)*BX1(7) T34 = CVALUS(2)*BX1(8) C C C COMPUTE COEFFICIENTS AND RIGHT SIDE OF THE COLLOCATION C EQUATIONS CORRESPONDING TO THE OPERATOR L(UC) = L(UC) + CUXY*UXY C COEF(NROW,1) = BY1(5)*T31 + COEF(NROW,1) COEF(NROW,2) = BY1(7)*T31 + COEF(NROW,2) COEF(NROW,3) = BY1(5)*T33 + COEF(NROW,3) COEF(NROW,4) = BY1(7)*T33 + COEF(NROW,4) C COEF(NROW,5) = BY1(5)*T32 + COEF(NROW,5) COEF(NROW,6) = BY1(7)*T32 + COEF(NROW,6) COEF(NROW,7) = BY1(5)*T34 + COEF(NROW,7) COEF(NROW,8) = BY1(7)*T34 + COEF(NROW,8) C COEF(NROW,9) = BY1(6)*T32 + COEF(NROW,9) COEF(NROW,10) = BY1(8)*T32 + COEF(NROW,10) COEF(NROW,11) = BY1(6)*T34 + COEF(NROW,11) COEF(NROW,12) = BY1(8)*T34 + COEF(NROW,12) C COEF(NROW,13) = BY1(6)*T31 + COEF(NROW,13) COEF(NROW,14) = BY1(8)*T31 + COEF(NROW,14) COEF(NROW,15) = BY1(6)*T33 + COEF(NROW,15) COEF(NROW,16) = BY1(8)*T33 + COEF(NROW,16) 130 CONTINUE C ***** COMPUTE THE RIGHT SIDE ***** C COEF(NROW,17) = PDERHS(X(I),Y(J)) C INCREASE EQUATION COUNTER NROW = NROW + 1 140 CONTINUE RETURN END SUBROUTINE HBASVD(BSVD,Z,H) C C ***** EVALUATE THE NON ZERO BASIS FUNCTIONS OF THE 1-D HERMITE CUBIC C PIECEWISE POLYNOMIAL AND THEIR FIRST AND SECOND DERIVATIVES C AT A POINT 'XPOINT' IN THE SUBINTERVAL (X(RIGHT),X(LEFT)) OF C LENGTH H ***** C C --- INPUT ARGUMENTS --- H: SIZE OF THE INTERVAL (X(RIGHT),X(LEFT)), C Z: Z = XPOINT - X(LEFT); C --- OUTPUT ARGUMENTS --- BSVDÕ1:12å : VALUES OF BASIS FUNCTIONS AND C AND THEIR DERIVATIVES AT 'XPOINT'; C --- DECLARATIONS --- REAL BSVD(12) C C --- LOCAL VARIABLES --- C RH = 1./H S = Z*RH T = 1. - S STT = S*T ST2 = S*2. TT2 = T*2. TEMP = S - TT2 C C --- BASIS FUNCTIONS --- C BSVD(1) = S*S* (ST2-3.) + 1. BSVD(2) = 1. - BSVD(1) BSVD(3) = STT* (H-Z) BSVD(4) = -STT*Z C C --- 1-ST DERIVATIVE OF BASIS FUNCTIONS --- C BSVD(5) = 6.* (S-1.)*RH*S BSVD(6) = -BSVD(5) BSVD(7) = T* (T-ST2) BSVD(8) = S*TEMP C C --- 2-OD DERIVATIVE OF BASIS FUNCTIONS --- C BSVD(9) = 6.* (ST2-1.)*RH*RH BSVD(10) = -BSVD(9) BSVD(11) = 2.*TEMP*RH BSVD(12) = 6.* (S-T)*RH - BSVD(11) RETURN END SUBROUTINE INTPTS(T,KNOTS,R,NRD,LEVEL) C C **** GIVEN A PARTITION TÕ1:KNOTSå OF AN INTERVAL ÕA,Bå C (I) DETERMINES THE TWO GAUSS POINTS OF ÕT(I),T(I+1)å C (II) EVALUATES THE BASIS FUNCTIONS AND THEIR 1-ST AND C 2-ND DERIVATIVES AT THE ABOVE POINTS ***** C C ---- INPUT ARGUMENTS :: TÕ1:KNOTSå :ARRAY OF PARTITION POINTS, C KNOTS : NUMBER OF PARTITION POINTS, C LEVEL : OUTPUT LEVEL ; ----- C C ---- OUTPUT ARGUMENTS :: RÕ1:NRDå: R(1),R(2) ARE THE GAUSS POINTS C OF THE INTERVAL ÕT(1),T(2)å,R(3)-R(6) C VALUES OF BASIS FUNCTIONS AT R(1), C R(7)-R(10) VALUES OF THE FIRST C DERIVATIVE AT R(1), R(11)-R(14) VALUES C OF THE SECOND DERIVATIVE AT C R(1),R(15)-R(26) VALUES OF BASIS C FUNCTIONS AND THEIR DERIVATIVES AT C R(2),..., NRD : DIMENSION OF R IS C 26*(KNOTS-1); --------- C ----- DECLARATIONS ------------------ C REAL T(KNOTS),R(NRD) C C ----- CONSTANTS ------ C SQRT3 = SQRT(3.) FACTOR = .5/SQRT3 C ----- LOCAL VARIABLE ----- C NINT = KNOTS - 1 C DO 10 I = 1,NINT C C COMPUTES FIRST GAUSS POINT AND STORES AT R(IBLOK1) C IP1 = I + 1 H = (T(IP1)-T(I)) TMDL = (T(IP1)+T(I))*.5 FCTR = H*FACTOR C IBLOK1 = 26*I - 25 R(IBLOK1) = TMDL - FCTR TAU = R(IBLOK1) - T(I) C C COMPUTES BASIS FUNCTIONS AND THEIR DERIVATIVES AT R(IBLOK1) AND C STORES THEM AT R(IBLOK3),...,R(IBLOKOCK3+12) C IBLOK3 = IBLOK1 + 2 CALL HBASVD(R(IBLOK3),TAU,H) C C COMPUTES SECOND GAUSS POINT AND STORES AT R(IBLOK2) C IBLOK2 = IBLOK1 + 1 R(IBLOK2) = TMDL + FCTR TAU = R(IBLOK2) - T(I) C C COMPUTES BASIS FUNCTIONS AND THEIR DERIVATIVES AT R(IBLOK2) AND C STORES THEM AT R(IBLOK4),...,R(IBLOK4+12) C IBLOK4 = IBLOK1 + 14 CALL HBASVD(R(IBLOK4),TAU,H) C 10 CONTINUE C C IF LEVEL = 4 PRINTS ARRAY RÕ1:NRDå C IF (LEVEL.LT.4) GO TO 20 WRITE (6,9001) WRITE (6,9011) (R(IR),IR=1,NRD) 20 CONTINUE RETURN C 9001 FORMAT (39H PRINTS 1-D BASIS , 1-ST AND 2-ND DERIV) 9011 FORMAT (1X,8F10.4) END SUBROUTINE RCTCOL(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ, . UNKVCT,NUMUNK,WRKSP,HOMOBC,BCTYPE,P1,P2,OUTFNC, . OUTTYP,NOUT,TABX,TABY,NTABX,NTABY) C C /* PROGRAM DESCRIPTION */ C C PROGRAM 'RCTCOL' IS THE REALIZATION OF THE COLLOCATION C METHOD DESCRIBED IN ÕHOUSTIS, MITCHELL, RICE (1983), C "COLLOCATION SOFTWARE FOR SECOND ORDER ELLIPTIC PARTIAL C DIFFERENTIAL EQUATIONS"å FOR APPROXIMATING THE SOLUTION C OF THE SECOND ORDER 2-DIMENSIONAL ELLIPTIC PDE : C C (1) CUXX*DDXU+CUXY*DXDYU+CUYY*DDYU+CUX*DXU+CUY*DYU+CU*U = F C C DEFINED ON A RECTANGULAR REGION ÕAX,BXåXÕAY,BYå C AND SUBJECT TO BOUNDARY CONDITIONS C C (2) AU+BUX+CUY = G ON THE BOUNDARY OF THE REGION C C /* NOTE 1 */ C THERE ARE TWO DIFFERENT IMPLEMENTATIONS OF THE C DISCRETIZATION OF THE PDE PROBLEM WHICH DEPEND ON C THE NATURE OF THE BOUNDARY CONDITIONS. C C CASE BOUND-COND-TYPE OF : C /* A&B#0 OR A&C#0 MIGHT HOLD FOR SOME (X,Y) */ C MIXED : APPLY 'HERMITE COLLOCATION' DISCRETIZATION C ÕALGORITHM HERMCOLå C /* ALL BUT ONE OF A,B,C IS ZERO AT EACH (X,Y) */ C UNCOUPLED : APPLY 'INTERIOR COLLOCATION' DISCRETIZATION C ÕALGORITHM INTCOLå C /* B OR C # 0 AT EACH (X,Y) */ C NEUMANN : APPLY 'INTERIOR COLLOCATION' C C C /* NOTE 2 */ C THE 'INTERIOR COLLOCATION' DISCRETIZATION EXPECTS THAT C (I) THE BOUNDARY CONDITIONS ARE DEFINED BETWEEN BOUNDARY C NODES C (II) IN CASE OF 'NEUMAN' BOUNDARY CONDITIONS THE VALUE C OF THE PDE SOLUTION IS DEFINED AT A BOUNDARY NODE C BY THE USER SUPPLIED FUNCTION 'UNQU(UNQX,UNQY)'. C C C /* PROBLEM DEFINITION */ C (A) USER SUPPLIED FUNCTIONS FOR THE COEFFICIENTS AND RIGHT C SIDE OF THE PDE (1) C C REAL FUNCTION PDE(X,Y) C REAL CVALUS(7) C CVALUS(1) = CUXX(X,Y) C CVALUS(2) = CUXY(X,Y) C CVALUS(3) = CUYY(X,Y) C CVALUS(4) = CUX(X,Y) C CVALUS(5) = CUY(X,Y) C CVALUS(6) = CU(X,Y) C CVALUS(7) = F(X,Y) C RETURN C END C C(B)USER SUPPLIED FUNCTION FOR THE BOUNDARY CONDITIONS (2) C C REAL FUNCTION BCOND(ISIDE,X,Y,BVALUS) C VALUES OF THE BOUNDARY CONDITION ON SIDE I C AT (X,Y) C 9 C I = 4 C ----------------- C I I C I I C I=3 I I I=1 C I I C I I C ----------------- C I = 2 C REAL BVALUS(4) C GO TO (10,20,30,40),ISIDE C10 BVALUS(1) = A C BVALUS(2) = B C BVALUS(3) = C C BVALUS(4) = G(X,Y) C BCOND = BVALUS(4) C GO TO 100 C20 BVALUS(1) = A C BVALUS(2) = B C BVALUS(3) = C C BVALUS(4) = G(X,Y) C BCOND = BVALUS(4) C GO TO 100 C30 BVALUS(1) = A C BVALUS(2) = B C BVALUS(3) = C C BVALUS(4) = G(X,Y) C BCOND = BVALUS(4) C GO TO 100 C40 BVALUS(1) = A C BVALUS(2) = B C BVALUS(3) = C C BVALUS(4) = G(X,Y) C BCOND = BVALUS(4) C RETURN C END C(C)USER SUPPLIED FUNCTION FOR THE TRUE SOLUTION IF KNOWN C C REAL FUNCTION TRUE(X,Y) C TRUE = ... C RETURN C END C C(D) USER SUPPLIED FUNCTION FOR OBTAINING A UNIQUE SOLUTION C OF THE PDE PROBLEM WITH NEUMANN BOUNDARY CONDITIONS. C IT IS USED ONLY BY 'INTERIOR COLLOCATION' DISCRETIZATION C MODULE. C C REAL FUNCTION UNQU(UNQX,UNQY) C UNQX = /* THE X-COORDINATE OF THE BOUNDARY NODE */ C UNQY = /* THE Y-COORDINATE OF THE BOUNDARY NODE */ C UNQU = /* THE VALUE OF 'U' AT (UNQX,UNQY) */ C RETURN C END C C C /* OUTPUT CONTROL */ C C........ C C ********************************************************* C /* COLLOCATION ALGORITHM */ C *********************************************************** C C ----- DEFINITION ------ C ALGORITHM RCTCOL (GRIDX,NGRIDX,GRIDY,NGRIDY,COEF, C A IDCOEF, MXNEQ,UNKVCT,NUMUNK,WRKSP, C B HOMOBC,BCTYPE,LEVEL) C C ------ DECLARATIONS --------- C REAL GRIDX(*),GRIDY(*),COEF(MXNEQ,*),WRKSP(*) REAL UNKVCT(4,*) REAL TABX(*),TABY(*) INTEGER IDCOEF(MXNEQ,*) INTEGER NUMUNK(4,*) INTEGER OUTFNC(*),OUTTYP(*) LOGICAL HOMOBC CHARACTER *4 BCTYPE,IUNCU,INEUM COMMON /INTEGS/NUMBEQ,NUMCOE,LEVEL,MOUTPT,NROW C ------ INITIALIZATIONS ------------- DATA IUNCU,INEUM/'UNCU','NEUM'/ DATA NBANDU,NBANDL/0,0/ C C ------ INPUT ARGUMENTS ---------- C C GRIDX REAL ARRAY GRIDXÕ1:NGRIDXå C GRIDY REAL ARRAY GRIDYÕ1:NGRIDYå C (X,Y) - COORDINATES OF GRID POINTS C BCTYPE CHARACTER VARIABLE C DEFINES THE TYPE OF BOUNDARY C CONDITIONS. THE ALLOWED VALUES ARE C 'MIXD','UNCU','NEUM'. C HOMOBC LOGICAL VARIABLE C INDICATES WHETHER THE RIGHT SIDE OF C BOUNDARY CONDITION IS ZERO OR NOT. C NOUT NUMBER OF OUTPUTS TO PERFORM C OUTFNC(I) WHAT FUNCTION TO USE FOR THE ITH OUTPUT C =1 APPROXIMATE SOLUTION C =2 ERROR = U - TRUE (TRUE MUST BE SUPPLIED) C =3 RESIDUAL = LU - F C OUTTYP(I) WHAT INFORMATION IS TO BE PRINTED C IN THE OUTPUT C =1 MAX, L1, L2 NORMS BASED ON C DISCRETIZATION GRID C =2 MAX, L1, L2 NORMS BASED ON C USER SUPPLIED GRID (TABX,TABY) C =3 TABLE OF FUNCTION ON DISC. GRID C =4 TABLE OF FUNCTION ON USER SUPPLIED GRID C TABX,TABY X AND Y COORDINATES OF GRID FOR OUTTYP 2 AND 4 C NTABX,NTABY NUMBER OF VALUES IN TABX,TABY C C ----- OUTPUT ARGUMENTS ------ C C COEF REAL ARRAY COEFÕ1:MXNEQ;1:17å C EACH ROW CONTAINS THE 16 NONZERO C COEFFICIENTS OF EACH COLLOCATION C EQUATION AND ITS RIDE SIDE. C IDCOEF INTEGER ARRAY IDCOEFÕ1:MXNEQ;1:16å C THE INDECIES OF THE ACTIVE DEGREES OF C FREEDOM (DOFS) ASSOCIATED WITH NONZERO C COEFFICIENTS AT EACH COLLOCATION EQUA- C TION. A GLOBAL NUMBERING OF THE ACTIVE C DOFS IS STORED IN 'NUMUNK' C UNKVCT REAL ARRAY UNKVCTÕ1:4;1:NUMBEQå C THE VALUES OF THE 4*NGRIDX*NGRIDY C DOFS OF THE APPROXIMATE SOLUTION. C NUMUNK INTEGER ARRAY NUMUNKÕ1:4;1:NUMBEQå C A GLOBAL NUMBERING OF ACTIVE DOFS. C THE INDECIES OF INACTIVE DOFS ARE ZERO. C WRKSP REAL ARRAY WRKSPÕ1:NWRKSPå C WORKING SPACE OF LENGTH AT MOST C 32*NGRIDX*(NGRIDY)**2+64*NGRIDX*NGRIDY C FOR 'HRMCOL' OR NWRKSP/2 FOR 'INTCOL'. C C /* DOMAIN DISCRETIZATION */ C C IMPLICITLY DEFINED BY THE USER SUPPLIED ARRAYS C GRIDX(I) , I=1,NGRIDX C GRIDY(J) , J=1,NGRIDY C IT IS ASSUMED THAT AX=GRIDX(1),BX=GRIDX(NGRIDX) C AY=GRIDY(1),BY=GRIDY(NGRIDY) C C NUMBER OF NODES IN THE FINITE ELEMENT PARTITION C NNODES = NGRIDX*NGRIDY C C /* PDE AND BOUNDARY CONDITIONS DISCRETIZATION */ C C SET 'HERMITE COLLOCATION' AS THE DEFUALT OPTION MODULE = 2 C EXAMINE THE TYPE OF BC TO DETERMINE THE APPROPRIATE C DISCRETIZATION MODULE IF (BCTYPE.EQ.IUNCU) MODULE = 1 IF (BCTYPE.EQ.INEUM) MODULE = 1 C CALL THE DISCRETIZATION MODULE MODULE = 2 IF (MODULE.NE.1) GO TO 10 WRITE (MOUTPT,9001) CALL INTCOL(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ,UNKVCT, . NUMUNK,WRKSP,HOMOBC,BCTYPE,LEVEL) NUMBEQ = 4* (NGRIDX-1)* (NGRIDY-1) GO TO 20 10 CONTINUE WRITE (MOUTPT,9011) CALL HRMCOL(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ,UNKVCT, . NUMUNK,WRKSP,HOMOBC,BCTYPE,P1,P2,LEVEL) NUMBEQ = 4*NNODES 20 CONTINUE IF (LEVEL.LT.4) GO TO 40 WRITE (MOUTPT,9021) DO 30 NR = 1,NUMBEQ WRITE (MOUTPT,9031) NR, (COEF(NR,NC),NC=1,16),COEF(NR,17) WRITE (MOUTPT,9041) (IDCOEF(NR,NC),NC=1,16),IDCOEF(NR,17) 30 CONTINUE 40 CONTINUE C C /* REFORMAT LINEAR SYSTEM FOR LINEAR EQUATION SOLVERS OF C BAND MATRIX FORM WITH BANDWIDTH '2*NGRIDY+4' */ C GO TO (50,60),MODULE C 50 CONTINUE C REFORMAT THE COLLOCATION EQUATIONS OF 'INTCOL' NUMCOE = 17 IBLOK1 = 1 LDA = 6*NGRIDY + 10 IBLOK2 = IBLOK1 + LDA*NUMBEQ + 1 MXNCOE = 17 C CALL SETUP(COEF,IDCOEF,MXNEQ,MXNCOE,NUMBEQ,NUMCOE,WRKSP(IBLOK1), . LDA,WRKSP(IBLOK2),NBANDU,NBANDL) GO TO 70 C C REFORMAT THE COLLOCATION EQUATIONS OF 'HRMCOL' C 60 MXNCOE = 17 NUMCOE = 17 LDA = 12*NGRIDY + 22 IBLOK1 = 1 IBLOK2 = IBLOK1 + LDA*NUMBEQ + 1 C CALL SETUP(COEF,IDCOEF,MXNEQ,MXNCOE,NUMBEQ,NUMCOE,WRKSP(1),LDA, . UNKVCT,NBANDU,NBANDL) 70 CONTINUE C C /* SOLVE LINEAR SYSTEM OF COLLOCATION EQUATIONS BY C GAUSS ELIMINATION WITH SCALLING */ C GO TO (80,90),MODULE C 80 CONTINUE IBLOK3 = IBLOK2 + NUMBEQ + 1 CALL SOLVE(WRKSP(1),LDA,NUMBEQ,WRKSP(IBLOK2),WRKSP(IBLOK3),NBANDU, . NBANDL) GO TO 100 C 90 CONTINUE CALL SOLVE(WRKSP(1),LDA,NUMBEQ,UNKVCT,WRKSP(IBLOK2),NBANDU,NBANDL) 100 CONTINUE C C /* DETERMINE COLLOCATION APPROXIMATION */ C IF (MODULE.NE.1) GO TO 130 IDUNK = IBLOK2 - 1 DO 120 NOD = 1,NNODES DO 110 NDF = 1,4 INDX = NUMUNK(NDF,NOD) IF (INDX.EQ.0) GO TO 110 IDUNK = IDUNK + 1 UNKVCT(NDF,NOD) = WRKSP(IDUNK) 110 CONTINUE 120 CONTINUE 130 CONTINUE C C /* GENERATE OUTPUT * / C DO 140 I = 1,NOUT C CALL OUTPUT(OUTFNC(I),OUTTYP(I),TABX,NTABX,TABY,NTABY,WRKSP, . UNKVCT) C 140 CONTINUE C RETURN C 9001 FORMAT (39H DISCRETIZATION BY INTERIOR COLLOCATION) 9011 FORMAT (38H DISCRETIZATION BY HERMITE COLLOCATION) 9021 FORMAT (23H0MATRIX COEF AND IDCOEF/) 9031 FORMAT (1X,I5,10F10.3/ (6X,10F10.3)) 9041 FORMAT (6X,10I10) END C////////////////////////////////////////////////////////////////// C****************** ALGORITHMS INTCOL & HERMCOL ******************* C////////////////////////////////////////////////////////////////// C>>>>>>>>>>>>>>> LOGICAL FILE 4: GENERAL UTILITIES <<<<<<<<<<<<< C///////////////////////////////////////////////////////////////// C FUNCTION UNQU(X,Y) C X = 0.0 Y = 0.0 UNQU = 1.0 RETURN END SUBROUTINE BASE(B1,B2,B3,B4,Z,H) C C C PURPOSE C C DEFINE THE CUBIC HERMITE POLYNOMIALS C C PARAMETERS C C INPUT C Z - DIFFERENCE BETWEEN THE POINT WHERE THE POLYNOMIALS C ARE TO BE EVALUATED AND THE LOWER ENDPOINT OF THE C INTERVAL THE POINT IS CONTAINED IN C H - LENGTH OF INTERVAL C C OUTPUT C B1,B2,B3,B4 - VALUES OF THE FOUR POLYNOMIALS C C S = Z/H T = 1. - S B1 = S*S* (2.*S-3.) + 1. B2 = 1. - B1 B3 = S*T* (H-Z) B4 = -S*T*Z RETURN END FUNCTION COLAPR(X,Y,UNKVCT,DERVSL) C C PURPOSE C C EVALUATE THE APPROXIMATE SOLUTION C AND ITS DERIVATIVES AT (X,Y) C C PARAMETERS C C X,Y - POINT AT WHICH TO EVALUATE SOLUTION C UNKVCT - SOLUTION VECTOR. CONTAINS APPROXIMATION OF C U,UX,UY AND UXY AT THE NODES C DERVSL - RETURN DERIVATIVES AND SOLUTION C DERVSL(1)=UXX C DERVSL(2)=UXY C DERVSL(3)=UYY C DERVSL(4)=UX C DERVSL(5)=UY C DERVSL(6)=U C COLAPR =U C C COMMON /GRIDXZ/GRIDX(1) COMMON /GRIDYZ/GRIDY(1) COMMON /PROBI/NGRIDX,NGRIDY COMMON /COLNUM/NODELM(4,1) C REAL UNKVCT(4,1),DERVSL(6),Q(16) INTEGER UR,UL C C ZERO OUT VECTOR OF COEFFICIENT VALUES C DO 10 I = 1,16 Q(I) = 0. 10 CONTINUE C C C FIND THE ELEMENT TO WHICH (X,Y) BELONGS C I = 1 20 CONTINUE I = I + 1 IF (GRIDX(I).LT.X) GO TO 20 I = I - 1 C J = 1 30 CONTINUE J = J + 1 IF (GRIDY(J).LT.Y) GO TO 30 J = J - 1 C C COMPUTE SIZES OF ELEMENT (I,J) C DX = GRIDX(I+1) - GRIDX(I) DY = GRIDY(J+1) - GRIDY(J) C C COMPUTE ELEMENT INDICES C LL = J + (I-1)*NGRIDY LR = LL + NGRIDY UR = LR + 1 UL = LL + 1 C C LOCATE DEGREES OF FREEDOM ASSOCIATED WITH ELEMENT (I,J) C DO 40 NDF = 1,4 IF (LL.NE.0) Q(NDF) = UNKVCT(NDF,LL) NDF1 = NDF + 4 IF (LR.NE.0) Q(NDF1) = UNKVCT(NDF,LR) NDF2 = NDF1 + 4 IF (UR.NE.0) Q(NDF2) = UNKVCT(NDF,UR) NDF3 = NDF2 + 4 IF (UL.NE.0) Q(NDF3) = UNKVCT(NDF,UL) 40 CONTINUE C C COMPUTE BASIS FUNCTIONS AT (X,Y) C XP = X - GRIDX(I) YP = Y - GRIDY(J) C CALL BASE(BX1,BX2,BX3,BX4,XP,DX) CALL BASE(BY1,BY2,BY3,BY4,YP,DY) CALL DBASE(DBX1,DBX2,DBX3,DBX4,XP,DX) CALL DBASE(DBY1,DBY2,DBY3,DBY4,YP,DY) CALL DDBASE(DDBX1,DDBX2,DDBX3,DDBX4,XP,DX) CALL DDBASE(DDBY1,DDBY2,DDBY3,DDBY4,YP,DY) C C TEMP1 = Q(1)*BX1*BY1 + Q(2)*BX1*BY3 + Q(3)*BX3*BY1 + Q(4)*BX3*BY3 TEMP2 = Q(5)*BX2*BY1 + Q(6)*BX2*BY3 + Q(7)*BX4*BY1 + Q(8)*BX4*BY3 TEMP3 = Q(9)*BX2*BY2 + Q(10)*BX2*BY4 + Q(11)*BX4*BY2 + . Q(12)*BX4*BY4 TEMP4 = Q(13)*BX1*BY2 + Q(14)*BX1*BY4 + Q(15)*BX3*BY2 + . Q(16)*BX3*BY4 C COLAPR = TEMP1 + TEMP2 + TEMP3 + TEMP4 DERVSL(6) = COLAPR C C CALCULATE THE DY - DERIVATIVE OF THE COLLOCATION APPROXIMATION C C C TEMP1 = Q(1)*BX1*DBY1 + Q(2)*BX1*DBY3 + Q(3)*BX3*DBY1 + . Q(4)*BX3*DBY3 TEMP2 = Q(5)*BX2*DBY1 + Q(6)*BX2*DBY3 + Q(7)*BX4*DBY1 + . Q(8)*BX4*DBY3 TEMP3 = Q(9)*BX2*DBY2 + Q(10)*BX2*DBY4 + Q(11)*BX4*DBY2 + . Q(12)*BX4*DBY4 TEMP4 = Q(13)*BX1*DBY2 + Q(14)*BX1*DBY4 + Q(15)*BX3*DBY2 + . Q(16)*BX3*DBY4 C DERVSL(5) = TEMP1 + TEMP2 + TEMP3 + TEMP4 C C CALCULTE THE DX - DERIVATIVE OF THE APPROXIMATION C C C TEMP1 = Q(1)*DBX1*BY1 + Q(2)*DBX1*BY3 + Q(3)*DBX3*BY1 + . Q(4)*DBX3*BY3 TEMP2 = Q(5)*DBX2*BY1 + Q(6)*DBX2*BY3 + Q(7)*DBX4*BY1 + . Q(8)*DBX4*BY3 TEMP3 = Q(9)*DBX2*BY2 + Q(10)*DBX2*BY4 + Q(11)*DBX4*BY2 + . Q(12)*DBX4*BY4 TEMP4 = Q(13)*DBX1*BY2 + Q(14)*DBX1*BY4 + Q(15)*DBX3*BY2 + . Q(16)*DBX3*BY4 C DERVSL(4) = TEMP1 + TEMP2 + TEMP3 + TEMP4 C C CALCULATE THE DYY - DERIVATIVE OF THE APPROXIMATION C C C TEMP1 = Q(1)*BX1*DDBY1 + Q(2)*BX1*DDBY3 + Q(3)*BX3*DDBY1 + . Q(4)*BX3*DDBY3 TEMP2 = Q(5)*BX2*DDBY1 + Q(6)*BX2*DDBY3 + Q(7)*BX4*DDBY1 + . Q(8)*BX4*DDBY3 TEMP3 = Q(9)*BX2*DDBY2 + Q(10)*BX2*DDBY4 + Q(11)*BX4*DDBY2 + . Q(12)*BX4*DDBY4 TEMP4 = Q(13)*BX1*DDBY2 + Q(14)*BX1*DDBY4 + Q(15)*BX3*DDBY2 + . Q(16)*BX3*DDBY4 C DERVSL(3) = TEMP1 + TEMP2 + TEMP3 + TEMP4 C C CALCULATE THE DXY - DERIVATIVE OF THE APPROXIMATION C C C TEMP1 = Q(1)*DBX1*DBY1 + Q(2)*DBX1*DBY3 + Q(3)*DBX3*DBY1 + . Q(4)*DBX3*DBY3 TEMP2 = Q(5)*DBX2*DBY1 + Q(6)*DBX2*DBY3 + Q(7)*DBX4*DBY1 + . Q(8)*DBX4*DBY3 TEMP3 = Q(9)*DBX2*DBY2 + Q(10)*DBX2*DBY4 + Q(11)*DBX4*DBY2 + . Q(12)*DBX4*DBY4 TEMP4 = Q(13)*DBX1*DBY2 + Q(14)*DBX1*DBY4 + Q(15)*DBX3*DBY2 + . Q(16)*DBX3*DBY4 C DERVSL(2) = TEMP1 + TEMP2 + TEMP3 + TEMP4 C C CALCULATE THE DXX - DERIVATIVE OF THE APPROXIMATION C C TEMP1 = Q(1)*DDBX1*BY1 + Q(2)*DDBX1*BY3 + Q(3)*DDBX3*BY1 + . Q(4)*DDBX3*BY3 TEMP2 = Q(5)*DDBX2*BY1 + Q(6)*DDBX2*BY3 + Q(7)*DDBX4*BY1 + . Q(8)*DDBX4*BY3 TEMP3 = Q(9)*DDBX2*BY2 + Q(10)*DDBX2*BY4 + Q(11)*DDBX4*BY2 + . Q(12)*DDBX4*BY4 TEMP4 = Q(13)*DDBX1*BY2 + Q(14)*DDBX1*BY4 + Q(15)*DDBX3*BY2 + . Q(16)*DDBX3*BY4 C DERVSL(1) = TEMP1 + TEMP2 + TEMP3 + TEMP4 RETURN END SUBROUTINE DBASE(DB1,DB2,DB3,DB4,Z,H) C C C PURPOSE C C DEFINE THE FIRST DERIVATIVE OF THE C CUBIC HERMITE POLYNOMIALS C C PARAMETERS C C INPUT C Z - DIFFERENCE BETWEEN THE POINT WHERE THE POLYNOMIALS C ARE TO BE EVALUATED AND THE LOWER ENDPOINT OF THE C INTERVAL THE POINT IS CONTAINED IN C H - LENGTH OF INTERVAL C C OUTPUT C DB1,DB2,DB3,DB4 - VALUES OF THE FIRST DERIVATIVES OF C THE FOUR POLYNOMIALS C C RH = 1./H S = Z*RH T = 1. - S DB1 = 6.* (S-1.)*S*RH DB2 = -DB1 DB3 = T* (T-2.*S) DB4 = S* (S-2.*T) RETURN END SUBROUTINE DDBASE(DDB1,DDB2,DDB3,DDB4,Z,H) C C C PURPOSE C C DEFINE THE SECOND DERIVATIVES OF THE C CUBIC HERMITE POLYNOMIALS C C PARAMETERS C C INPUT C Z - DIFFERENCE BETWEEN THE POINT WHERE THE POLYNOMIALS C ARE TO BE EVALUATED AND THE LOWER ENDPOINT OF THE C INTERVAL THE POINT IS CONTAINED IN C H - LENGTH OF INTERVAL C C OUTPUT C DDB1,DDB2,DDB3,DDB4 - VALUES OF THE SECOND DERIVATIVES C OF THE FOUR POLYNOMIALS C C RH = 1./H S = Z*RH T = 1. - S DDB1 = 6.* (2.*S-1.)*RH*RH DDB2 = -DDB1 DDB3 = 2.* (S-2.*T)*RH DDB4 = 6.* (S-T)*RH - DDB3 RETURN END C C////////////////////////////////////////////////////////////////// C******************** ALGORITHMS INTCOL & HERMCOL ***************** C///////////////////////////////////////////////////////////////// C>>>>>>>>>>>>>>>> LOGICAL FILE 5 : ALGORITHM HERMCOL <<<<<<<<<<<<< C///////////////////////////////////////////////////////////////// SUBROUTINE HRMCOL(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ, . UNKVCT,NUMUNK,WRKSP,HOMOBC,BCTYPE,P1,P2,LEVEL) C C -- PURPOSE -- C C APPLY COLLOCATION METHOD BASED ON HERMITE BICUBICS C TO DISCRETIZE THE PDE PROBLEM (1),(2) IN CASE B.C C ARE OF MIXED TYPE. C C -- DESCRIPTION -- C C HRMCOL GENERATES ,FOR EACH ELEMENT, EXPLICITLY THE C BOUNDARY AND INTERIOR COLLOCATION EQUATIONS. THE C PARAMETERS P1,P2, DETERMINE THE LOCATION OF TWO C BOUNDARY POINTS AT EACH ELEMENT BOUNDARY SIDE. IN CASE C P1=P2=0 (DEFAULT CASE) THE BOUNDARY COLLOCATION POINTS C ARE SELECTED TO BE THE TWO GAUSS POINTS. C C -- REFERENCES -- C C E.N.HOUSTIS,W.F.MITCHELL AND J.R.RICE,COLLOCATION C SOFTWARE FOR SECOND ORDER ELLIPTIC PARTIAL DIFFERENTIAL C TRANSACTIONS OF MATH. SOFT. TO APPEAR. C C -- DECLARATIONS -- C COMMON /INTEGS/NUMBEQ,NUMCOE,LEVELZ,MOUTPT,NROW REAL GRIDX(*),GRIDY(*),COEF(MXNEQ,*),UNKVCT(4,*),WRKSP(*),P1,P2 INTEGER IDCOEF(MXNEQ,*),NUMUNK(4,*),LEVEL CHARACTER *4 BCTYPE LOGICAL HOMOBC C C -- INPUT ARGUMENTS -- C C GRIDX,NGRIDX,GRIDY,NGRIDY,BCTYPE,HOMOBC C (DEFINED IN S/R RCTCOL) C P1,P2 REAL VARIABLES C SPECIFY THE BOUNDARY COLLOCATION C POINTS C C -- OUTPUT ARGUMENTS -- C C COEF,IDCOEF,UNKVCT,NUMUNK,WRKSP C (DEFINED IN S/R RCTCOL) C C C -- INITIALIZATIONS -- C NNODES = NGRIDX*NGRIDY DO 20 NOD = 1,NNODES DO 10 NDF = 1,4 UNKVCT(NDF,NOD) = 0. NUMUNK(NDF,NOD) = 4* (NOD-1) + NDF 10 CONTINUE 20 CONTINUE C PRINT*,'UNKVCT',UNKVCT,NUMUNK C C /* ALGORITHM HERMITE COLLOCATION */ C C PROCESS P1: GENERATE DETAIL GRID INFORMATION C C *NOTE* IN THIS CASE THE USER SUPPLIES THIS INFORMATION C C ENDP1 C PROCESS P2: FIND INTERIOR X - COLLOCATION POINTS, EVALUATE C BASIS FUNCTIONS AND DERIVATIVES AT THESE POINTS C NRDX = 26* (NGRIDX-1) C CALL INTPTS(GRIDX,NGRIDX,WRKSP(1),NRDX,LEVEL) C ENDP2 C C PROCESS P3: FIND INTERIOR Y - COLLOCATION POINTS, EVALUATE C BASIS FUNCTIONS AND DERIVATIVES AT THESE POINTS C NRDY = 26* (NGRIDY-1) LRS = NRDX + 1 LWK1 = LRS + NRDY LWK2 = LWK1 + NRDX C CALL INTPTS(GRIDY,NGRIDY,WRKSP(LRS),NRDY,LEVEL) C ENDP3 C C PROCESS P4: LOOP OVER ELEMENTS , K=1 TO NX*NY , OF DOMAIN C AND COMPUTE THE COEFFICIENTS FROM COLLOCATION C OF PDE AT 4 INTERIOR POINTS OF EACH ELEMENT C AND THE COEFFICIENTS FROM COLLOCATION OF B.C C AND THE ASSOCIATED BOUNDARY COLLOCATION POINTS. C NUMBEQ = 4*NGRIDX*NGRIDY MXNCOE = 17 NUMCOE = 17 C C CHECK FOR DEFAULT CASE IN SPECIFYING THE BOUNDARY C COLLOCATION POINTS : IF ICASE = 1 THEN BOUNDARY C COLLOCATION POINTS ARE SPECIFIED ELSE USE DEFAILT C OPTION. IF (P1.NE.0. .OR. P2.NE.0.) GO TO 30 C C USE DEFAULT OPTION C CALL COLEQT(COEF,MXNCOE,MXNEQ,GRIDX,NGRIDX,GRIDY,NGRIDY,WRKSP(1), . WRKSP(1),NRDX,WRKSP(LRS),WRKSP(LRS),NRDY) C GO TO 40 30 CONTINUE C USE SPECIFIED BOUNDARY COLLOCATION POINTS IN X-VARIABLE CALL BNDPTS(GRIDX,NGRIDX,NGRIDX,WRKSP(LWK1),NRDX,P1,P2,LEVEL) C USE SPECIFIED BOUNDARY COLLOCATION POINTS IN Y-VARIABLE C CALL BNDPTS(GRIDY,NGRIDY,NGRIDY,WRKSP(LWK2),NRDY,P1,P2,LEVEL) C PROCESS P5: C GENERATE COLLOCATION EQUATIONS C CALL COLEQT(COEF,MXNCOE,MXNEQ,GRIDX,NGRIDX,GRIDY,NGRIDY,WRKSP(1), . WRKSP(LWK1),NRDX,WRKSP(LRS),WRKSP(LWK2),NRDY) C 40 CONTINUE C ENDP5 C C PROCESS P6 : ORDER THE COLLOCATION EQUATIONS C CALL ORDEQT(IDCOEF,MXNCOE,MXNEQ,LEVEL) C C ENDP6 RETURN END SUBROUTINE BNDEQN(ISIDE,T,BS,COEF,MXNEQ,MXNCOE,HXX,HYY,RHXHY) C C***** DISCRETIZES THE BOUNDARY CONDITIONS BY FORCING THE C COLLOCATION APPROXIMATION 'UC' TO SATISFY EXACTLY C THE BOUNDARY CONDITIONS AT CERTAIN BOUNDARY POINTS ***** C C***** BOUNDARY CONDITIONS ASSUMED: A(X,Y)U + B(X,Y)UX ON SIDES C 1,3 AND A(X,Y)U + B(X,Y)UY ON SIDES 2,4 ****** C C***** COLLOCATION DISCRETIZATION OF B.C: C C ---- BOUNDARY SIDE 1---- C C COLLOCATION APPROX.: UC = Q(5)BY1+Q(6)BY3+Q(9)BY2+Q(10)BY4 C DUC/DX = Q(7)BY1+Q(8)BY3+Q(11)BY2+Q(12)BY4 C C DISCRETIZED EQUATION : A(BX,.)UC +B(BX,.)DUC/DX = G(BX,.) C C ---- BOUNDARY SIDE 3---- C C COLLOCATION APPROX.: UC = Q(1)BY1+Q(2)BY3+Q(13)BY2+Q(14)BY4 C DUC/DX = Q(3)BY1+Q(4)BY3+Q(15)BY2+Q(16)BY4 C C DISCRETIZED EQUATION : A(AX,.)UC +B(AX,.)DUC/DX = G(AX,.) C C ---- BOUNDARY SIDE 2---- C C COLLOCATION APPROX.: UC = Q(1)BX1+Q(2)BX1+Q(5)BX2+Q(6)BX2 C DUC/DY = Q(3)BX3+Q(4)BX3+Q(7)BX4+Q(8)BX4 C C DISCRETIZED EQUATION : A(AY,.)UC +B(AY,.)DUC/DY = G(AY,.) C C ---- BOUNDARY SIDE 4---- C C COLLOCATION APPROX.: UC = Q(13)BX1+Q(14)BX1+Q( 9)BX2+Q(10)BX2 C DUC/DY = Q(15)BX3+Q(16)BX3+Q(11)BX4+Q(12)BY4 C C DISCRETIZED EQUATION : A(BY,.)UC +B(BY,.)DUC/DY = G(BY,.) C C ****** C C ----- INPUT ARGUMENTS :: ISIDE : BOUNDARY SIDE , TÕ1:2å : BOUNDARY C COLLOCATION POINTS ASSOCIATED ON THIS SIDE C BSÕ1:24å : VALUES OF BASES FUNCTIONS AND C THEIR DERIVATIVES AT THESE POINTS, C ----- OUTPUT ARGUMENTS :: C COEFÕNROW:NROW+1,1:16å : COEFFICIENTS OF C BOUNDARY COLLOCATION EQUATIONS, C COEFÕNROW:NROW+1,17å : RIGHT SIDE OF THE C THE EQUATIONS, MXNEQ: NUMBER OF COLLOCATION C EQUATIONS C C ----- COMMON VARIABLES ---- C C COMMON / SCALE / HXX,HYY,RHXHY C ----- DECLARATIONS ----- C REAL COEF(MXNEQ,*),T(2),BS(24),BVALUS(4) COMMON /INTEGS/NUMBEQ,NUMCOE,LEVEL,MOUTPT,NROW COMMON /PROBR/AX,BX,AY,BY C C ----- LOCAL VARIABLES ----- BS1 = BS(1) BS2 = BS(2) BS3 = BS(3) BS4 = BS(4) C C ----- PROCEDURE FOR GENERATING THE BOUNDARY COLLOCATION EQS -- C C SELECT BOUNDARY SIDE C GO TO (10,30,50,70),ISIDE C GENERATE BOUNDARY COLLOCATION EQUATIONS ON SIDE 1 C 10 CONTINUE DO 20 IC = 1,2 C COMPUTE FACTOR TO SCALE THE EQUATIONS RS = BCOND(1,BX,T(IC),BVALUS) SCAL = RHXHY IF (BVALUS(2).NE.0.) SCAL = SCAL*HXX COEF(NROW,MXNCOE) = RS*SCAL C BV1 = BVALUS(1)*SCAL BV2 = BVALUS(2)*SCAL C COEF(NROW,5) = BV1*BS1 COEF(NROW,6) = BV1*BS3 COEF(NROW,7) = BV2*BS1 COEF(NROW,8) = BV2*BS3 C COEF(NROW,9) = BV1*BS2 COEF(NROW,10) = BV1*BS4 COEF(NROW,11) = BV2*BS2 COEF(NROW,12) = BV2*BS4 C C REDEFINES TEMPORARIES FOR SECOND PASS THROUGH THE LOOP C BS1 = BS(13) BS2 = BS(14) BS3 = BS(15) BS4 = BS(16) C INCREASE EQUATION COUNTER NROW = NROW + 1 20 CONTINUE GO TO 90 C 30 CONTINUE C C GENERATE BOUNDARY COLLOCATION EQUATIONS ON SIDE 2 C DO 40 IC = 1,2 C COMPUTE FACTOR TO SCALE THE EQUATIONS RS = BCOND(2,T(IC),AY,BVALUS) SCAL = RHXHY IF (BVALUS(3).NE.0.) SCAL = SCAL*HYY COEF(NROW,MXNCOE) = RS*SCAL C BV1 = BVALUS(1)*SCAL BV3 = BVALUS(3)*SCAL COEF(NROW,1) = BV1*BS1 COEF(NROW,2) = BV3*BS1 COEF(NROW,3) = BV1*BS3 COEF(NROW,4) = BV3*BS3 C COEF(NROW,5) = BV1*BS2 COEF(NROW,6) = BV3*BS2 COEF(NROW,7) = BV1*BS4 COEF(NROW,8) = BV3*BS4 C C REDEFINES TEMPORARIES FOR SECOND PASS THROUGH THE LOOP C BS1 = BS(13) BS3 = BS(15) BS2 = BS(14) BS4 = BS(16) C INCREASE EQUATION COUNTER NROW = NROW + 1 40 CONTINUE GO TO 90 C 50 CONTINUE C C GENERATE BOUNDARY COLLOCATION EQUATIONS ON SIDE 3 C DO 60 IC = 1,2 C COMPUTE FACTOR TO SCALE THE EQUATIONS RS = BCOND(3,AX,T(IC),BVALUS) SCAL = RHXHY IF (BVALUS(2).NE.0.) SCAL = SCAL*HXX COEF(NROW,MXNCOE) = RS*SCAL C BV1 = BVALUS(1)*SCAL BV2 = BVALUS(2)*SCAL C C COEF(NROW,1) = BV1*BS1 COEF(NROW,2) = BV1*BS3 COEF(NROW,3) = BV2*BS1 COEF(NROW,4) = BV2*BS3 C COEF(NROW,13) = BV1*BS2 COEF(NROW,14) = BV1*BS4 COEF(NROW,15) = BV2*BS2 COEF(NROW,16) = BV2*BS4 C C REDEFINES TEMPORARIES FOR SECOND PASS THROUGH THE LOOP C BS1 = BS(13) BS3 = BS(15) BS2 = BS(14) BS4 = BS(16) C INCREASE EQUATION COUNTER NROW = NROW + 1 60 CONTINUE GO TO 90 C 70 CONTINUE C C GENERATE BOUNDARY COLLOCATION EQUATIONS ON SIDE 4 C DO 80 IC = 1,2 C COMPUTE FACTOR TO SCALE THE EQUATIONS RS = BCOND(4,T(IC),BY,BVALUS) SCAL = RHXHY IF (BVALUS(3).NE.0.) SCAL = SCAL*HYY COEF(NROW,MXNCOE) = RS*SCAL C BV1 = BVALUS(1)*SCAL BV3 = BVALUS(3)*SCAL COEF(NROW,13) = BV1*BS1 COEF(NROW,14) = BV3*BS1 COEF(NROW,15) = BV1*BS3 COEF(NROW,16) = BV3*BS3 C COEF(NROW,9) = BV1*BS2 COEF(NROW,10) = BV3*BS2 COEF(NROW,11) = BV1*BS4 COEF(NROW,12) = BV3*BS4 C C REDEFINES TEMPORARIES FOR SECOND PASS THROUGH THE LOOP C BS1 = BS(13) BS3 = BS(15) BS2 = BS(14) BS4 = BS(16) C INCREASE EQUATION COUNTER NROW = NROW + 1 80 CONTINUE 90 RETURN END SUBROUTINE BNDPTS(X,NDKNTS,KNTS,R,NRD,P1,P2,LEVEL) C REAL X(NDKNTS),R(NRD) C C SET 1-D COLLOCATION POINTS FOR BOUNDARY, BASED ON PARAMETERS P1, P2 C AND COMPUTE BASIS FUNCTIONS AND THEIR DERIVATIVES C C3 = 1./3. NINT = KNTS - 1 IBL1 = 1 C C GENERATE TWO POINTS IN EACH INTERVAL C DO 10 I = 1,NINT IP1 = I + 1 XLEFT = X(I) H = X(IP1) - XLEFT C TDUM = P1*H R(IBL1) = XLEFT + TDUM C IBL3 = IBL1 + 2 CALL HBASVD(R(IBL3),TDUM,H) C IBL2 = IBL1 + 1 TDUM = P2*H R(IBL2) = XLEFT + TDUM C IBL4 = IBL1 + 14 CALL HBASVD(R(IBL4),TDUM,H) C IBL1 = IBL1 + 26 10 CONTINUE C C IF P1=0, THE FIRST INTERVAL MUST BE ADJUSTED C IF (P1.NE.0.) GO TO 20 C H = X(2) - X(1) H3 = H*C3 C R(1) = X(1) + H3 CALL HBASVD(R(3),H3,H) C R(2) = R(1) + H3 TDUM = 2.*H3 CALL HBASVD(R(15),TDUM,H) C 20 CONTINUE C C IF P2=1, THE LAST INTERVAL MUST BE ADJUSTED C IF (P2.NE.1.) GO TO 30 C H3 = H*C3 IBL1 = IBL1 - 26 R(IBL1) = X(NINT) + H3 C IBL3 = IBL1 + 2 CALL HBASVD(R(IBL3),H3,H) C TDUM = 2.*H3 IBL2 = IBL1 + 1 R(IBL2) = R(IBL1) + H3 C IBL4 = IBL1 + 14 CALL HBASVD(R(IBL4),TDUM,H) C 30 CONTINUE IF (LEVEL.LT.4) GO TO 40 WRITE (6,9001) WRITE (6,9011) (R(IR),IR=1,NRD) 40 CONTINUE RETURN C 9001 FORMAT (39H PRINTS 1-D BASES , 1-ST AND 2-ND DERIV) 9011 FORMAT (1X,8F10.4) END SUBROUTINE COLEQT(COEF,MXNCOE,MXNEQ,GRIDX,NGRDXD,GRIDY,NGRDYD,RS1, . WK1,NDRSX,RS2,WK2,NDRSY) C C ***** FORMS THE DISCRETE COLLOCATION EQUATIONS IN THE C CASE OF NON-HOMOGENEOUS BOUNDARY CONDITIONS **** C C DECLARATIONS C REAL COEF(MXNEQ,MXNCOE),BVALUS(4),GRIDX(NGRDXD),GRIDY(NGRDYD), . RS1(NDRSX),RS2(NDRSY),WK1(NDRSX),WK2(NDRSY) REAL AX,BX,AY,BY C COMMON /PROBR/AX,BX,AY,BY COMMON /PROBI/NGRIDX,NGRIDY COMMON /INTEGS/NUMBEQ,NUMCOE,LEVEL,MOUTPT,NROW COMMON /BNDRY/IPIECE,NBOUND,NBNDPT C INITIALIZATION DO 20 IR = 1,MXNEQ DO 10 IC = 1,MXNCOE COEF(IR,IC) = 0. 10 CONTINUE 20 CONTINUE C NX = NGRIDX - 1 NY = NGRIDY - 1 NYM1 = NY - 1 NXM1 = NX - 1 C C C **** CASE I **** C NROW = 1 HXX = GRIDX(2) - GRIDX(1) HYY = GRIDY(2) - GRIDY(1) RHXHY = (1.+HYY)* (1.+1./HXX)/HXX + (1.+HXX)* (1.+1./HYY)/HYY SCAL = RHXHY RS = BCOND(2,AX,AY,BVALUS) IF (BVALUS(2).NE.0.) SCAL = SCAL*HXX IF (BVALUS(3).NE.0.) SCAL = SCAL*HYY C COEF(1,MXNCOE) = RS*SCAL COEF(1,1) = BVALUS(1)*SCAL COEF(1,2) = BVALUS(3)*SCAL COEF(1,3) = BVALUS(2)*SCAL C NROW = NROW + 1 CALL BNDEQN(2,WK1(1),WK1(3),COEF,MXNEQ,MXNCOE,HXX,HYY,RHXHY) CALL BNDEQN(3,WK2(1),WK2(3),COEF,MXNEQ,MXNCOE,HXX,HYY,RHXHY) C CALL EQBLCK(NROW,RS1(1),RS1(3),RS1(15),RS2(1),RS2(3),RS2(15),COEF, . MXNEQ,MXNCOE) C C CASE II C IF (NGRIDY.EQ.3) GO TO 40 DO 30 JY = 2,NYM1 HYY = GRIDY(JY+1) - GRIDY(JY) RHXHY = (1.+HYY)* (1.+1./HXX)/HXX + (1.+HXX)* (1.+1./HYY)/HYY IBLY1 = 26*JY - 25 IBLY2 = IBLY1 + 2 CALL BNDEQN(3,WK2(IBLY1),WK2(IBLY2),COEF,MXNEQ,MXNCOE,HXX,HYY, . RHXHY) IBLY3 = IBLY2 + 12 CALL EQBLCK(NROW,RS1(1),RS1(3),RS1(15),RS2(IBLY1),RS2(IBLY2), . RS2(IBLY3),COEF,MXNEQ,MXNCOE) 30 CONTINUE 40 CONTINUE C C CASE III C C HYY = GRIDY(NGRIDY) - GRIDY(NY) RHXHY = (1.+HYY)* (1.+1./HXX)/HXX + (1.+HXX)* (1.+1./HYY)/HYY IBLY1 = 26*NY - 25 IBLY2 = IBLY1 + 2 CALL BNDEQN(3,WK2(IBLY1),WK2(IBLY2),COEF,MXNEQ,MXNCOE,HXX,HYY, . RHXHY) C RS = BCOND(3,AX,BY,BVALUS) SCAL = RHXHY IF (BVALUS(2).NE.0.) SCAL = SCAL*HXX IF (BVALUS(3).NE.0.) SCAL = SCAL*HYY COEF(NROW,MXNCOE) = RS*SCAL COEF(NROW,13) = BVALUS(1)*SCAL COEF(NROW,14) = BVALUS(3)*SCAL COEF(NROW,15) = BVALUS(2)*SCAL NROW = NROW + 1 C CALL BNDEQN(4,WK1(1),WK1(3),COEF,MXNEQ,MXNCOE,HXX,HYY,RHXHY) C IBLY3 = IBLY2 + 12 CALL EQBLCK(NROW,RS1(1),RS1(3),RS1(15),RS2(IBLY1),RS2(IBLY2), . RS2(IBLY3),COEF,MXNEQ,MXNCOE) C IF (NGRIDX.EQ.3) GO TO 80 DO 70 IX = 2,NXM1 C C **** CASE IV **** C IBLX1 = 26*IX - 25 IBLX2 = IBLX1 + 2 HXX = GRIDX(IX+1) - GRIDX(IX) HYY = GRIDY(2) - GRIDY(1) RHXHY = (1.+HYY)* (1.+1./HXX)/HXX + (1.+HXX)* (1.+1./HYY)/HYY CALL BNDEQN(2,WK1(IBLX1),WK1(IBLX2),COEF,MXNEQ,MXNCOE,HXX,HYY, . RHXHY) C IBLX3 = IBLX2 + 12 CALL EQBLCK(NROW,RS1(IBLX1),RS1(IBLX2),RS1(IBLX3),RS2(1), . RS2(3),RS2(15),COEF,MXNEQ,MXNCOE) C C ***** CASE V **** C IF (NGRIDY.EQ.3) GO TO 60 DO 50 JY = 2,NYM1 IBLY1 = 26*JY - 25 IBLY2 = IBLY1 + 2 IBLY3 = IBLY2 + 12 CALL EQBLCK(NROW,RS1(IBLX1),RS1(IBLX2),RS1(IBLX3), . RS2(IBLY1),RS2(IBLY2),RS2(IBLY3),COEF,MXNEQ, . MXNCOE) 50 CONTINUE 60 CONTINUE C C **** CASE VI ***** C IBLY1 = 26*NY - 25 IBLY2 = IBLY1 + 2 HYY = GRIDY(NGRIDY) - GRIDY(NY) RHXHY = (1.+HYY)* (1.+1./HXX)/HXX + (1.+HXX)* (1.+1./HYY)/HYY CALL BNDEQN(4,WK1(IBLX1),WK1(IBLX2),COEF,MXNEQ,MXNCOE,HXX,HYY, . RHXHY) C IBLY3 = IBLY2 + 12 CALL EQBLCK(NROW,RS1(IBLX1),RS1(IBLX2),RS1(IBLX3),RS2(IBLY1), . RS2(IBLY2),RS2(IBLY3),COEF,MXNEQ,MXNCOE) 70 CONTINUE 80 CONTINUE C C **** CASE VII **** C C IBLX1 = 26*NX - 25 IBLX2 = IBLX1 + 2 HXX = GRIDX(NGRIDX) - GRIDX(NX) HYY = GRIDY(2) - GRIDY(1) RHXHY = (1.+HYY)* (1.+1./HXX)/HXX + (1.+HXX)* (1.+1./HYY)/HYY CALL BNDEQN(2,WK1(IBLX1),WK1(IBLX2),COEF,MXNEQ,MXNCOE,HXX,HYY, . RHXHY) C RS = BCOND(2,BX,AY,BVALUS) SCAL = RHXHY IF (BVALUS(2).NE.0.) SCAL = SCAL*HXX IF (BVALUS(3).NE.0.) SCAL = SCAL*HYY COEF(NROW,MXNCOE) = RS*SCAL COEF(NROW,5) = BVALUS(1)*SCAL COEF(NROW,6) = BVALUS(3)*SCAL COEF(NROW,7) = BVALUS(2)*SCAL NROW = NROW + 1 C CALL BNDEQN(1,WK2(1),WK2(3),COEF,MXNEQ,MXNCOE,HXX,HYY,RHXHY) C IBLX3 = IBLX2 + 12 CALL EQBLCK(NROW,RS1(IBLX1),RS1(IBLX2),RS1(IBLX3),RS2(1),RS2(3), . RS2(15),COEF,MXNEQ,MXNCOE) C C CASE VIII C IF (NGRIDY.EQ.3) GO TO 100 DO 90 JY = 2,NYM1 IBLY1 = 26*JY - 25 IBLY2 = IBLY1 + 2 HYY = GRIDY(JY+1) - GRIDY(JY) RHXHY = (1.+HYY)* (1.+1./HXX)/HXX + (1.+HXX)* (1.+1./HYY)/HYY CALL BNDEQN(1,WK2(IBLY1),WK2(IBLY2),COEF,MXNEQ,MXNCOE,HXX,HYY, . RHXHY) IBLY3 = IBLY2 + 12 CALL EQBLCK(NROW,RS1(IBLX1),RS1(IBLX2),RS1(IBLX3),RS2(IBLY1), . RS2(IBLY2),RS2(IBLY3),COEF,MXNEQ,MXNCOE) 90 CONTINUE 100 CONTINUE C C CASE IX C C IBLY1 = 26*NY - 25 IBLY2 = IBLY1 + 2 HYY = GRIDY(NGRIDY) - GRIDY(NY) RHXHY = (1.+HYY)* (1.+1./HXX)/HXX + (1.+HXX)* (1.+1./HYY)/HYY CALL BNDEQN(1,WK2(IBLY1),WK2(IBLY2),COEF,MXNEQ,MXNCOE,HXX,HYY, . RHXHY) C RS = BCOND(1,BX,BY,BVALUS) SCAL = RHXHY IF (BVALUS(2).NE.0.) SCAL = SCAL*HXX IF (BVALUS(3).NE.0.) SCAL = SCAL*HYY COEF(NROW,MXNCOE) = RS*SCAL COEF(NROW,9) = BVALUS(1)*SCAL COEF(NROW,10) = BVALUS(3)*SCAL COEF(NROW,11) = BVALUS(2)*SCAL NROW = NROW + 1 C CALL BNDEQN(4,WK1(IBLX1),WK1(IBLX2),COEF,MXNEQ,MXNCOE,HXX,HYY, . RHXHY) C IBLY3 = IBLY2 + 12 CALL EQBLCK(NROW,RS1(IBLX1),RS1(IBLX2),RS1(IBLX3),RS2(IBLY1), . RS2(IBLY2),RS2(IBLY3),COEF,MXNEQ,MXNCOE) C RETURN END C C SUBROUTINE INTCOL(GRIDX,NGRIDX,GRIDY,NGRIDY,COEF,IDCOEF,MXNEQ, C . UNKVCT,NUMUNK,WRKSP,HOMOBC,BCTYPE,LEVEL) C C DUMMY FOR USE IN DIRECTORY FOR EXTERIOR COLLOCATION C C RETURN C END SUBROUTINE ORDEQT(IDCOEF,MXNCOE,MXNEQ,LEVEL) C C **** NUMBERS THE UNKNOWNS IN CASE OF NON-HOMOGENEOUS BOUNDARY COND**** C C C C DECLARATIONS INTEGER IDCOEF(MXNEQ,MXNCOE),UL,UR,COL(16) C C REAL BVALUS(4) C COMMON /PROBI/NGRIDX,NGRIDY C C INITIALIZATIONS C DO 20 IQ = 1,MXNEQ DO 10 JQ = 1,MXNCOE IDCOEF(IQ,JQ) = 0 10 CONTINUE 20 CONTINUE C C TEMPORARY VARIABLES C NX = NGRIDX - 1 NY = NGRIDY - 1 NY4 = 4*NGRIDY C C SET COLUMN INDICIES FOR FIRST ROW COL(1) = 1 COL(2) = 2 COL(3) = 3 COL(4) = 4 COL(5) = 1 + NY4 COL(6) = 2 + NY4 COL(7) = 3 + NY4 COL(8) = 4 + NY4 COL(9) = 5 + NY4 COL(10) = 6 + NY4 COL(11) = 7 + NY4 COL(12) = 8 + NY4 COL(13) = 5 COL(14) = 6 COL(15) = 7 COL(16) = 8 C C C FORM ARRAY " IDCOEF" C NROW = 0 C DO 100 IX = 1,NX C C IF NOT THE FIRST COLUMN, SHIFT COLUMNS EIGHT C IF (IX.EQ.1) GO TO 40 DO 30 I = 1,16 COL(I) = COL(I) + 8 30 CONTINUE 40 CONTINUE DO 90 JY = 1,NY C C IF NOT FIRST ROW, SHIFT COLUMNS FOUR C IF (JY.EQ.1) GO TO 60 DO 50 I = 1,16 COL(I) = COL(I) + 4 50 CONTINUE 60 CONTINUE C C NELCLP = 4 IF (IX.EQ.1 .OR. IX.EQ.NX) NELCLP = 6 IF (JY.EQ.1 .OR. JY.EQ.NY) NELCLP = 6 IF (IX.EQ.1 .AND. ((JY.EQ.NY).OR. (JY.EQ.1))) NELCLP = 9 IF (IX.EQ.NX .AND. ((JY.EQ.NY).OR. (JY.EQ.1))) NELCLP = 9 C DO 80 K = 1,NELCLP NROW = NROW + 1 DO 70 ICOL = 1,16 IDCOEF(NROW,ICOL) = COL(ICOL) 70 CONTINUE 80 CONTINUE C 90 CONTINUE 100 CONTINUE RETURN END C/////////////////////////////////////////////////////////////// C................ END OF ALGORITHMS INTCOL & HERMCOL .......... C////////////////////////////////////////////////////////////// C//////////////////////////////////////////////////////////////////// C******************ALGORITHM INTCOL & HERMCOL ********************** C//////////////////////////////////////////////////////////////////// C >>>>>>>>>>>>>>>>>> LOGICAL FILE 6 : EXAMPLE 1 <<<<<<<<<<<<<<<<<<<< C////////////////////////////////////////////////////////////////////// REAL FUNCTION BCOND(I,X,Y,BVALUS) C REAL BVALUS(4) COMMON /PARAMS/ALPHA,RL C C GO TO (10,20,30,40),I 10 BVALUS(1) = 1. BVALUS(2) = 0. BVALUS(3) = 0. BVALUS(4) = 0. BCOND = BVALUS(4) C RETURN C C 20 CONTINUE BVALUS(1) = 0. BVALUS(2) = 0. BVALUS(3) = 1. BVALUS(4) = 0. BCOND = BVALUS(4) RETURN 30 BVALUS(1) = 0. BVALUS(2) = 1. BVALUS(3) = 0. BVALUS(4) = 0. BCOND = BVALUS(4) RETURN 40 IF (X.GE..5) GO TO 50 BVALUS(1) = 0. BVALUS(2) = 0. BVALUS(3) = 1. BVALUS(4) = 0. BCOND = BVALUS(4) RETURN 50 BVALUS(1) = 1. BVALUS(2) = 0. BVALUS(3) = 0. BVALUS(4) = 0. BCOND = BVALUS(4) RETURN END SUBROUTINE BCOORD(P,X,Y,I) C C DUMMY TO SATISFY LOADING C USED WITH INTERIOR AND HERMITE ONLY IF RESIDUAL IS COMPUTED C RETURN END SUBROUTINE PDE(X,Y,CVALUS) C REAL CVALUS(6) C C 10 CONTINUE CVALUS(1) = 1. CVALUS(2) = 0. CVALUS(3) = 1./ (X*X) CVALUS(4) = 1./X CVALUS(5) = 0. CVALUS(6) = 0. RETURN END REAL FUNCTION PDERHS(X,Y) C C PDERHS = -1. RETURN END REAL FUNCTION TRUE(X,Y) C C TRUE SOLUTION IS UNKNOWN -- THIS IS A DUMMY C TRUE = 0. RETURN END C////////////////////////////////////////////////////////////////// C///////////////// END OF FILE 6 ////////////////////////////////// C/////////////////////////////////////////////////////////////////// C/////////////////////////////////////////////////////////////// C************ALGORITHMS INTCOL & HERMCOL************************* C/////////////////////////////////////////////////////////////// C>>>>>>>>>>>>>>>> LOGICAL FILE 7 : EXAMPLE 2 <<<<<<<<<<<<<<<<<<< C/////////////////////////////////////////////////////////////// REAL FUNCTION BCOND(I,X,Y,BVALUS) C REAL BVALUS(4) C C BVALUS(1) = 1. BVALUS(2) = 0. BVALUS(3) = 0. BVALUS(4) = TRUE(X,Y) BCOND = BVALUS(4) C RETURN END SUBROUTINE BCOORD(P,X,Y,I) C C DUMMY TO SATISFY LOADING C USED WITH INTERIOR AND HERMITE ONLY IF RESIDUAL IS COMPUTED C RETURN END SUBROUTINE PDE(X,Y,CVALUS) C REAL CVALUS(6) C C 10 CONTINUE CVALUS(1) = 1. CVALUS(2) = 0. CVALUS(3) = 1./ (X*X) CVALUS(4) = 2./X CVALUS(5) = 1./ (X*TAN(Y)) CVALUS(6) = 0. RETURN END REAL FUNCTION PDERHS(X,Y) C C PDERHS = (1.+1./ (X*X)+2./X+1./ (X*TAN(Y)))*EXP(X+Y) RETURN END REAL FUNCTION TRUE(X,Y) TRUE = EXP(X+Y) RETURN C////////////////////////////////////////////////////////////////// C//////////////////////// END OF FILE 7 /////////////////////////// C////////////////////////////////////////////////////////////////// C***************** END OF ALGORITHMS INTCOL & HERCOL ************* END