C FIGURE 5.6.1 IMPLICIT DOUBLE PRECISION(A-H,O-Z) C NTRI = NUMBER OF TRIANGLES C M = NUMBER OF UNKNOWNS PARAMETER (NTRI=16, M=5) DIMENSION TRIX(3,NTRI),TRIY(3,NTRI),NODES(3,NTRI),B(M),A(M) C AMAT DIMENSIONED (M,-L:L), BUT HALF C BANDWIDTH L UNKNOWN UNTIL EXECUTION TIME ALLOCATABLE AMAT(:,:) 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 FUND(X,Y) = 1.0 FUNA(X,Y) = 0.0 FUNF(X,Y) = 1.0 FUNR(X,Y) = 0.0 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 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 ALLOCATE (AMAT(M,-L:L)) C ZERO BANDED COEFFICIENT MATRIX AMAT AND C RIGHT-HAND SIDE VECTOR B. DO 30 K=1,M B(K) = 0.0 DO 25 ISUB=-L,L AMAT(K,ISUB) = 0.0 25 CONTINUE 30 CONTINUE C BEGIN ELEMENT-BY-ELEMENT ASSEMBLY OF C MATRIX AND RIGHT-HAND SIDE VECTOR. 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,OMEGAX,OMEGAY REPRESENT THE C VALUE, X-DERIVATIVE AND Y-DERIVATIVE C OF OMEGA AT THE TRIANGLE MIDPOINT. OMEGA = 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)) OMEGA = OMEGA + R*PHI (J) OMEGAX = OMEGAX + R*PHIX(J) OMEGAY = OMEGAY + R*PHIY(J) ENDIF 35 CONTINUE C DMID,AMID,FMID ARE THE VALUES OF THE C COEFFICIENT FUNCTIONS D(X,Y),A(X,Y), C F(X,Y) AT THE TRIANGLE MIDPOINT. XMID = (X1+X2+X3)/3.0 YMID = (Y1+Y2+Y3)/3.0 DMID = FUND(XMID,YMID) AMID = FUNA(XMID,YMID) FMID = FUNF(XMID,YMID) AREA = ABS(DET) C CALCULATE CONTRIBUTIONS OF TRIANGLE KTRI C TO MATRIX AND RIGHT-HAND SIDE INTEGRALS, 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) + & (FMID*PHI(KX) & -DMID*(OMEGAX*PHIX(KX) + OMEGAY*PHIY(KX)) & -AMID*OMEGA*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 ENDIF 40 CONTINUE ENDIF 45 CONTINUE 50 CONTINUE C CALL LBAND TO SOLVE BANDED LINEAR SYSTEM CALL LBAND(AMAT,A,B,M,L) PRINT 55, (A(I),I=1,M) 55 FORMAT (5F10.5) STOP END