C FIGURE 5.9.1 IMPLICIT DOUBLE PRECISION(A-H,O-Z) C NTRI = NUMBER OF TRIANGLES C M = NUMBER OF UNKNOWNS C NSTEP = NUMBER OF TIME STEPS PARAMETER (NTRI=16, M=5, NSTEP=10) DIMENSION TRIX(3,NTRI),TRIY(3,NTRI),NODES(3,NTRI),A(M),B(M) C AMAT,BMAT DIMENSIONED (M,-L:L), BUT C HALF BANDWIDTH L UNKNOWN UNTIL EXECUTION C TIME ALLOCATABLE AMAT(:,:),BMAT(:,:) DIMENSION PHI(3),PHIX(3),PHIY(3) C ARRAYS TRIX,TRIY HOLD THE VERTICES OF THE C TRIANGLES. THE THREE VERTICES OF TRIANGLE C KTRI ARE C (TRIX(J,KTRI) , TRIY(J,KTRI)) (J=1,2,3) C AND THE CORRESPONDING NODE NUMBERS ARE C NODES(J,KTRI). C THE NUMBER OF A VERTEX ON THE BOUNDARY IS 0 DATA ((TRIX(I,L),TRIY(I,L),NODES(I,L),I=1,3),L=1,NTRI)/ & 0.0, 0.0, 0, 0.0,-1.0, 0, 0.5,-0.5, 1, & 0.0,-1.0, 0, 1.0,-1.0, 0, 0.5,-0.5, 1, & 1.0,-1.0, 0, 1.0, 0.0, 0, 0.5,-0.5, 1, & 1.0, 0.0, 0, 0.5, 0.0, 2, 0.5,-0.5, 1, & 0.5, 0.0, 2, 0.0, 0.0, 0, 0.5,-0.5, 1, & 0.0, 0.0, 0, 0.5, 0.0, 2, 0.5, 0.5, 3, & 0.5, 0.0, 2, 1.0, 0.0, 0, 0.5, 0.5, 3, & 1.0, 0.0, 0, 1.0, 1.0, 0, 0.5, 0.5, 3, & 1.0, 1.0, 0, 0.0, 1.0, 0, 0.5, 0.5, 3, & 0.0, 1.0, 0, 0.0, 0.5, 4, 0.5, 0.5, 3, & 0.0, 0.5, 4, 0.0, 0.0, 0, 0.5, 0.5, 3, & 0.0, 0.0, 0, 0.0, 0.5, 4, -0.5, 0.5, 5, & 0.0, 0.5, 4, 0.0, 1.0, 0, -0.5, 0.5, 5, & 0.0, 1.0, 0, -1.0, 1.0, 0, -0.5, 0.5, 5, & -1.0, 1.0, 0, -1.0, 0.0, 0, -0.5, 0.5, 5, & -1.0, 0.0, 0, 0.0, 0.0, 0, -0.5, 0.5, 5/ C STATEMENT FUNCTIONS DEFINING PDE C COEFFICIENTS FUNC(X,Y,T) = 1.0 FUND(X,Y,T) = 1.0 FUNA(X,Y,T) = 0.0 FUNF(X,Y,T) = -(1.0+X+Y) FUNR(X,Y,T) = (1.0-T)*(1.0+X+Y) FUNRT(X,Y,T) = -(1.0+X+Y) FUNH(X,Y) = 1.0+X+Y C CALCULATE HALF-BANDWIDTH, L L = 0 DO 20 KTRI=1,NTRI DO 15 KX=1,3 K = NODES(KX,KTRI) IF (K.GT.0) THEN C ASSIGN INITIAL VALUES A(K) = FUNH(TRIX(KX,KTRI) , TRIY(KX,KTRI)) DO 10 IX=1,3 I = NODES(IX,KTRI) IF (I.GT.0) L = MAX(L,ABS(I-K)) 10 CONTINUE ENDIF 15 CONTINUE 20 CONTINUE C NOW L IS KNOWN, SO ALLOCATE AMAT,BMAT ALLOCATE (AMAT(M,-L:L),BMAT(M,-L:L)) C TAKE NSTEP STEPS OF SIZE DT DT = 1.0D0/NSTEP DO 70 KT=1,NSTEP TKP1 = KT*DT C ZERO BAND MATRICES AMAT,BMAT AND VECTOR B. DO 30 K=1,M B(K) = 0.0 DO 25 ISUB=-L,L AMAT(K,ISUB) = 0.0 BMAT(K,ISUB) = 0.0 25 CONTINUE 30 CONTINUE C BEGIN ELEMENT-BY-ELEMENT ASSEMBLY OF C AMAT,BMAT AND B. DO 50 KTRI=1,NTRI X1 = TRIX(1,KTRI) Y1 = TRIY(1,KTRI) X2 = TRIX(2,KTRI) Y2 = TRIY(2,KTRI) X3 = TRIX(3,KTRI) Y3 = TRIY(3,KTRI) DET = (X2-X1)*(Y3-Y1) - (X3-X1)*(Y2-Y1) C PHI(J),PHIX(J),PHIY(J) REPRESENT THE C VALUE, X-DERIVATIVE AND Y-DERIVATIVE C OF PHI-J AT THE TRIANGLE MIDPOINT. PHI (1) = 1.0/3.0D0 PHIX(1) = (Y2-Y3)/DET PHIY(1) = (X3-X2)/DET PHI (2) = 1.0/3.0D0 PHIX(2) = (Y3-Y1)/DET PHIY(2) = (X1-X3)/DET PHI (3) = 1.0/3.0D0 PHIX(3) = (Y1-Y2)/DET PHIY(3) = (X2-X1)/DET C OMEGA,OMEGAT,OMEGAX,OMEGAY REPRESENT THE C VALUE, T-DERIVATIVE, X-DERIVATIVE, C AND Y-DERIVATIVE OF OMEGA AT THE TRIANGLE C MIDPOINT. OMEGA = 0.0 OMEGAT = 0.0 OMEGAX = 0.0 OMEGAY = 0.0 DO 35 J=1,3 IF (NODES(J,KTRI).LE.0) THEN R = FUNR(TRIX(J,KTRI),TRIY(J,KTRI),TKP1) RT=FUNRT(TRIX(J,KTRI),TRIY(J,KTRI),TKP1) OMEGA = OMEGA + R*PHI (J) OMEGAT = OMEGAT + RT*PHI(J) OMEGAX = OMEGAX + R*PHIX(J) OMEGAY = OMEGAY + R*PHIY(J) ENDIF 35 CONTINUE C CMID,DMID,AMID,FMID ARE THE VALUES OF THE C COEFFICIENT FUNCTIONS C(X,Y,T),D(X,Y,T), C A(X,Y,T) AND F(X,Y,T) AT THE TRIANGLE C MIDPOINT. XMID = (X1+X2+X3)/3.0 YMID = (Y1+Y2+Y3)/3.0 CMID = FUNC(XMID,YMID,TKP1) DMID = FUND(XMID,YMID,TKP1) AMID = FUNA(XMID,YMID,TKP1) FMID = FUNF(XMID,YMID,TKP1) AREA = ABS(DET) C CALCULATE CONTRIBUTIONS OF TRIANGLE KTRI C TO MATRICES A AND B AND VECTOR B, C USING THE MIDPOINT RULE TO APPROXIMATE C ALL INTEGRALS. DO 45 KX=1,3 K = NODES(KX,KTRI) IF (K.GT.0) THEN B(K) = B(K) + & (-CMID*OMEGAT*PHI(KX) & -DMID*(OMEGAX*PHIX(KX) + OMEGAY*PHIY(KX)) & -AMID*OMEGA*PHI(KX) & +FMID*PHI(KX)) *AREA DO 40 IX=1,3 I = NODES(IX,KTRI) IF (I.GT.0) THEN AMAT(K,I-K) = AMAT(K,I-K) + & (DMID*(PHIX(KX)*PHIX(IX) + PHIY(KX)*PHIY(IX)) & +AMID*PHI(KX)*PHI(IX)) *AREA BMAT(K,I-K) = BMAT(K,I-K) + & CMID*PHI(KX)*PHI(IX)*AREA ENDIF 40 CONTINUE ENDIF 45 CONTINUE 50 CONTINUE C BACKWARD DIFFERENCE METHOD REQUIRES C SOLUTION OF (BMAT + DT*AMAT)*a(k+1) C = DT*b + BMAT*a(k) EACH TIME STEP DO 60 K=1,M B(K) = DT*B(K) DO 55 I=MAX(K-L,1),MIN(K+L,M) AMAT(K,I-K) = BMAT(K,I-K) + DT*AMAT(K,I-K) B(K) = B(K) + BMAT(K,I-K)*A(I) 55 CONTINUE 60 CONTINUE C CALL LBAND TO SOLVE BANDED LINEAR SYSTEM CALL LBAND(AMAT,A,B,M,L) PRINT 65, TKP1,(A(I),I=1,M) 65 FORMAT (6F10.5) 70 CONTINUE STOP END